On finding geometric median with ternary searches

#	User	Rating
1	jiangly	3810
2	Benq	3676
3	Kevin114514	3655
4	maroonrk	3463
5	strapple	3447
6	Um_nik	3387
7	heuristica	3322
8	turmax	3317
9	tourist	3307
10	jiangbowen	3291

#	User	Contrib.
1	Qingyu	156
2	nik_exists	150
2	maspy	150
4	Um_nik	141
5	Errichto	139
6	adamant	137
7	AmShZ	136
8	BledDest	132
9	maroonrk	131
10	qwexd	129

Hello community! Please consider the following two problems.

Given $$$N$$$ points on 2D plane, we need to calculate:

A) Coordinates of a new point such that the sum of Euclidean distances from that point to all of given points is minimal (Geometric median finding problem);

B) The index of one of initial points such that the sum of Euclidean distances from that point to all other given points is minimal (Medoid finding problem).

My first hypothesis claims that problem (A) could be solved in $$$O(N\log^2(C))$$$ time (and $$$C$$$ is the absolute value of the maximal coordinate of point) with two nested ternary searches. It is so because of Euclidean distance between given point and some other point is a convex function, and the sum of convex functions is still a convex function, so we can find the global extremum point with ternary search.

I believe I saw that idea a time or two as a suggested solution for Fermat point problem, but I still can't find any mention of this algorithm in papers considering problem (A). So I have some doubts that it may give wrong answer, but I can't find a counterexample.

For the problem (B), my idea is that we can use the solution of problem (A) and then find the point from given array that is closest to the calculated one (or any of such points if there are multiple). But, again, I could not find any mention of this approach in papers.

So I tried to make a stress-test via Polygon to find a counterexample. I use the code that implements an approach described above and a bruteforce code that juct checks all possible points in $$$O(N^2)$$$ time.

Ternary searches code

vector<pair<double, double>> p;

double dist(double x1, double y1, double x2, double y2) {
    return hypot(x1 - x2, y1 - y2);
}

double dist(double x, double y) {
    double res = 0;
    for (int i = 0; i < p.size(); i++)
        res += dist(x, y, p[i].first, p[i].second);
    return res;
}

pair<double, double> findY(double x) {
    double yl = -1e9, yr = 1e9;
    for (int i = 0; i < 60; i++) {
        double ym1 = yl + (yr - yl) / 3;
        double ym2 = yr - (yr - yl) / 3;
        double d1 = dist(x, ym1);
        double d2 = dist(x, ym2);
        if (d1 < d2)
            yr = ym2;
        else
            yl = ym1;
    }
    return { yl, dist(x, yl) };
}

pair<double, double> findXY() {
    double xl = -1e9, xr = 1e9;
    for (int i = 0; i < 60; i++) {
        double xm1 = xl + (xr - xl) / 3;
        double xm2 = xr - (xr - xl) / 3;
        auto [y1, d1] = findY(xm1);
        auto [y2, d2] = findY(xm2);
        if (d1 < d2)
            xr = xm2;
        else
            xl = xm1;
    }
    return { xl, findY(xl).first };
}

int main() {
    int n;
    scanf("%d", &n);

    p.resize(n);
    for (int i = 0; i < n; i++)
        scanf("%lf%lf", &p[i].first, &p[i].second);

    auto [x, y] = findXY();

    int bestI = 0;
    double bestD = dist(x, y, p[0].first, p[0].second);
    for (int i = 1; i < n; i++) {
        double d = dist(x, y, p[i].first, p[i].second);
        if (d < bestD) {
            bestD = d;
            bestI = i;
        }
    }

    printf("%.10lf", dist(p[bestI].first, p[bestI].second));
}

Brute force code

vector<pair<double, double>> p;

double dist(double x1, double y1, double x2, double y2) {
    return hypot(x1 - x2, y1 - y2);
}

double dist(double x, double y) {
    double res = 0;
    for (int i = 0; i < p.size(); i++)
        res += dist(x, y, p[i].first, p[i].second);
    return res;
}

int main() {
    int n;
    scanf("%d", &n);

    p.resize(n);
    for (int i = 0; i < n; i++)
        scanf("%lf%lf", &p[i].first, &p[i].second);


    int bestI = 0;
    double bestD = dist(p[0].first, p[0].second);
    for (int i = 1; i < n; i++) {
        double d = dist(p[i].first, p[i].second);
        if (d < bestD) {
            bestD = d;
            bestI = i;
        }
    }

    printf("%.10lf", bestD);
}

Unfortunately stress-testing gave nothing, so I come with two questions for geometry-skilled colleagues:

Is my approach for problem (A) right? (I believe it is)
Is my approach for problem (B) right? (I believe it is not). How may the counterexample look like?

Thanks in advance!

Comments (4)

Write comment?

Um_nik

7 years ago, hide # |

Yes, and you proved it. The reason you cannot find any papers about it is probably because nested ternary search is not a go-to method for optimization problems on convex functions.
I'm too lazy to check, but basically you need to construct some test where near optimum you have different derivatives along different directions, then repeat it many times so that effect of added points (real optimum and what your algo will find) will be negligible. This should work:

repeated 1000 times {
-1000000 -1
-1000000 1
1000000 -1
1000000 1
}
2 0
-2 0
0 1
0 -1

→ Reply

shortHorse

7 years ago, hide # ^ |

Oh, thank you, you are right. When we change 2 and 1 in the last 4 lines to 10 and 5 respectively, the error increases so that it becomes noticeable by checker.

Thank you again for your help!

insert_cool_handle

6 years ago, hide # |

You seem to be using the following theorem(?):

If $$$\displaystyle f:R^n \rightarrow R$$$ is a convex function, then the function $$$g(x_2, x_3, \ldots, x_n) = \displaystyle \min_{x_1} f(x_1, x_2, \ldots x_n)$$$ is also convex.

I was unable to prove/disprove this. Is it well known?

SYury

6 years ago, hide # ^ |

The proof can be found in this book, pages 87-88

shortHorse's blog