algorithm - promedio - Cómo averiguar la mediana geométrica
media geometrica para datos agrupados (6)
La pregunta es:
Dados los puntos N (en 2D) con las coordenadas x e y, encuentre un punto P (en N puntos dados) tal que la suma de las distancias de otros puntos (N-1) a P sea mínima.
Este punto es comúnmente conocido como Mediana geométrica . ¿Hay algún algoritmo eficiente para resolver este problema, aparte del ingenuo O(N^2) ?
Como se mencionó anteriormente, el tipo de algoritmo a utilizar depende de la forma en que mida la distancia. Como su pregunta no especifica esta medida, aquí hay implementaciones en C tanto para la distancia de Manhattan como para la distancia euclidiana al cuadrado . Use dim = 2 para los puntos 2D. Complejidad O(n log n) .
Distancia de Manhattan
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);
}
Breve explicación: Podemos sumar la distancia por dimensión, 2 en su caso. Digamos que tenemos N puntos y los valores en una dimensión son v_0 , .., v_(N-1) y T = v_0 + .. + v_(N-1) . Luego, para cada valor v_i tenemos S_i = v_0 .. v_(i-1) . Ahora podemos expresar la distancia de Manhattan para este valor sumando las del lado izquierdo: i * v_i - S_i y el lado derecho: T - S_i - (N - i) * v_i , lo que resulta en (2 * i - N) * v_i - 2 * S_i + T Agregar T a todos los elementos no cambia el orden, así que lo dejamos fuera. Y S_i se puede calcular sobre la marcha.
Aquí está el resto del código que lo convierte en un programa de C real:
#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;
}
Distancia euclidiana al cuadrado
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);
}
Explicación más breve: Más o menos el mismo enfoque que el anterior, pero con una derivación un poco más complicada. Diga TT = v_0^2 + .. + v_(N-1)^2 obtenemos TT + N * v_i^2 - 2 * v_i^2 * T De nuevo, TT se agrega a todos para que se pueda omitir. Más explicación a petición.
Implementé el método de Weiszfeld (sé que no es lo que está buscando, pero puede ayudar a aproximar su Punto), la complejidad es O (N * M / k) donde N es el número de puntos, M la dimensión de la puntos (en su caso es 2) y k es el error deseado:
Parece que el problema es difícil de resolver en un tiempo mejor que O(n^2) cuando se usan distancias euclidianas. Sin embargo, el punto que minimiza la suma de las distancias de Manhattan a otros puntos o el punto que minimiza la suma de los cuadrados de las distancias euclidianas a otros puntos se puede encontrar en tiempo O(n log n) . (Suponiendo que multiplicar dos números es O(1) ). Permítame copiar y pegar descaradamente mi solución para las distancias de Manhattan desde una post reciente:
Cree una matriz ordenada de coordenadas x y, para cada elemento de la matriz, calcule el costo "horizontal" de elegir esa coordenada. El costo horizontal de un elemento es la suma de las distancias a todos los puntos proyectados en el eje X. Esto se puede calcular en tiempo lineal al escanear la matriz dos veces (una vez de izquierda a derecha y otra en la dirección inversa). De manera similar, cree una matriz ordenada de coordenadas y, y para cada elemento de la matriz, calcule el costo "vertical" de elegir esa coordenada.
Ahora, para cada punto de la matriz original, podemos calcular el costo total de todos los demás puntos en O (1) al sumar los costos horizontales y verticales. Entonces podemos calcular el punto óptimo en O (n). Por lo tanto, el tiempo total de ejecución es O (n log n).
Podemos seguir un enfoque similar para calcular el punto que minimiza la suma de los cuadrados de las distancias euclidianas a otros puntos. Sean las coordenadas x ordenadas: x 1 , x 2 , x 3 , ..., x n . Escaneamos esta lista de izquierda a derecha y para cada punto x i calculamos:
l i = suma de distancias a todos los elementos a la izquierda de x i = (x i -x 1 ) + (x i -x 2 ) + .... + (x i -x i-1 ), y
sl i = suma de cuadrados de distancias a todos los elementos a la izquierda de x i = (x i -x 1 ) ^ 2 + (x i -x 2 ) ^ 2 + .... + (x i -x i -1 ) ^ 2
Tenga en cuenta que, dado i i y sl i , podemos calcular l i + 1 y sl i + 1 en O(1) siguiente manera:
Sea d = x i + 1 -x i . Entonces:
l i + 1 = l i + i d y sl i + 1 = sl i + i d ^ 2 + 2 * i * d
Por lo tanto, podemos calcular todos los l i y sl i en tiempo lineal escaneando de izquierda a derecha. De manera similar, para cada elemento podemos calcular la r i : suma de distancias a todos los elementos a la derecha y sr i : suma de cuadrados de distancias a todos los elementos a la derecha en tiempo lineal. Sumando sr i y sl i para cada i, se obtiene la suma de los cuadrados de distancias horizontales a todos los elementos, en tiempo lineal. Del mismo modo, calcula la suma de cuadrados de distancias verticales a todos los elementos.
Luego podemos escanear a través de la matriz de puntos original y encontrar el punto que minimiza la suma de cuadrados de distancias verticales y horizontales como antes.
Puede resolver el problema como una programación convexa (la función objetivo no siempre es convexa). El programa convexo se puede resolver utilizando un iterativo como L-BFGS. El costo para cada iteración es O (N) y, por lo general, el número de iteraciones requeridas no es grande. Un punto importante para reducir el número de iteraciones requeridas es que sabemos que la respuesta óptima es uno de los puntos en la entrada. Por lo tanto, la optimización se puede detener cuando su respuesta se acerque a uno de los puntos de entrada.
Resolví algo similar para un juez local en línea una vez que usé recocido simulado . Esa fue la solución oficial también y el programa obtuvo AC.
La única diferencia era que el punto que tenía que encontrar no tenía que ser parte de los N puntos dados.
Este era mi código C ++, y N podría ser tan grande como 50000 . El programa se ejecuta en 0.1s en un pentium 4 de 2ghz.
// 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;
}
Entonces creo que es correcto elegir el de tu lista que esté más cercano al (x, y) devuelto por este algoritmo.
Este algoritmo aprovecha lo que dice este párrafo de wikipedia sobre la mediana geométrica:
Sin embargo, es sencillo calcular una aproximación a la mediana geométrica utilizando un procedimiento iterativo en el que cada paso produce una aproximación más precisa. Los procedimientos de este tipo pueden derivarse del hecho de que la suma de las distancias a los puntos de muestra es una función convexa, ya que la distancia a cada punto de la muestra es convexa y la suma de las funciones convexas sigue siendo convexa. Por lo tanto, los procedimientos que disminuyen la suma de las distancias en cada paso no pueden quedar atrapados en un óptimo local.
Un enfoque común de este tipo, llamado algoritmo de Weiszfeld después del trabajo de Endre Weiszfeld, [4] es una forma de mínimos cuadrados re-ponderados iterativamente. Este algoritmo define un conjunto de ponderaciones que son inversamente proporcionales a las distancias desde la estimación actual a las muestras, y crea una nueva estimación que es el promedio ponderado de las muestras de acuerdo con estas ponderaciones. Es decir,
El primer párrafo anterior explica por qué esto funciona: porque la función que intentamos optimizar no tiene mínimos locales, por lo que puede encontrar el mínimo con avidez al mejorarlo iterativamente.
Piense en esto como una especie de búsqueda binaria. Primero, aproximas el resultado. Una buena aproximación será el centro de gravedad, que mi código calcula al leer la entrada. Luego, verás si los puntos adyacentes a esto te dan una mejor solución. En este caso, un punto se considera adyacente si está a una distancia de un step desde su punto actual. Si es mejor, entonces está bien descartar su punto actual, porque, como dije, esto no lo atrapará en un mínimo local debido a la naturaleza de la función que está tratando de minimizar.
Después de esto, tiene la mitad del tamaño del paso, al igual que en la búsqueda binaria, y continúa hasta que tenga lo que considere una aproximación lo suficientemente buena (controlada por la constante eps ).
La complejidad del algoritmo, por lo tanto, depende de la precisión con la que desee que sea el resultado.
Paso 1: ordene la colección de puntos por dimensión x (nlogn)
Paso 2: Calcule la distancia x entre cada punto y todos los puntos HASTA LA IZQUIERDA :
xLDist[0] := 0
for i := 1 to n - 1
xLDist[i] := xLDist[i-1] + ( ( p[i].x - p[i-1].x ) * i)
Paso 3: Calcule la distancia x entre cada punto y todos los puntos A LA DERECHA :
xRDist[n - 1] := 0
for i := n - 2 to 0
xRDist[i] := xRDist[i+1] + ( ( p[i+1].x - p[i].x ) * i)
Paso 4: Suma ambos para obtener la distancia x total de cada punto a los otros puntos N-1
for i := 0 to n - 1
p[i].xDist = xLDist[i] + xRDist[i]
Repita el paso 1,2,3,4 con la dimensión y para obtener
p[i].yDist
El punto con la suma más pequeña de xDist y yDist es la respuesta
Complejidad Total O (nlogn)
Explicación adicional:
La idea es reutilizar la distancia total ya calculada del punto anterior.
Digamos que tenemos 3 puntos ABCD ordenados, vemos que la distancia total izquierda de D a los demás antes de que sea:
AD + BD + CD = (AC + CD) + (BC + CD) + CD = AC + BC + 3CD
En la cual (AC + BC) es la distancia total izquierda de C a las otras anteriores, aprovechamos esto y solo necesitamos calcular ldist(C) + 3CD