如何找出几何中位数

Question

问题是:

给定具有x和y坐标的N个点(在2D中),找到点P(在N个给定点中),使得从其他(N-1)个点到P的距离之和最小.

这一点通常称为几何中位数.有没有有效的算法来解决这个问题,除了天真的算法O(N^2)？

Answer 1

一旦使用模拟退火,我为当地的在线评委解决了类似问题.这也是官方解决方案,该计划得到了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的工作之后称为Weiszfeld的算法,[4]是迭代重新加权最小二乘的一种形式.该算法定义一组权重,其与从当前估计到样本的距离成反比,并根据这些权重创建新的估计,即样本的加权平均.那是,

上面的第一段解释了为什么这样做:因为我们试图优化的函数没有任何局部最小值,所以你可以通过迭代地改进它来贪婪地找到最小值.

可以把它想象成一种二元搜索.首先,您估算结果.一个很好的近似将是重心,我的代码在读取输入时计算.然后,您会看到相邻的点是否为您提供了更好的解决方案.在这种情况下,如果一个点距离step当前点的距离,则认为该点是相邻的.如果它更好,那么丢弃你当前的点是好的,因为,正如我所说,由于你试图最小化的功能的性质,这不会陷入局部最小值.

在此之后,您将步长的一半,就像二进制搜索一样,并继续直到您认为是足够好的近似值(由eps常量控制).

因此,算法的复杂性取决于您希望结果的准确程度.

Answer 2

O(n^2)在使用欧几里德距离时,似乎难以解决问题.然而,可以及时找到最小化曼哈顿到其他点的距离的总和或最小化欧几里德到其他点的距离的平方的点的点O(n log n).(假设乘以两个数字O(1)).让我无耻地从最近的帖子中复制/粘贴我的曼哈顿距离解决方案:

创建一个排序的x坐标数组,并为数组中的每个元素计算选择该坐标的"水平"成本.元素的水平成本是投影到X轴上的所有点的距离之和.这可以通过扫描阵列两次(一次从左到右,一次在反方向)以线性时间计算.类似地,创建一个y坐标的排序数组,并为数组中的每个元素计算选择该坐标的"垂直"成本.

现在,对于原始数组中的每个点,我们可以通过添加水平和垂直成本来计算O(1)时间内所有其他点的总成本.因此我们可以计算O(n)中的最佳点.因此总运行时间为O(n log n).

我们可以遵循类似的方法来计算最小化欧几里德到其他点的距离的平方和的点.令排序的x坐标为:x ₁,x ₂,x ₃,...,x _n.我们从左到右扫描这个列表,对于我们计算的每个点x _i:

l _i =到x _i =(x _i -x ₁)+(x _i -x ₂)+ .... +(x _i -x _i-1)左边所有元素的距离之和,以及

sl _i =到x左边所有元素的距离的平方和_{i i} =(x _i -x ₁)^ 2 +(x _i -x ₂)^ 2 + .... +(x _i -x _{i -1})^ 2

需要注意的是动量l _我和SL _我,我们可以计算升_{I + 1}和SL _{I + 1}的O(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:到右边的所有元素的距离之和以及sr _i:线性时间内到右边所有元素的距离的平方和.为每个i 添加sr _i和sl _i,以线性时间给出所有元素的水平距离的平方和.同样,计算所有元素的垂直距离的平方和.

然后我们可以扫描原始点阵列并找到最小化垂直和水平距离的平方和的点,如前所述.

Answer 3

如前所述,要使用的算法类型取决于您测量距离的方式.由于您的问题未指定此度量,因此这里是曼哈顿距离和平方欧几里德距离的 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添加到所有,所以它可以被省略.更多解释请求.