Introducción
El problema del par de puntos más cercanos consiste en dados \(n\) puntos, encontrar el par con menor distancia euclidiana. Aunque parezca trivial, no lo es tanto.
Por ejemplo, dados los siguientes puntos, qué par es más cercano?
\[\begin{align*} (72, 34)&& (24, 45)&& (73, 39) &&(29, 13) \\ (67, 27)&& (85, 27)&& (74, 88) &&(50, 37) \\ (31, 41)&& (5, 44) &&(79, 14) &&(86, 47) \\ (61, 15)&& (5, 90) &&(26, 18) &&(12, 13) \end{align*}\]La respuesta correcta es \((72, 34), (73, 39).\) La única forma (de primeras) de comprobarlo es probar todas las parejas e ir calculando las distancias para quedarse con la menor. Esto son \(\mathcal{O}(n^2)\) operaciones, podemos hacerlo mejor?
Lo cierto es que sí. Existe un algoritmo que usa la técnica divide y vencerás, dividiendo el problema por la mitad recursivamente y que lo acaba resolviendo en un total de \(\mathcal{O}(n \log n)\) operaciones, lo que resulta una mejora significativa.
Podemos mejorarlo todavía más? Lo cierto es que sí, pero hasta cierto punto. Existe una cota inferior de \(\mathcal{O}(n \log \log n)\), por lo que no podemos mejorarlo más allá de eso.
Entonces, hasta aquí esta entrada de blog? Pues lo cierto es que esta cota es para algoritmos deterministas, por lo que no podremos tener un algoritmo que siempre sea más rápido y funcione, pero nada nos impide tener un algoritmo con una cierta probabilidad de fallo, o una cierta (pero arbitariamente pequeña) probabilidad de acabar siendo tan lento como \(\mathcal{O}(n \log \log n)\).
Un algoritmo aleatorio
-
En primer lugar, tomaremos \(\mathcal{O}(n)\) pares de puntos al azar (todas las elecciones independientemente y uniformemente). Calcularemos las distancias de esos pares y nos quedaremos con la mínima, que la llamaremos \(d\).
-
Ahora dividiremos el plano en bloques de tamaño \(d \times d\), y agruparemos todos los puntos en sus respectivos bloques (usando, por eficiencia, tablas hash).
-
Calculamos distancias solo entre pares en un mismo bloque o en bloques adyacentes y cogemos la mínima de todas estas y de \(d\).
-
Fin del algoritmo.
Implementación
Enlace a la implementación en C++.
Análisis
Según esta entrada en Wikipedia, el tiempo esperado de este algoritmo será \(\mathcal{O}(n)\) y siempre dará la respuesta correcta. A continuación vamos a hacer nuestro análisis.
La respuesta siempre será correcta
En el algoritmo habíamos encontrado un par de puntos con distancia \(d\). Al hacer bloques de \(d \times d\), no queremos mirar puntos separados por un bloque entre medio porque entonces ya sabemos que la distancia será \(d\) o más, y buscamos puntos más cercanos. Es por esto que ya solo nos hace falta mirar dentro de un bloque o entre bloques adyacentes.
El tiempo depende en ensencia del tamaño de los bloques individualmente
Como hay \(n\) puntos, como mucho hay \(n\) bloques no vacíos. Digamos que el primero tiene \(b_1\) puntos, el segundo \(b_2\), etc… Entonces el tiempo de ejecución es \(T(n) = \mathcal{O}\bigl(\sum b_i^2\bigr)\). Para un bloque \(i\) solo estamos mirando sus propios pares, de los que hay \(\mathcal{O}\bigl(b_i^2\bigr)\). Pares con bloques adyacentes hay \(\mathcal{O}\bigl(b_i b_j\bigr)\). Hay como mucho 8 bloques adyacentes, y para cada uno de estos, \(\mathcal{O}\bigl(b_i b_j\bigr) \le\) \(\mathcal{O}\bigl(\max(b_i, b_j)^2\bigr) \le\) \(\mathcal{O}\bigl(b_i^2 + b_j^2\bigr),\) por lo que la suma de los bloques adyacentes es menor o igual a la de los bloques individualmente y solo tenemos que fijarnos cada bloque individualmente.
Bloques grandes
Supongamos que hemos ejecutado la fase aleatoria y hemos calculado una cierta \(d\) y partido el plano en bloques \(d \times d\).
Para cualquier \(\varepsilon > 0\) previamente fijado (tan pequeño como se quiera), llamamos bloque grande a uno que tenga \(\sqrt{n}^{1+\varepsilon}\) o más puntos. Los bloques grandes son problemáticos para la complejidad lineal, ya que si el tamaño del bloque es \(\sqrt{n}^{1+\varepsilon}\) nuestro algoritmo sería como mínimo \(\mathcal{O}(n^{1+\varepsilon})\), y podría ser \(\mathcal{O}(n^2)\) en el peor caso (si tenemos bloques de tamaño \(n/2\), por ejemplo).
Digamos que incialmente en el algoritmo fijamos una constante \(C\) tal que cogemos \(Cn\) puntos iniciales con los que calcular \(d\). Calculemos la probabilidad de tener un bloque grande, definiendo para ello distintos eventos:
- \(A\): Se elige un punto del bloque grande
- \(B\): Se elige un par del bloque grande
- \(W\): Se elige algun par del bloque grande en las \(Cn\) rondas
- \(\overline{B}\): No se elige un par del bloque grande (opuesto a \(B\))
- \(\overline{W}\): No se elige ningun par del bloque grande en las \(Cn\) rondas (opuesto a \(W\))
Vemos que cuando \(n \to \infty\), tenemos que \(\Pr(W) \to 1\). Esto nos dice que cuando \(n\) crece tenemos esencialmente la certeza de que saldrán pares de cualquier bloque grande. De hecho, aunque no lo analizamos más allá, que la probabilidad sea \(1\) nos dice que seguramente saldrán muchas veces pares de un mismo bloque. Hay muy pocos pares posibles dentro de un bloque que no tengan una distancia \(< d\), por lo que lo esperable es que acabemos encontrando un par de puntos con distancia \(< d\), pero esto contradice el hecho de que \(d\) es la distancia más pequeña encontrada. Es decir, al presuponer que habrá bloques grandes, tenemos un escenario en el que hay una contradicción con probabilidad que tiende a \(1.\) Esto nos dice que con probabilidad \(1\) no habrá ningún bloque grande cuando \(n\) tiende a infinito.
Tiempo prácticamente lineal
Para cualquier \(\varepsilon > 0\) que queramos, podemos decir que el algoritmo será \(\mathcal{O}\bigl(n^{1+\varepsilon}\bigr)\). Tenemos una cota prácticamente lineal.
Un refinamiento probabilístico de la cota
Si en el desarrollo anterior (donde los bloques grandes) tomamos \(\varepsilon = 0,\) considerando bloques de tamaño \(\sqrt{n}\) tenemos que \(\Pr(W) = 1 - \bigl(1 - \frac{1}{n} \bigr)^{Cn}\), y cuando \(n \to \infty,\) tenemos que \(\Pr(W) \to 1 - e^{-C}.\)
Esto nos dice que, en general, podemos coger una \(C\) que hace relativamente probable que se coja un par de un bloque \(\sqrt{n}\), lo que genera una contradicción y por tanto, con una \(C\) suficientemente grande, es muy probable (aunque no con probabilidad \(1\)) que no tengamos bloques \(\sqrt{n}\).
Como con probabilidad \(1 - e^{-C}\) no tenemos bloques \(\sqrt{n}\) y el algoritmo es \(\mathcal{O}(n)\), si aceptamos una probabilidad \(\large \epsilon\) de que el algoritmo no sea puramente lineal, podemos decidir qué \(C\) tomar. Por ejemplo, podemos coger \(\epsilon = 1/1000\) y \(C = -\log(\epsilon) \approx 7\) si aceptamos que no sea lineal una de cada $1000$ veces. Aún así ya hemos visto que con probabilidad \(1\) será prácticamente lineal, por lo que en la práctica exhibe un comportamiento indistinguible del lineal.
Es el tiempo esperado de este algoritmo lineal?
Viendo que es prácticamente lineal y que podemos hacer arbitrariamente alta la probabilidad de que sea verdaderamente lineal, no puede ser que estas dos cosas impliquen que el tiempo esperado es lineal?
Lo cierto es que según Wikipedia está demostrado que es lineal, pero no he encontrado esa demostración. Como veremos a continuación, que en este caso sea lineal no implica que estas dos propiedades que hemos visto previamente sean suficientes para decir esto.
Las observaciones previas no son suficientes para decir que la esperanza (en el sentido preciso de esperanza matemática) sea lineal. Hablando sobre esto con Felix Moreno, se le ocurrió un ejemplo de algoritmo que cumple lo mismo que el nuesto, pero la esperanza no es lineal:
Imaginemos un algoritmo que tiene un tiempo de ejecución:
\[T(n) = \begin{cases} n, & \text{ con probabilidad } 1-\frac{1}{\log\log n} \\ n \log n, & \text{ con probabilidad } \frac{1}{\log\log n} \\ \end{cases}\]Entonces,
\[\mathbb{E}[T(n)] = n \Bigl( 1-\frac{1}{\log\log n} \Bigr) + n \log n \Bigl(\frac{1}{\log\log n} \Bigr)\]Este algoritmo sería \(\mathcal{O}(n^{1+\varepsilon})\) para cualquier \(\varepsilon>0\) y \(\mathcal{O}(n)\) con una probabilidad que tiende a \(1\). Sin embargo, la esperanza sigue siendo \(\mathcal{O}\bigl(n \frac{\log n}{\log \log n}\bigr) >\) \(\mathcal{O}(n).\)
Un algoritmo más sencillo de analizar
Veremos ahora un algoritmo no tan sencillo y que en la práctica es menos eficiente, pero que nos permite un análisis más sencillo. [Kleinberg, Tardos. Algorithm Design, 12.7]
- Permutamos los puntos aleatoriamente.
- Definimos $\delta := \text{dist}(p_1, p_2).$
- Particionamos el plano en cuadrados de lado $\delta / 2$.
- Para $i = 1,2,\dots,n:$
- Encontramos el cuadrado que corresponde a $p_i.$
- Miramos los 25 cuadrados cercanos al cuadrado de $p_i,$ los que están entre $\pm2$ posiciones alrededor como máximo (sin almacenar puntos inicialmente en los cuadrados).
- Si hay algun $p_j$ almacenado en esos cuadrados tal que $\text{dist}(p_i, p_j) < \delta,$ entonces:
- Redefinimos $\delta := \text{dist}(p_i, p_j).$
- Recalculamos la partición del plano con la nueva $\delta.$
- Almacenamos los puntos $p_1,\dots,p_i$ en los cuadrados que les corresponden.
- En caso contrario:
- Almacenamos $p_i$ en el cuadrado que le corresponde.
Análisis
Aunque inicialmente parece que este segundo algoritmo debería ser muy lento porque podría estar recalculando los cuadrados con puntos constantemente, vamos a ver que la probabilidad de que eso pase decrece lo suficientemente rápido para que el tiempo esperado sea lineal.
Llamemos $X_i$ a la variable aleatoria binaria que es $1$ cuando el punto $p_i$ provoca un cambio en $\delta$ y por tanto que se recalculen los cuadrados. El número total de veces que se calcula para cada punto su cuadrado es
\[n + \sum_{i=0}^{n-1} i \, X_i\]Pero resulta que \(\Pr(X_i = 1) \le 2/i.\) Esto se debe a que los puntos están ordenados aleatoriamente y la probabilidad de que la distancia mínima caiga en un par con el punto $i$ pero en ningún par anterior es como mucho $2/i$. Usando esta cota tenemos que
\[\mathbb{E}\Bigl[ n + \sum_{i=0}^{n-1} i \, X_i \Bigr] = n + \sum_{i=0}^{n-1} i \Pr(X_i = 1) \le n + \sum_{i=0}^{n-1} i \, \frac{2}{i} = 3n\]Por lo que el tiempo esperado es \(\mathcal{O}(n).\)
Conclusión
Tenemos un algoritmo que es esencialmente lineal, mientras que los métodos deterministas tenían una cota inferior de \(\mathcal{O}(n\log \log n)\). No solo hemos ganado eficiencia, sino también simplicidad tanto para entender el algoritmo como para implementarlo.
Luego hemos visto un algoritmo que aunque parecería más lento (y en la práctica lo es, aunque no asintóticamente), resulta que simplifica mucho su análisis, por lo que a veces podemos encontrar cierta distancia entre la teoría y la implementación.
Por último que, en ocasiones, probar cosas aleatoriamente puede darnos una ventaja imposible de conseguir con métodos deterministas.