DIMENSIONALITY REDUCTION

Think Globally, Fit Locally; LLE

Locally linear embedding (o Local linear embedding ) es otro método de reducción de la dimensión no-lineal. Propuesto por Roweis y Saul (2000), LLE se postula en aquel entonces como un método bastante diferente a lo que se venía proponiendo en el campo de los métodos no lineales de reducción de la dimensionalidad. Sin embargo, las matemáticas usadas eran bastante conocidas desde los años 50. Anteriores métodos se basaban en el escalado multidimensional para preservar las distancias entre observaciones. Estas distancias pueden ser euclídeas o distancias más sofisticadas como la aproximación geodésica en ISOMAP. LLE presenta un enfoque radicalmente diferente (aunque más tarde se demostraría que LLE es de hecho, otra clase de kernel) ya que no intentará preservar las distancias entre observaciones sino la estructura local.

Intuición

LLE se basa en la idea de que un manifold se puede aproximar mediante local linear patches. Asume que una observación $x_i$ junto a sus vecinos se encuentran sobre un subespacio lineal. Por tanto, existirá alguna combinación lineal de los vecinos de $x_i$ que reconstruya $x_i$. Es decir, $x_i$ se podrá aproximar mediante sus vecinos. Aprendiendo las relaciones lineales de los vecinos de las observaciones, LLE puede llegar a “aprender” la estructura global del manifold.

Fuente figura: Karasuyama, M., Mamitsuka, H. Adaptive edge weighting for graph-based learning algorithms. Mach Learn 106, 307–335 (2017)

Dada una matriz $X \in \mathbb{R}^{n \times m}$ el objetivo de LLE será encontrar un embedding de $X$ en una dimensión menor; $Y \in \mathbb{R}^{n \times q}$ donde $q$ < $m$. En general, LLE tiene tres fases, primero se construye un grafo de los k-vecinos más cercanos; se encuentran los pesos que aproximan de forma lineal una observación $x_i$ mediante sus vecinos; finalmente se encuentra una representación en un espacio de una menor dimensión que la original que tome en cuenta los pesos de la anterior etapa. Es decir, la reducción de la dimensión se hace de forma que se preservan las relaciones locales con los vecinos para cada $x_i$.

Fuente figura: S. T. Roweis and L. K. Saul. Think Globally, Fit Locally: Unsupervised Learning of Low Dimensional Manifolds. Journal of Machine Learning Research 4 (2003) 119-155

Algoritmo LLE

K-vecinos más cercanos

La primera fase del algoritmo es encontrar los vecinos más cercanos de cada observación. Una forma muy simple de hacerlo es encontrando los K-vecinos más cercanos usando la distancia euclídea.

\[N(x_i) = \{x_j \in D : d(x_j, x_i) \leq d(x_{L}, x_i )\}\]

Donde $x_L$ es la observación k más cercana de $x_i$. Otra forma de encontrar los vecinos, como ya vimos en DBSCAN, es considerar como vecinos a todas aquellas observaciones que se encuentren a un radio $\epsilon$.

\[N_{Eps}(x_i) = \{x_j \in D : d(x_j, x_i) \leq \epsilon \}\]

Una de las desventajas de esta aproximación es que el grafo resultante puede ser un grafo no-conexo y por tanto tendríamos componentes desconectadas [2]. En este caso, deberíamos considerar a cada componente como un manifold distinto y aplicar LLE por separado en cada componente.

La idea es que el número de vecinos no tiene por qué ser el mismo en cada observación. Pero sigue siendo un gran tema de debate cómo elegir estos vecinos. En nuestro caso, asumiremos que los escogemos según los k-vecinos más cercanos (para evitarnos problemas ya que el grafo será conexo).

Pesos w

Como ya se ha comentado, LLE asume que una observación $x_i$ junto a sus vecinos se encuentran sobre un subespacio lineal. Entonces, existirá alguna combinación de vecinos que reconstruya $x_i$, es decir:

\[x_i = \sum_{j=1}^k w_{ij}x_{N_i(j)}\]

Donde $x_{N_i (j)}$ es el vecino $j$ de la observación $i$. Los pesos $w$ nos servirán para reconstruir $x_i$ en el espacio original y en la reducción de la dimensión. Son los pesos $w$ los que caracterizarán al manifold.

Como queremos reconstruir la observación $x_i$, querremos que el error de reconstrucción sea el menor posible. Podemos formular el problema como:

\[\underset{w}{\text{min}} \; \sum_{i=1}^n \lVert x_i - \sum_{j=1}^k w_{ij}x_{N_i (j)} \rVert ^2\]

Donde $x_{N_i (j)}$ es el vecino $j$ de la observación $i$. La función a optimizar es una suma de las distancias euclídeas al cuadrado entre una observación y su propia reconstrucción.

Invariancia a traslaciones

Como nos centramos únicamente en la forma geométrica del manifold, nos da igual cómo lo movamos en el espacio, por ello, la invariancia a traslaciones es una restricción que queremos imponer al problema. Para que los pesos sean invariantes a traslaciones y reflejen la geometría local, estos deben sumar 1. Ahora, si sumamos un vector $c$ a $x_i$ y a todos sus vecinos, no le va a ocurrir nada a la función que queremos minimizar.

\[x_i + c - \sum_{j=1}^k w_{ij}(x_{N_i (j)} + c) = x_i + c - (\sum_{j=1}^k w_{ij}x_{N_i (j)}) - c =\] \[= x_i - \sum_{j=1}^k w_{ij}x_{N_i (j)}\]

Nótese que si podemos extraer a $c$ del sumatorio es porque $\sum_{j=1}^kw_{ij} = 1$.

Encontrando los pesos

Sea $w_i$ un vector de dimensiones $k \times 1$ que contiene los pesos $w$ para reconstruir la observación $x_i$;

\[w_i = \underbrace{\left[ \begin{array}{ccc} w_{i1}\\ w_{i2}\\ \vdots\\ w_{ik} \end{array} \right] }_{k \times 1}\]

Donde $w_{ik}$ es el peso del vecino $k$ de la observación $i$.

Consideremos ahora el error individual de una observación $x_i$;

\[\varepsilon_i = \lVert x_i - \sum_{j=1}^k w_{ij}x_{N_i (j)} \rVert^2 = \lVert (w_{i1} + \cdots + w_{ik}) x_i - \sum_{j=1}^k w_{ij}x_{N_i (j)} \rVert ^2 =\] \[= \lVert \sum_{j=1}^k w_{ij}x_i - \sum_{j=1}^k w_{ij}x_{N_i (j)} \rVert^2 = \lVert \sum_{j=1}^k w_{ij}(x_i - x_{N_i(j)}) \rVert^2\]

Definimos la matriz $Z$ donde las filas de la matriz serán;

\[z_j = x_i - x_{N_i(j)}\]

Recordemos que estamos considerando el error individual de una observación, por tanto, $Z$ será diferente para cada $x_i$. Podemos reescribir lo anterior como;

\[\varepsilon_i = \lVert \sum_{j=1}^k w_{ij}z_j \rVert^2 = \lVert w_{i1}z_1 + \cdots + w_{ik}z_k \rVert^2 = \lVert w_i^t Z \rVert^2\]

Ahora, teniendo en cuenta que la norma euclídea de cualquier vector fila $m$, es $m m^t$;

\[\varepsilon_i = (w_i^tZ)(w_i^tZ)^t = w_i^t Z Z^t w_i\]

Nótese que $Z Z^t$ es una matriz de dimensiones $k \times k$ de productos escalares de los vecinos, centrada en $x_i$. A esta matriz se la denomina Gram Matrix o matriz de productos escalares (ya la vimos en MDS). La llamaremos $G^{(i)}$ (uso el superíndice $i$ para hacer hincapié en que esta matriz depende de la observación $i$). Los elementos de la matriz serán;

\[G^{(i)}_{jk} = \langle (x_i - x_j), (x_i - x_k) \rangle\]

Entonces, podemos reescribir el error individual de una observación como;

\[\varepsilon_i = w_i^tG^{(i)}w_i\]
Optimización

Queremos minimizar $\varepsilon_i$, teniendo en cuenta además que $\sum_{j=1}^kw_{ij} = 1$. Para resolver este problema de optimización se usan los multiplicadores de Lagrange. Representamos la restricción en forma matricial como;

\[1^tw_i = 1 \rightarrow 1^tw_i - 1 = 0\]

Donde $1$ es una matriz de dimensiones $k \times 1$ formada por unos. Ahora, podemos escribir el lagrangiano como;

\[L(w_i, \lambda) = w_i^tG^{(i)}w_i - \lambda (1^tw_i - 1)\]

Tomando derivadas parciales;

\[\frac{\partial L}{\partial w_i} = 2 G^{(i)}w_i - \lambda1 = 0\] \[\frac{ \partial L}{\partial\lambda} = 1^tw_i - 1 = 0\]

De la primera derivada tenemos que;

\[G^{(i)}w_i = \frac{\lambda}{2}1 \rightarrow w_i = \frac{\lambda}{2}(G^{(i)})^{-1}1\]

Donde $\lambda$ puede tomar un valor aleatorio y luego normalizaríamos $w_i$ para que su suma sea igual a 1.  Y ya tendríamos una expresión para calcular los pesos de $x_i$💃💃 Pero hay veces en las que encontrar la inversa de $G^{(i)}$ es muy costoso. También se puede resolver el siguiente sistema lineal [1] y luego reescalar $w$ para que sume 1;

\[G^{(i)}w_i = 1\]

Donde las k ecuaciones serán;

\[\sum_{k=1}^k G^{(i)}_{jk}w_{ik}\]

Los pesos para cada observación se obtienen de forma individual.

Reducción de la dimensión

Una vez conocemos los pesos $w$, solo nos queda calcular $Y$. Para simplificar, vamos a asumir que $Y$ es una matriz de dimensiones $n \times 1$, o sea, un espacio 1-dimensional. En este espacio 1-dimensional querremos tener una representación de las observaciones que tengan en cuenta los pesos de los vecinos para cada observación ($w_i$). Básicamente, estamos preservando la estructura local de los datos al mantener los pesos que reconstruyen a una observación cualquiera. La función objetivo será la siguiente:

\[\phi(Y) = \sum_{i=1}^n \lVert y_i - \sum_{j=1}^k w_{ij}y_j\rVert^2\]

La función es igual que la anterior con la diferencia que ahora las $w$ sí son conocidas pero las $y$ no lo son.

Restricciones

Al igual que en el caso de la función $\varepsilon$, se deben añadir algunas restricciones a las $y$.

Translational degree of freedom 🏃 Sumar un vector $c$ a $y_i$ no afectará a $\phi$ (la demostración es muy similar a la vista en el caso anterior). Esta invariancia a traslaciones se puede conseguir haciendo que la matriz $Y$ sea una matriz centrada en el origen de coordenadas, esto es;

\[\sum_{i=1}^n y_i = 0 \rightarrow 1^tY = 0\]

Rotational degree of freedom 🌌 La matriz $Y$ también debería poder ser rotada sin que ello afectase a $\phi(Y)$. Para ello debemos añadir la restricción de que las observaciones tengan una matriz de covarianzas igual a la matriz identidad;

\[\frac{1}{n}Y^tY = \underbrace{I}_{q \times q}\]

Donde $q$ en nuestro caso sería igual a 1. Esta condición implica que las componentes principales (sé que no se llaman así) estarán incorreladas, o sea, $Y$ será una matriz ortogonal. ¿No recuerda un poco a PCA? ¿No huele a vectores propios? ¿Será que LLE es algún tipo de kernel? 🤷

Teniendo esto en cuenta, podemos reescribir $\phi(Y)$ como;

\[\phi(Y) = \sum_{i=1}^n \lVert y_i - \sum_{j=1}^k w_{ij}y_j\rVert^2 = \sum_{i=1}^n (y_i - \sum_{j=1}^k w_{ij}y_j)^2 =\] \[= \sum_{i=1}^n (y_i^2 - (2y_i\sum_{j=1}^k w_{ij}y_j) + (\sum_{j=1}^k w_{ij}y_j)^2) =\] \[= \sum_{i=1}^n (y_i^2 - y_i(\sum_{j=1}^k w_{ij}y_j) - (\sum_{j=1}^k w_{ij}y_j)y_i + (\sum_{j=1}^k w_{ij}y_j)^2 )\]

Para escribir la expresión anterior de forma matricial iremos por partes (como hubiese dicho Jack the Ripper);

(1)  $\sum_{i=1}^n y_i^2$

Sea $Y$ una matriz de dimensiones $n \times 1$, podemos reescribir la expresión como $Y^tY$

\[\sum_{i=1}^n y_i^2 = Y^tY = \pmatrix{y_1 & \cdots y_n} \pmatrix{y_1 \cr \vdots \cr y_n} = y_1^2 + \cdots + y_n^2\]
(2)  $\sum_{i=1}^n\sum_{j=1}^k w_{ij}y_j$

Sea $\tilde{W}$ una matriz de dimensiones $n \times n$ con los pesos $w_i$ en las filas;

\[\tilde{W} = \left[ \begin{array}{ccc} w_1\\ w_2\\ \vdots\\ w_n \end{array} \right]\]

Donde $w_i$ serán los pesos de todas las observaciones que reconstruyen $x_i$. Pero, solo habrán valores en las observaciones vecinas, para las demás el valor $w_{ij}$ será cero, esto es;

\[w_i =\left[ \begin{array}{ccc} 0, 0, 0, 0, \cdots, 0, w_{i1}, \cdots, w_{ik},\cdots, 0 \end{array}\right]\]

Donde $w_{ik}$ es el peso del vecino $k$ de la observación $x_i$. Teniendo esto en cuenta, podemos reescribir $\sum_{i=1}^n\sum_{j=1}^k w_{ij}y_j$ como $\tilde{W}Y$. Es decir;

\[\sum_{i=1}^n\sum_{j=1}^k w_{ij}y_j = \tilde{W}Y = \underbrace{ \left[ \begin{array}{ccc} w_1\\ w_2\\ \vdots\\ w_n \end{array} \right] }_{n \times n} \underbrace{ \left[ \begin{array}{ccc} y_1\\ y_2\\ \vdots\\ y_n \end{array} \right] }_{n \times 1}\]

Y tomando (1) y (2), podemos escribir $\phi(Y)$ en forma matricial;

\[\phi(Y) = Y^tY - Y^t(\tilde{W}Y) - (\tilde{W}Y)^tY + (\tilde{W}Y)^t(\tilde{W}Y) =\] \[= Y^t(Y - \tilde{W}Y) - (\tilde{W}Y)^t(Y + \tilde{W}Y) =\] \[= (Y^t - Y^t \tilde{W}^t)(Y - \tilde{W}Y) =\] \[= Y^t(I - \tilde{W}^t)(I - \tilde{W}) Y =\] \[= Y^t(I - \tilde{W})^t(I - \tilde{W}) Y\]

Sea $M = (I - \tilde{W})^t(I- \tilde{W})$. Podemos expresar $\phi(Y)$ como;

\[\phi(Y) = Y^tMY\]
Optimización

Queremos minimizar la función $\phi(Y)$, teniendo en cuenta además las dos restricciones comentadas. Para resolver este problema de optimización se usan los multiplicadores de Lagrange.

\[L(Y, \mu, \psi) = Y^tMY - \mu(n^{-1}Y^tY - 1) - \psi(1^tY)\] \[\frac{\partial \; L(Y, \mu, \psi)}{\partial \; Y} = 2MY - 2\mu n^{-1} Y - \psi 1^t = 0\]

Como la restricción de la invariancia a traslaciones es lineal en $Y$, podemos olvidarnos de $\psi 1^t$ y tenemos que;

\[MY = \frac{\mu}{n}Y\]

De donde se observa que $Y$ es un vector propio de $M$ (quién lo iba a decir… 😱 ). La matriz $M$ es semidefinida positiva de dimensiones $n \times n$. Por tanto, tendrá $n$ valores propios y $n$ vectores propios (de dimensiones $n \times 1 $ cada uno). Como es semidefinida positiva, los valores propios serán reales y no nulos.

Entonces, escogeremos los $q$ vectores propios asociados a los $q$ valores propios más pequeños de la matriz $M$ (ya que estamos minimizando $\phi(Y)$) para construir la matriz $Y$. Sin tener en cuenta el vector propio de unos asociado al valor propio 0 (trivial).

Algunas consideraciones

1. Cuando el número de vecinos $k$ es mayor que $m$ (dimensión de los vectores $x_i$), no hay una única solución al sistema de ecuaciones ya que tenemos más incógnitas que ecuaciones. Por ejemplo, para un $m = 2$ y un $k = 4$, tendríamos que resolver el siguiente sistema lineal;

\[\left\{ \begin{array}{ccccc} w_{i1} x_{11} + w_{i2} x_{21} +w_{i3} x_{31} + w_{i4} x_{41} = & x_{i1} \\ w_{i1} x_{12} + w_{i2} x_{22} +w_{i3} x_{32} + w_{i4} x_{42} = & x_{i2} \\ \end{array} \right.\]

En este caso, el LLE clásico aplica un parámetro de regularización.

2. LLE mantendrá cerca a aquellas observaciones que se encuentren cerca en la dimensión original. Sin embargo, para observaciones que estén lejos, LLE no nos asegura que también lo estén en el espacio q-dimensional.

3. No hay forma de estimar la dimensionalidad de los datos. En PCA podíamos observarla en los valores propios de la matriz de covarianzas. En cambio, en LLE, los valores propios no reflejan la dimensionalidad de los datos.

Referencias