How to find out Geometric Median

Question

The question is:

Given N points(in 2D) with x and y coordinates, find a point P (in N given points) such that the sum of distances from other(N-1) points to P is minimum.

This point is commonly known as Geometric Median. Is there any efficient algorithm to solve this problem, other than the naive O(N^2) one?

Answer 1

As mentioned earlier, the type of algorithm to use depends on the way you measure distance. Since your question does not specify this measure, here are C implementations for both the Manhattan distance and the Squared Euclidean distance. Use dim = 2 for 2D points. Complexity O(n log n).

Manhattan distance

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);
}

Short explanation: We can sum the distance per dimension, 2 in your case. Say we have N points and the values in one dimension are v_0, .., v_(N-1) and T = v_0 + .. + v_(N-1). Then for each value v_i we have S_i = v_0 .. v_(i-1). Now we can express the Manhattan distance for this value by summing those on the left side: i * v_i - S_i and the right side: T - S_i - (N - i) * v_i, which results in (2 * i - N) * v_i - 2 * S_i + T. Adding T to all elements does not change the order, so we leave that out. And S_i can be computed on the fly.

Here is the rest of the code that makes it into an actual C program:

#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;
}

Squared Euclidean distance

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);
}

Shorter explanation: Pretty much the same approach as the previous, but with a slightly more complicated derivation. Say TT = v_0^2 + .. + v_(N-1)^2 we get TT + N * v_i^2 - 2 * v_i^2 * T. Again TT is added to all so it can be left out. More explanation on request.

Answer 2

You can solve the problem as a convex programming (The objective function is not always convex). The convex program can be solved using an iterative such as L-BFGS. The cost for each iteration is O(N) and usually the number of required iteration is not large. One important point to reduce the number of required iterations is that we know the optimum answer is one of the point in the input. So, the optimization can be stopped when its answer become close to one of the input points.

Answer 3

It appears that the problem is difficult to solve in better than O(n^2) time when using Euclidean distances. However the point that minimizes the sum of Manhattan distances to other points or the point that minimizes the sum of squares of Euclidean distances to other points can be found in O(n log n) time. (Assuming multiplying two numbers is O(1)). Let me shamelessly copy/paste my solution for Manhattan distances from a recent post:

Create a sorted array of x-coordinates and for each element in the array compute the "horizontal" cost of choosing that coordinate. The horizontal cost of an element is the sum of distances to all the points projected onto the X-axis. This can be computed in linear time by scanning the array twice (once from left to right and once in the reverse direction). Similarly create a sorted array of y-coordinates and for each element in the array compute the "vertical" cost of choosing that coordinate.

Now for each point in the original array, we can compute the total cost to all other points in O(1) time by adding the horizontal and vertical costs. So we can compute the optimal point in O(n). Thus the total running time is O(n log n).

We can follow a similar approach for computing the point that minimizes the sum of squares of Euclidean distances to other points. Let the sorted x-coordinates be: x₁, x₂, x₃, ..., x_n. We scan this list from left to right and for each point x_i we compute:

l_i = sum of distances to all the elements to the left of x_i = (x_i-x₁) + (x_i-x₂) + .... + (x_i-x_i-1) , and

sl_i = sum of squares of distances to all the elements to the left of x_i = (x_i-x₁)^2 + (x_i-x₂)^2 + .... + (x_i-x_i-1)^2

Note that given l_i and sl_i we can compute l_i+1 and sl_i+1 in O(1) time as follows:

Let d = x_i+1-x_i. Then:

l_i+1 = l_i + id and sl_i+1 = sl_i + id^2 + 2*i*d

Thus we can compute all the l_i and sl_i in linear time by scanning from left to right. Similarly, for every element we can compute the r_i: sum of distances to all elements to the right and the sr_i: sum of squares of distances to all the elements to the right in linear time. Adding sr_i and sl_i for each i, gives the sum of squares of horizontal distances to all the elements, in linear time. Similarly, compute the sum of squares of vertical distances to all the elements.

Then we can scan through the original points array and find the point that minimizes the sum of squares of vertical and horizontal distances as before.

Answer 4

I implemented the Weiszfeld method (I know it's not what you are looking for, but it may help to aproximate your Point), the complexity is O(N*M/k) where N is the number of points, M the dimension of the points (in your case is 2) and k is the error desired:

https://github.com/j05u3/weiszfeld-implementation

Answer 5

I solved something similar for a local online judge once using simulated annealing. That was the official solution as well and the program got AC.

The only difference was that the point I had to find did not have to be part of the N given points.

This was my C++ code, and N could be as large as 50000. The program executes in 0.1s on a 2ghz pentium 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;
}

Then I think It's correct to pick the one from your list that is closest to the (x, y) returned by this algorithm.

This algorithm takes advantage of what this wikipedia paragraph on the geometric median says:

However, it is straightforward to calculate an approximation to the geometric median using an iterative procedure in which each step produces a more accurate approximation. Procedures of this type can be derived from the fact that the sum of distances to the sample points is a convex function, since the distance to each sample point is convex and the sum of convex functions remains convex. Therefore, procedures that decrease the sum of distances at each step cannot get trapped in a local optimum.

One common approach of this type, called Weiszfeld's algorithm after the work of Endre Weiszfeld,[4] is a form of iteratively re-weighted least squares. This algorithm defines a set of weights that are inversely proportional to the distances from the current estimate to the samples, and creates a new estimate that is the weighted average of the samples according to these weights. That is,

The first paragraph above explains why this works: because the function we are trying to optimize does not have any local minimums, so you can greedily find the minimum by iteratively improving it.

Think of this as a sort of binary search. First, you approximate the result. A good approximation will be the center of gravity, which my code computes when reading the input. Then, you see if adjacent points to this give you a better solution. In this case, a point is considered adjacent if it as a distance of step away from your current point. If it is better, then it is fine to discard your current point, because, as I said, this will not trap you into a local minimum because of the nature of the function you are trying to minimize.

After this, you half the step size, just like in binary search, and continue until you have what you consider to be a good enough approximation (controlled by the eps constant).

The complexity of the algorithm therefore depends on how accurate you want the result to be.

Answer 6

Step 1: Sort the points collection by x-dimension (nlogn)
Step 2: Calculate the x-distance between each point and all points TO THE LEFT of it:

xLDist[0] := 0
for i := 1 to n - 1
       xLDist[i] := xLDist[i-1] + ( ( p[i].x - p[i-1].x ) * i)

Step 3: Calculate the x-distance between each point and all points TO THE RIGHT of it:

xRDist[n - 1] := 0
for i := n - 2 to 0
       xRDist[i] := xRDist[i+1] + ( ( p[i+1].x - p[i].x ) * i)  

Step 4: Sum both up you'll get the total x-distance from each point to the other N-1 points

for i := 0 to n - 1
       p[i].xDist = xLDist[i] + xRDist[i]

Repeat Step 1,2,3,4 with the y-dimension to get p[i].yDist

The point with the smallest sum of xDist and yDist is the answer

Total Complexity O(nlogn)

Answer in C++

Further explanation:
The idea is to reuse the already computed total distance of the previous point.
Lets say we have 3 point ABCD sorted, we see that the total left distance of D to the others before it are:

AD + BD + CD = (AC + CD) + (BC + CD) + CD = AC + BC + 3CD

In which (AC + BC) is the total left distance of C to the others before it, we took advantage of this and only need to compute ldist(C) + 3CD