问题是：
给定具有 x 和 y 坐标的 N 个点（在 2D 中），找到一个点 P（在 N 个给定点中），使得从其他（N-1）个点到 P 的距离之和最小。
这个点通常被称为几何中位数。除了天真的算法之外，还有什么有效的算法可以解决这个问题O(N^2)吗？
问问题
32227 次
7 回答
24
我曾经使用模拟退火为当地在线法官解决了类似的问题。这也是官方的解决方案，程序得到了 AC。
唯一的区别是我必须找到的点不必是N给定点的一部分。
这是我的 C++ 代码，N可能与50000. 该程序在0.1s2ghz 奔腾 4 上执行。
// header files for IO functions and math
#include <cstdio>
#include <cmath>
// the maximul value n can take
const int maxn = 50001;
// given a point (x, y) on a grid, we can find its left/right/up/down neighbors
// by using these constants: (x + dx[0], y + dy[0]) = upper neighbor etc.
const int dx[] = {-1, 0, 1, 0};
const int dy[] = {0, 1, 0, -1};
// controls the precision - this should give you an answer accurate to 3 decimals
const double eps = 0.001;
// input and output files
FILE *in = fopen("adapost2.in","r"), *out = fopen("adapost2.out","w");
// stores a point in 2d space
struct punct
{
double x, y;
};
// how many points are in the input file
int n;
// stores the points in the input file
punct a[maxn];
// stores the answer to the question
double x, y;
// finds the sum of (euclidean) distances from each input point to (x, y)
double dist(double x, double y)
{
double ret = 0;
for ( int i = 1; i <= n; ++i )
{
double dx = a[i].x - x;
double dy = a[i].y - y;
ret += sqrt(dx*dx + dy*dy); // classical distance formula
}
return ret;
}
// reads the input
void read()
{
fscanf(in, "%d", &n); // read n from the first
// read n points next, one on each line
for ( int i = 1; i <= n; ++i )
fscanf(in, "%lf %lf", &a[i].x, &a[i].y), // reads a point
x += a[i].x,
y += a[i].y; // we add the x and y at first, because we will start by approximating the answer as the center of gravity
// divide by the number of points (n) to get the center of gravity
x /= n;
y /= n;
}
// implements the solving algorithm
void go()
{
// start by finding the sum of distances to the center of gravity
double d = dist(x, y);
// our step value, chosen by experimentation
double step = 100.0;
// done is used to keep track of updates: if none of the neighbors of the current
// point that are *step* steps away improve the solution, then *step* is too big
// and we need to look closer to the current point, so we must half *step*.
int done = 0;
// while we still need a more precise answer
while ( step > eps )
{
done = 0;
for ( int i = 0; i < 4; ++i )
{
// check the neighbors in all 4 directions.
double nx = (double)x + step*dx[i];
double ny = (double)y + step*dy[i];
// find the sum of distances to each neighbor
double t = dist(nx, ny);
// if a neighbor offers a better sum of distances
if ( t < d )
{
update the current minimum
d = t;
x = nx;
y = ny;
// an improvement has been made, so
// don't half step in the next iteration, because we might need
// to jump the same amount again
done = 1;
break;
}
}
// half the step size, because no update has been made, so we might have
// jumped too much, and now we need to head back some.
if ( !done )
step /= 2;
}
}
int main()
{
read();
go();
// print the answer with 4 decimal points
fprintf(out, "%.4lf %.4lf\n", x, y);
return 0;
}
然后我认为从您的列表中选择最接近(x, y)该算法返回的一个是正确的。
这个算法利用了这个关于几何中位数的维基百科段落所说的：
然而，使用迭代过程计算几何中值的近似值很简单，其中每一步都会产生更准确的近似值。这种类型的过程可以从到样本点的距离之和是一个凸函数这一事实推导出来，因为到每个样本点的距离是凸的并且凸函数的总和仍然是凸的。因此，减少每一步距离总和的程序不能陷入局部最优。
这种类型的一种常见方法，在 Endre Weiszfeld [4] 的工作之后称为 Weiszfeld 算法，是一种迭代重新加权最小二乘的形式。该算法定义了一组与当前估计值到样本的距离成反比的权重，并根据这些权重创建一个新的估计值，即样本的加权平均值。那是，
上面的第一段解释了为什么会这样：因为我们试图优化的函数没有任何局部最小值，所以你可以通过迭代改进来贪婪地找到最小值。
将其视为一种二分搜索。首先，你近似结果。一个很好的近似值是我的代码在读取输入时计算的重心。然后，您会查看与此相邻的点是否为您提供了更好的解决方案。step在这种情况下，如果一个点与您的当前点相距一定距离，则该点被认为是相邻的。如果它更好，那么丢弃您当前的点就可以了，因为正如我所说，由于您试图最小化的函数的性质，这不会使您陷入局部最小值。
在此之后，步长减半，就像二分搜索一样，并继续，直到你得到你认为足够好的近似值（由eps常数控制）。
因此，算法的复杂性取决于您希望结果有多准确。
于 2012-10-17T12:38:03.027 回答
10
O(n^2)使用欧几里得距离时，似乎比时间更难解决问题。但是，可以及时找到使曼哈顿到其他点的距离之和最小化的点或使到其他点的欧几里得距离的平方和O(n log n)最小化的点。（假设两个数字相乘是O(1)）。让我从最近的一篇文章中无耻地复制/粘贴我的曼哈顿距离解决方案：
创建一个排序的 x 坐标数组，并为数组中的每个元素计算选择该坐标的“水平”成本。元素的水平成本是到投影到 X 轴上的所有点的距离之和。这可以通过扫描阵列两次（一次从左到右，一次在相反方向）以线性时间计算。类似地，创建一个 y 坐标的排序数组，并为数组中的每个元素计算选择该坐标的“垂直”成本。
现在对于原始数组中的每个点，我们可以通过添加水平和垂直成本来计算 O(1) 时间内所有其他点的总成本。所以我们可以计算 O(n) 中的最佳点。因此总运行时间为 O(n log n)。
我们可以采用类似的方法来计算使欧几里得距离平方和最小化到其他点的点。让排序后的 x 坐标为： x 1 , x 2 , x 3 , ..., x n。我们从左到右扫描这个列表，并为每个点 x i我们计算：
l i = 到 x i = (x i -x 1 ) + (x i -x 2 ) + .... + (x i -x i-1 )左侧所有元素的距离总和，并且
sl i = 到 x i左侧所有元素的距离平方和= (x i -x 1 )^2 + (x i -x 2 )^2 + .... + (x i -x i -1 )^2
请注意，给定 l i和 sl i我们可以及时计算 l i+1和sl i+1O(1)如下：
令 d = x i+1 -x i。然后：
l i+1 = l i + i d 和 sl i+1 = sl i + i d^2 + 2*i*d
因此，我们可以通过从左到右扫描，在线性时间内计算所有的 l i和 sl i 。类似地，对于每个元素，我们可以计算 r i：到右侧所有元素的距离之和，以及 sri ：在线性时间内到右侧所有元素的距离平方和。为每个 i添加 sr i和 sl i，以线性时间给出到所有元素的水平距离平方和。同样，计算到所有元素的垂直距离的平方和。
然后我们可以像以前一样扫描原始点数组，找到使垂直和水平距离的平方和最小的点。
于 2012-10-17T14:39:31.177 回答
5
如前所述，要使用的算法类型取决于您测量距离的方式。由于您的问题未指定此度量，因此这里是曼哈顿距离和平方欧几里得距离的 C 实现。用于dim = 2二维点。复杂性O(n log n)。
曼哈顿距离
double * geometric_median_with_manhattan(double **points, int N, int dim) {
for (d = 0; d < dim; d++) {
qsort(points, N, sizeof(double *), compare);
double S = 0;
for (int i = 0; i < N; i++) {
double v = points[i][d];
points[i][dim] += (2 * i - N) * v - 2 * S;
S += v;
}
}
return min(points, N, dim);
}
简短说明：我们可以将每个维度的距离相加，在您的情况下为 2。假设我们有点N，一维中的值是v_0、 ..v_(N-1)和T = v_0 + .. + v_(N-1)。v_i然后对于我们拥有的每个值S_i = v_0 .. v_(i-1)。现在我们可以通过将左侧的值相加来表示该值的曼哈顿距离：i * v_i - S_i和右侧：T - S_i - (N - i) * v_i，结果为(2 * i - N) * v_i - 2 * S_i + T。添加T到所有元素不会更改顺序，因此我们将其省略。并且S_i可以即时计算。
以下是使其成为实际 C 程序的其余代码：
#include <stdio.h>
#include <stdlib.h>
int d = 0;
int compare(const void *a, const void *b) {
return (*(double **)a)[d] - (*(double **)b)[d];
}
double * min(double **points, int N, int dim) {
double *min = points[0];
for (int i = 0; i < N; i++) {
if (min[dim] > points[i][dim]) {
min = points[i];
}
}
return min;
}
int main(int argc, const char * argv[])
{
// example 2D coordinates with an additional 0 value
double a[][3] = {{1.0, 1.0, 0.0}, {3.0, 1.0, 0.0}, {3.0, 2.0, 0.0}, {0.0, 5.0, 0.0}};
double *b[] = {a[0], a[1], a[2], a[3]};
double *min = geometric_median_with_manhattan(b, 4, 2);
printf("geometric median at {%.1f, %.1f}\n", min[0], min[1]);
return 0;
}
平方欧几里得距离
double * geometric_median_with_square(double **points, int N, int dim) {
for (d = 0; d < dim; d++) {
qsort(points, N, sizeof(double *), compare);
double T = 0;
for (int i = 0; i < N; i++) {
T += points[i][d];
}
for (int i = 0; i < N; i++) {
double v = points[i][d];
points[i][dim] += v * (N * v - 2 * T);
}
}
return min(points, N, dim);
}
更简短的解释：与前面的方法几乎相同，但推导稍微复杂一些。说TT = v_0^2 + .. + v_(N-1)^2我们得到TT + N * v_i^2 - 2 * v_i^2 * T. 再次将 TT 添加到所有内容中，因此可以将其排除在外。应要求提供更多解释。
于 2013-11-19T09:41:24.200 回答
2
我实现了 Weiszfeld 方法（我知道这不是你要找的，但它可能有助于近似你的点），复杂度是 O(N*M/k) 其中 N 是点的数量，M 的维度点（在您的情况下为 2），并且 k 是所需的错误：
于 2015-09-03T17:36:36.690 回答
2
第 1 步：按 x 维度 (nlogn) 对点集合进行排序
第 2 步：计算每个点与其左侧的所有点之间的 x距离：
xLDist[0] := 0
for i := 1 to n - 1
xLDist[i] := xLDist[i-1] + ( ( p[i].x - p[i-1].x ) * i)
第 3 步：计算每个点与其右侧的所有点之间的 x距离：
xRDist[n - 1] := 0
for i := n - 2 to 0
xRDist[i] := xRDist[i+1] + ( ( p[i+1].x - p[i].x ) * i)
第 4 步：将两者相加，您将得到从每个点到其他 N-1 个点的总 x 距离
for i := 0 to n - 1
p[i].xDist = xLDist[i] + xRDist[i]
使用 y 维度重复步骤 1、2、3、4 以获得
p[i].yDist
xDist和的和最小的点yDist就是答案
总复杂度 O(nlogn)
进一步解释：
这个想法是重用已经计算的前一点的总距离。
假设我们对 3 点 ABCD 进行了排序，我们看到 D 与之前其他点的总左距离为：
AD + BD + CD = (AC + CD) + (BC + CD) + CD = AC + BC + 3CD
其中(AC + BC)是 C 与之前其他人的总左距离，我们利用了这一点，只需要计算ldist(C) + 3CD
于 2016-04-05T01:23:48.403 回答
0
您可以将问题解决为凸编程（目标函数并不总是凸的）。凸程序可以使用诸如 L-BFGS 之类的迭代来求解。每次迭代的成本为 O(N)，通常所需的迭代次数并不多。减少所需迭代次数的重要一点是，我们知道最佳答案是输入中的点之一。因此，当其答案接近输入点之一时，可以停止优化。
于 2012-10-19T20:51:40.660 回答
-1
我们需要找到的答案是几何中位数
C++ 中的代码
#include <bits/stdc++.h>
using namespace std;
int main()
{
int n;
cin >> n;
int a[n],b[n];
for(int i=0;i<n;i++)
cin >> a[i] >> b[i];
int res = 0;
sort(a,a+n);
sort(b,b+n);
int m1 = a[n/2];
int m2 = b[n/2];
for(int i=0;i<n;i++)
res += abs(m1 - a[i]);
for(int i=0;i<n;i++)
res += abs(m2 - b[i]);
cout << res << '\n';
}
于 2021-01-27T18:18:23.980 回答