Возникает вопрос:
Учитывая N точек (в 2D) с координатами x и y, найдите точку P (в N заданные точки), так что сумма расстояний от других (N-1) точек P минимально.
Эта точка широко известна как Геометрическая медиана. Есть ли эффективный алгоритм для решения этой проблемы, кроме наивного O(N^2) один?
Я решил что-то подобное для локального онлайн-судьи, используя имитированный отжиг. Это было официальное решение, и программа получила AC.
Единственное различие заключалось в том, что точка, которую я должен был найти, не должна была быть частью N заданных точек.
Это мой код на С++, а N может быть размером 50000. Программа выполняется в 0.1s на 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;
}
Затем я думаю, что правильно выбрать тот из вашего списка, который ближе всего к (x, y), возвращаемому этим алгоритмом.
Этот алгоритм использует преимущества этого параграфа википедии в геометрической медианной области:
Однако, нетрудно вычислить приближение к геометрическая медиана с использованием итерационной процедуры, в которой каждый шаг дает более точное приближение. Процедуры такого типа могут быть вытекает из того факта, что сумма расстояний до точек выборки является выпуклой функцией, так как расстояние до каждой точки выборки равно выпуклой, а сумма выпуклых функций остается выпуклой. Следовательно, процедуры, которые уменьшают сумму расстояний на каждом шаге, не могут в ловушке локального оптимума.
Один общий подход этого типа, называемый Вайзфельда после работы Эндре Вайсфельда [4] является формой итерационно повторно взвешенных наименьших квадратов. Этот алгоритм определяет набор весов, которые обратно пропорциональны расстояниям от текущей оценки к образцам, и создает новую оценку, которая взвешенное среднее количество образцов в соответствии с этими весами. Что есть,
В первом абзаце выше объясняется, почему это работает: потому что функция, которую мы пытаемся оптимизировать, не имеет локальных минимумов, поэтому вы можете с жадностью находить минимум, итеративно улучшая ее.
Подумайте об этом как о бинарном поиске. Во-первых, вы приближаете результат. Хорошим приближением будет центр тяжести, который мой код вычисляет при чтении ввода. Затем вы видите, что соседние точки дают вам лучшее решение. В этом случае точка считается смежной, если она находится на расстоянии step от вашей текущей точки. Если это будет лучше, тогда хорошо удалить эту текущую точку, потому что, как я уже сказал, это не приведет вас к локальному минимуму из-за характера функции, которую вы пытаетесь свести к минимуму.
После этого вы получите половину размера шага, как в бинарном поиске, и продолжайте, пока не получите то, что считаете достаточно хорошим приближением (контролируется константой eps).
Таким образом, сложность алгоритма зависит от того, насколько точным вы хотите получить результат.
Похоже, что проблему трудно решить более чем за время 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 - х <суб > I-1суб > ) ^ 2
Обратите внимание, что при l i и sl i мы можем вычислить l я + 1 и sl я + 1 в O(1) времени следующим образом:
Пусть d = x я + 1 -x i. Тогда:
l я + 1= l i + id и sl я + 1= sl i + id ^ 2 + 2 * я * d
Таким образом, мы можем вычислить все l i и sl i в линейном времени путем сканирования слева направо. Аналогично, для каждого элемента мы можем вычислить r i: сумма расстояний до всех элементов справа и sr i: сумма квадратов расстояний ко всем элементам справа в линейном время. Добавление sr i и sl i для каждого я дает сумму квадратов горизонтальных расстояний ко всем элементам в линейном времени. По аналогии, вычислить сумму квадратов вертикальных расстояний ко всем элементам.
Затем мы можем сканировать исходный массив точек и найти точку, которая минимизирует сумму квадратов вертикальных и горизонтальных расстояний, как и раньше.
Как упоминалось ранее, тип используемого алгоритма зависит от способа измерения расстояния. Поскольку ваш вопрос не определяет эту меру, здесь представлены реализации 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 добавляется ко всем, поэтому его можно оставить без внимания. Дополнительные пояснения по запросу.
Я реализовал метод Weiszfeld (я знаю, что это не то, что вы ищете, но это может помочь aproximate вашей Point), сложность O (N * M/k), где N - количество точек, M размерность точек (в вашем случае 2), а k - желаемая ошибка:
Шаг 1: Сортировка коллекции точек по x-dimension (nlogn)
Шаг 2: Рассчитайте расстояние между каждой точкой и всеми точками ВЛЕВО:
xLDist[0] := 0
for i := 1 to n - 1
xLDist[i] := xLDist[i-1] + ( ( p[i].x - p[i-1].x ) * i)
Шаг 3: Рассчитайте расстояние между каждой точкой и всеми точками ВПРАВО:
xRDist[n - 1] := 0
for i := n - 2 to 0
xRDist[i] := xRDist[i+1] + ( ( p[i+1].x - p[i].x ) * i)
Шаг 4: Суммируйте оба значения: вы получите общее расстояние по x от каждой точки до других точек N-1.
for i := 0 to n - 1
p[i].xDist = xLDist[i] + xRDist[i]
Повторите шаг 1,2,3,4 с y-размерностью, чтобы получить
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
Вы можете решить проблему как выпуклое программирование (объектная функция не всегда выпукла). Выпуклая программа может быть решена с использованием итерации, такой как L-BFGS. Стоимость каждой итерации равна O (N), и обычно количество требуемой итерации невелико. Одним из важных моментов для сокращения числа требуемых итераций является то, что мы знаем, что оптимальный ответ является одной из точек ввода. Таким образом, оптимизация может быть остановлена, когда ее ответ приблизится к одной из входных точек.