问题是:
给定具有x和y坐标的N个点(在2D中),找到一个点P(在N个给定点中),以使从其他(N-1)个点到P的距离之和最小。
这一点通常称为“ 几何中值”。除了幼稚的算法以外,是否有任何有效的算法可以解决此问题O(N^2)?
O(N^2)
我曾经使用模拟退火为本地在线法官解决了类似问题。那也是官方的解决方案,程序获得了AC。
唯一的区别是,我必须找到的点不必一定是N给定点的一部分。
N
这是我的C ++代码,N可能和一样大50000。该程序在0.1s2GHz的奔腾4上执行。
50000
0.1s
// 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)此算法返回值的列表是正确的。
(x, y)
该算法利用了维基百科上有关几何中位数的内容:
但是,使用迭代过程来计算几何中值的近似值很简单,在该过程中,每个步骤都会产生更精确的近似值。这种类型的过程可以从以下事实得出:到采样点的距离之和是凸函数,因为到每个采样点的距离是凸的而凸函数之和保持凸。因此,减少每一步距离之和的过程不能陷入局部最优状态。 在Endre Weiszfeld [4]的工作之后,这种类型的常见方法称为Weiszfeld算法[4],是迭代地重新加权最小二乘的一种形式。该算法定义了一组权重,这些权重与从当前估计值到样本的距离成反比,并根据这些权重创建一个新的估计值,即样本的加权平均值。那是,
但是,使用迭代过程来计算几何中值的近似值很简单,在该过程中,每个步骤都会产生更精确的近似值。这种类型的过程可以从以下事实得出:到采样点的距离之和是凸函数,因为到每个采样点的距离是凸的而凸函数之和保持凸。因此,减少每一步距离之和的过程不能陷入局部最优状态。
在Endre Weiszfeld [4]的工作之后,这种类型的常见方法称为Weiszfeld算法[4],是迭代地重新加权最小二乘的一种形式。该算法定义了一组权重,这些权重与从当前估计值到样本的距离成反比,并根据这些权重创建一个新的估计值,即样本的加权平均值。那是,
上面的第一段说明了它的工作原理:由于我们要优化的函数没有任何局部最小值,因此您可以通过迭代地改进它来贪婪地找到最小值。
将此视为一种二进制搜索。首先,您估算结果。重心是一个很好的近似值,我的代码在读取输入时会计算重心。然后,您可以查看与此相邻的点是否为您提供更好的解决方案。在这种情况下,如果某点距step您当前点的距离为零,则该点被视为相邻点。如果更好,那么最好放弃当前点,因为正如我所说,由于您要最小化的函数的性质,这不会使您陷入局部最小值。
step
此后,将步长减小一半,就像在二进制搜索中一样,并继续进行操作,直到获得您认为足够好的近似值(由eps常数控制)为止。
eps
因此,算法的复杂度取决于您希望结果的准确性。