Problema del vector más cercano (CVP) en un retículo

Definición del problema: 

Sea $\varLambda$ un retículo completo sobre $\mathbb{R}^n$ y sea $\vec{w} \in \mathbb{R}^n$.

Queremos encontrar un vector $\vec{z} \in \varLambda$ tal que:

$$\|\vec{w}-\vec{z}\| = mín\{\|\vec{w}-\vec{b}\|, \forall \vec{b} \in \varLambda\}$$

Es decir, dado un vector del espacio queremos encontrar el vector del retículo más cercano al primero. Esto se conoce como el Problema del Vector más Cercano (CVP, por las siglas en inglés "Closest Vector Problem").

 

Teniendo en cuenta que la región fundamental de un retículo es un n-paralelepípedo, es inmediato ver la relación con las entradas anteriores (1 y 2) en las que se plantea encontrar el vértice del n-paralelepípedo más cercano a un punto.

Resumiendo, para dimensiones altas el método exhaustivo consistente en calcular la distancia con todos los vértices del n-paralelepípedo tiene un coste computacional demasiado elevado. 

Por otra parte, el algoritmo de Babai solo ofrece garantías de obtener la solución correcta si la base es lo "suficientemente ortogonal". Pero a esto además le añadimos que con bases que tienen mucho defecto de ortogonalidad (muy "lejos" de ser ortogonales) resulta que el vector más cercano ni tan siquiera tiene porqué ser vértice del n-paralelepípedo que lo contiene.

Puede verse que hay un punto naranja (del retículo) más próximo a $Q$ que los vértices del n-paralelepípedo que lo contiene

Con todo lo anterior, podemos concluir lo siguiente sobre la resolución del problema del vector más próximo en un retículo: una persona que conozca una base "bastante próxima" a ser ortogonal podrá realizar los cálculos de manera relativamente sencilla, mientras que otra persona que conozca una base "suficientemente lejos" de ser ortogonal tendrá mucha más dificultad para resolverlo.

Algoritmo de Babai con Python

 En la entrada anterior expliqué y ejemplifiqué el Algoritmo de Babai para encontrar el vértice del n-paralelepípedo más cercano a un punto.

Podemos realizar un script en Python para dicho algoritmo. Utilizando la librería NumPy queda muy sencillo:

 

import numpy as np

def babai(B, w):
    # Calculamos los coeficientes del vector w en la base B
    v = np.linalg.solve(B, w)
    # Devolvemos el producto matricial de B por las coordenadas de v redondeadas
    return B @ np.round(v).astype(int) 

 

By www.python.org - www.python.org, GPL,

 https://commons.wikimedia.org/w/index.php?curid=34991651

 

PD: Podéis encontrar en uno de mis repositorios de Github un script para probar el Algoritmo de Babai y la razón de Hadamard. 

Vértice de un n-paralelepípedo más cercano a un punto. Algoritmo de Babai.

 En una entrada anterior planteé el siguiente problema:

 Sea $\mathcal{F}$ un n-paralelepípedo (sólido) generado por el conjunto de vectores linealmente independientes $B = \{\vec{b}_1,...,\vec{b}_n\}$. Sea $P\in\mathcal{F}$ un punto de dicho n-paralelepípedo. ¿Cuál es el vértice del n-paralelepípedo más cercano a $P$?

Expliqué y ejemplifiqué una posible estrategia consistente en el método exhaustivo, es decir, comprobar todos los posibles vértices. Concluí la entrada con la observación de que el número de vértices crece de manera exponencial con la dimensión y esto supone un problema de coste computacional en dimensiones grandes.

En esta entrada vamos a trabajar otra posible estrategia, basada en el algoritmo de Babai, desarrollado por el matemático húngaro László Babai.

El fundamento detrás del procedimiento es relativamente sencillo.

Si $P\in\mathcal{F}$ entonces $P$ se puede escribir como combinación lineal de la base $B$:

$$P = \sum_{i=1}^{n}{a_i·\vec{b}_i}$$

Lo que hacemos es redondear al entero los coeficientes $0 \leq a_i \leq 1$ de la combinación lineal:

$$z_i = redondear(a_i)$$

Entonces, el vértice del n-paralelepípedo más próximo a $P$ viene dado por:

$$V = \sum_{i=1}^{n}{z_i·\vec{b}_i}$$

 

Ejemplo: Hagamos este método con el mismo ejemplo que utilizamos en el método exhaustivo, en $\mathbb{R}^2$.

$B=\{(30, 3), (5, 0)\}$ y el punto del paralelogramo $Q=(25.5, 2.5)$. Queremos hallar el vértice del paralelogramo más cercano a $Q$.

 Primero calcularemos la combinación lineal de la base $B$ que da el punto $Q$:

$$\vec{a} = B^{-1} · \overrightarrow{OQ} =\begin{pmatrix}30 & 5\\3 & 0\end{pmatrix}^{-1} · \begin{pmatrix}25.5 \\ 2.5\end{pmatrix} \simeq \begin{pmatrix}0.8333333333333333 \\ 0.10000000000000053\end{pmatrix}$$

Ahora redondeamos los coeficientes obtenidos:

$$\vec{z} = \begin{pmatrix}1 \\ 0\end{pmatrix}$$

Y calculamos a qué vértice corresponde dicha combinación lineal:

$$\overrightarrow{OV} = B · \vec{z} = \begin{pmatrix}30 & 5\\3 & 0\end{pmatrix} · \begin{pmatrix}1 \\ 0\end{pmatrix} = \begin{pmatrix}30 \\ 3\end{pmatrix}$$

Es decir, el vértice más cercano al punto $Q$ es (30, 3). (Coincide con el obtenido por el método exhaustivo)

Los cálculos del ejemplo han sido realizados con wxMaxima:

 

Desde un punto de vista geométrico e informal podríamos decir que el procedimiento anterior es equivalente a dividir con "hiperplanos medianeros" el n-paralelepípedo y seleccionar el vértice del "trozo" en el que queda el punto.

 

 El algoritmo de Babai rebaja considerablemente el coste computacional, dado que no necesita recorrer todos los vértices del n-paralelepípedo y se basa en operaciones bastante rápidas desde el punto de vista del cálculo computacional.

Sin embargo, el algoritmo de Babai presenta un gran problema: solo en determinadas condiciones la solución que devuelve es la correcta.

Veamos esto con un ejemplo. Seguimos utilizando la misma base que en el ejemplo anterior, pero cambiamos el punto a "aproximar" (encontrar el vértice más cercano) por $R=(30, 2.7)$.

Los cálculos del algoritmo devuelven el vértice $(35, 3)$ como el más cercano al punto $R$.

 

Sin embargo, esto no es cierto. El vértice más cercano es claramente el $(30, 3)$. 

 

Puede observarse en la representación gráfica que, tras trazar las medianas, el punto cae en el "trozo" del vértice $(35, 3)$ y por eso es la respuesta que devuelve el algoritmo de Babai. Pero no es el vértice más cercano.

Llegados a este punto lo razonable es descartar por completo este algoritmo porque ni tan siquiera devuelve la solución correcta.

Sin embargo, el algoritmo sí funciona con fiabilidad cuando la base es ortogonal (o "cerca" de ser ortogonal).

¿Y esto tiene alguna utilidad práctica? Habrá que esperar a las siguientes entradas para averiguarlo. 

Vértice de un n-paralelepípedo más cercano a un punto. Método exhaustivo.

 

Sea $\mathcal{F}$ un n-paralelepípedo (sólido) generado por el conjunto de vectores linealmente independientes $B = \{\vec{b}_1,...,\vec{b}_n\}$. Sea $P\in\mathcal{F}$ un punto de dicho n-paralelepípedo. ¿Cuál es el vértice del n-paralelepípedo más cercano a $P$?

 Estrategia 1. Método exhaustivo.

Podemos calcular las coordenadas de todos los vértices del n-paralelepípedo.

Por ejemplo, en $\mathbb{R}^3$, los 8 vértices del paralelepípedo son:

$V_1 = \vec{0} = (0,0,0)$

$V_2 = \vec{b}_1$

$V_3 = \vec{b}_2$

$V_4 = \vec{b}_3$

$V_5 = \vec{b}_1 + \vec{b}_2$

$V_6 = \vec{b}_1 + \vec{b}_3$

$V_7 = \vec{b}_2 + \vec{b}_3$

$V_8 = \vec{b}_1 + \vec{b}_2 + \vec{b}_3$ 

Después calculamos la distancia de cada uno de ellos al punto $P$

$$d_i = \|P-V_i\|$$ 

 Y nos quedamos con el que tenga menor distancia (puede haber más de uno).

 

Ejemplo: Hagamos este método con un caso muy sencillo en $\mathbb{R}^2$.

$B=\{(30, 3), (5, 0)\}$ y el punto del paralelogramo $Q=(25.5, 2.5)$. Queremos hallar el vértice del paralelogramo más cercano a $Q$.

Los vértices del paralelogramo son:

$V_1=(0, 0); V_2=(30, 3); V_3=(5, 0); V_4=(35, 3)$

Y las correspondientes distancias al punto $Q$:

$d_1 \simeq 25.622255950637914$

$d_2 \simeq 4.527692569068709$

$d_3 \simeq 20.65187642806338$

$d_4 \simeq 9.513148795220223$ 

Por tanto, el vértice con menor distancia al punto $Q$ es $V_2=(30, 3)$.

 

Los cálculos han sido realizados con wxMaxima:

 

Contras del método exhaustivo.

Para dimensiones pequeñas el método exhaustivo es rápido y muy eficiente, pero para dimensiones grandes el coste computacional se incrementa exponencialmente, dado que el número de vértices de un n-paralelogramo es $2^n$.

Por ejemplo, en $\mathbb{R}^{100}$ tenemos $2^{100} = 1267650600228229401496703205376$ vértices.

Medir la ortogonalidad de una base

En la entrada anterior hice un repaso rápido de qué son las bases ortogonales. Ahora voy a tratar el siguiente tema:

¿Hay alguna manera de cuantificar si una base está cerca o lejos de ser ortogonal?

La respuesta es sí.

Una manera de hacerlo es utilizando una aplicación de la desigualdad de Hadamard.

Dada una base $B = \{ \vec{b}_1, ..., \vec{b}_n \}$ del espacio euclídeo $\mathbb{R}^n$:

$$|det(B)| \leq \prod_{i=1}^{n}\|\vec{b}_i\|$$

Desde el punto de vista geométrico, la desigualdad anterior nos aporta una cota superior del volumen del n-paralelepídedo formado por $B$.

Pero para el caso que nos ocupa, es muy útil conocer que dicha desigualdad se convierte en igualdad cuando la base es ortogonal por lo que nos permite definir la siguiente "razón de ortogonalidad de Hadamart":

$$\frac{|det(B)|}{\prod_{i=1}^{n}\|\vec{b}_i\|} $$

 Dicho valor estará comprendido en el intervalo $(0,1]$ y será exactamente 1 cuando la base sea ortogonal. Cuanto más "cerca" esté una base de ser ortogonal, más cerca de 1 estará el valor de la razón anterior. Cuanto más "lejos" esté una base de ser ortogonal, más cerca de 0 estará dicho valor.

Desde un punto de vista intuitivo e informal, estamos diciendo que si fijamos el módulo de los vectores que generan un n-paralelepídedo entonces el volumen es máximo cuando estos son ortogonales (en el momento en el que variamos alguno de los ángulos el volumen del n-paralelepípedo baja). Por tanto, si comparamos el volumen del n-paralelepípedo generado por los vectores de la base con el volumen máximo (caso ortogonal) podemos tener una idea de si la base está "cerca" o "lejos" de ser ortogonal.

Vamos a calcular dicha razón en los ejemplos de la entrada anterior en $\mathbb{R}^3$.

Comenzamos con la base canónica $C = \{(1,0,0),(0,1,0),(0,0,1)\}$:

 

 $|det(C)|=1$

 $\|(1,0,0)\|=\sqrt{1^2+0^2+0^2}=1$ 

 $\|(0,1,0)\|=\sqrt{0^2+1^2+0^2}=1$ 

 $\|(0,0,1)\|=\sqrt{0^2+0^2+1^2}=1$

Por tanto,

$$\frac{|det(C)|}{\prod_{i=1}^{n}||\vec{c}_i||} = \frac{1}{1} = 1$$

Es decir, como el valor da 1 podemos afirmar que la base es ortogonal.

 

En el caso de la otra base utilizada, $D = \{(1,0,0),(1,1,1),(0,2,1)\}$ 

 $|det(D)|=1$

 $\|(1,0,0)\|=\sqrt{1^2+0^2+0^2}=1$ 

 $\|(1,1,1)\|=\sqrt{1^2+1^2+1^2}=\sqrt{3}$ 

 $\|(0,2,1)\|=\sqrt{0^2+2^2+1^2}=\sqrt{5}$

 Por tanto,

$$\frac{|det(D)|}{\prod_{i=1}^{n}\|\vec{d}_i\|} = \frac{1}{\sqrt{15}} = \frac{\sqrt{15}}{15} \simeq 0.2581$$

Da un valor alejado de 1 por lo que podemos afirmar que la base "está lejos" de ser ortogonal. 

 

En realidad resultaría de utilidad tener una "medida" normalizada y comparable entre dimensiones. Si añadimos una raíz n-ésima tendremos una especie de "media geométrica del grado de normalidad por vector". Así pues, es más frecuente utilizar la siguiente expresión cuando se habla de la razón de ortogonalidad de Hadamart:

$$H(B) = \sqrt[n]{\frac{|det(B)|}{\prod_{i=1}^{n}\|\vec{b}_i\|}} $$

Otros autores hablan del defecto de ortogonalidad de una base definiéndola de la siguiente manera:

$$\vartriangle(B) = \frac{\prod_{i=1}^{n}\|\vec{b}_i\|}{|det(B)|} $$

de manera que cuánto mayor es el cociente "menor" es la ortogonalidad de la base (que sigue siendo 1 para bases ortogonales).

PD: Entrada de Wikipedia sobre el matemático francés Jacques Hadamard. 

Región fundamental de un retículo

 Sea $\varLambda$ un retículo completo sobre $\mathbb{R}^n$. Sea $B$ una base que genera el retículo $\varLambda$.

Llamamos región fundamental del retículo $\varLambda(B)$ al conjunto:

$$\mathcal{F}(\varLambda(B)) = \{ \sum_{i=1}^{n}t_i·\vec{b_i} \mid t_i \in [0,1) \wedge \vec{b_i} \in B; \forall i \in \{1,...,n\} \}$$

Dicha región fundamental también es conocida como dominio fundamental o paralelepípedo fundamental.

A modo de ejemplo, retomamos el retículo completo sobre $\mathbb{R}^2$ generado por $B = \{ (5, 3), (5, 0) \}$

 

Su región fundamental tiene la siguiente representación gráfica (importante observar la frontera):

 

Una vez definida la región fundamental, podemos realizar las siguientes observaciones:

• El $\vec{0}$ siempre pertenece a la región fundamental, independientemente de la base que genera el retículo $\varLambda$.
• Excluyendo el $\vec{0}$, los vectores de la región fundamental NO pertenecen al retículo $\varLambda$.
• Si cambiamos la base generadora del retículo, la región fundamental asociada también cambia. 

Si en el mismo retículo completo del ejemplo anterior cogemos la base $C = \{(0,3), (5,0)\}$, la región fundamental queda representada de la siguiente manera:

Sin embargo,

• aunque cambie la región fundamental, el "volumen n-dimensional" de la región fundamental (n-paralelepípedo o paralelótopo) es invariante. ¿Por qué? Porque se calcula con el determinante y se explicó un resultado al respecto en la anterior entrada.
• Con la región fundamental podemos teselar todo el espacio $\mathbb{R}^n$.

Teselado de $\mathbb{R}^2$ con la región fundamental $\mathcal{F}(\varLambda(C))$.