¿Por qué Andrew Ng prefiere utilizar la SVD y no la EIG de la matriz de covarianza para hacer el ACP?

Question

Estoy estudiando PCA a partir del curso Coursera de Andrew Ng y otros materiales. En el curso de PNL de Stanford cs224n's primer encargo y en el Vídeo de la conferencia de Andrew Ng En el caso de las matrices de covarianza, se realiza la descomposición del valor singular en lugar de la descomposición de los vectores propios, y Ng incluso dice que la SVD es numéricamente más estable que la descomposición de los vectores propios.

Según tengo entendido, para el ACP debemos hacer la SVD de la matriz de datos de (m,n) tamaño, no de la matriz de covarianza de (n,n) tamaño. Y la descomposición de los vectores propios de la matriz de covarianza.

¿Por qué hacen la SVD de la matriz de covarianza y no de la matriz de datos?

Answer 1

Amoeba ya dio una buena respuesta en los comentarios, pero si quieres un argumento formal, aquí va.

La descomposición del valor singular de una matriz $A$ es $A=U\Sigma V^T$ donde las columnas de $V$ son vectores propios de $A^TA$ y las entradas diagonales de $\Sigma$ son los raíces cuadradas de sus valores propios, es decir $\sigma_{ii}=\sqrt{\lambda_i(A^TA)}$ .

Como sabes, los componentes principales son las proyecciones ortogonales de tus variables sobre el espacio de los vectores propios de la matriz de covarianza empírica $\frac{1}{n-1}A^TA$ . La varianza de los componentes viene dada por sus valores propios, $\lambda_i(\frac{1}{n-1}A^TA)$ .

Considere una matriz cuadrada cualquiera $B$ , $\alpha \in \mathbb R$ y un vector $v$ tal que $Bv=\lambda v$ . Entonces

1. $B^kv=\lambda^kv$
2. $\lambda(\alpha B) = \alpha\lambda( B)$

Definamos $S=\frac{1}{n-1}A^TA$ . El SVD de $S$ calculará la eigendecomposición de $S^TS=\frac{1}{(n-1)^2}A^TAA^TA$ para rendir

1. los vectores propios de $(A^TA)^TA^TA=A^TAA^TA$ que por la propiedad 1 son los de $A^TA$
2. el raíces cuadradas de los valores propios de $\frac{1}{(n-1)^2}A^TAA^TA$ que por la propiedad 2, luego 1, luego 2 de nuevo, son $\sqrt{\frac{1}{(n-1)^2} \lambda_i(A^TAA^TA)} = \sqrt{\frac{1}{(n-1)^2} \lambda_i^2(A^TA)} = \frac{1}{n-1}\lambda_i(A^TA) = \lambda_i(\frac{1}{n-1}A^TA)$ .

¡Voilà!

En cuanto a la estabilidad numérica, habría que averiguar cuáles son los alogrados empleados. Si te animas, creo que son las rutinas LAPACK que utiliza numpy:

• numpy.linalg.eig
• numpy.linalg.svd

Actualización: En cuanto a la estabilidad, la implementación de la SVD parece utilizar un enfoque de dividir y conquistar, mientras que la eigendecomposición utiliza un algoritmo QR simple. No puedo acceder a algunos documentos relevantes de SIAM de mi institución (culpa de los recortes de investigación), pero he encontrado algo que podría apoyar la evaluación de que la rutina SVD es más estable.

En

Nakatsukasa, Yuji, y Nicholas J. Higham. "Algoritmos espectrales estables y eficientes de divide y vencerás para la descomposición simétrica de valores propios y la SVD". SIAM Journal on Scientific Computing 35.3 (2013): A1325-A1349.

comparan la estabilidad de varios algoritmos de valores propios, y parece que el enfoque de dividir y conquistar (¡utilizan el mismo que numpy en uno de los experimentos!) es más estable que el algoritmo QR. Esto, junto con las afirmaciones en otros lugares de que los métodos de D&C son realmente más estables, apoya la elección de Ng.

Answer 2

@amoeba tuvo excelentes respuestas a las preguntas del PCA, incluyendo este sobre la relación de SVD con PCA. Respondiendo a su pregunta exacta voy a hacer tres puntos:

• matemáticamente no hay diferencia si se calcula el ACP sobre la matriz de datos directamente o sobre su matriz de covarianza
• la diferencia se debe puramente a la precisión y complejidad numérica. Aplicar la SVD directamente a la matriz de datos es numéricamente más estable que a la matriz de covarianza
• La SVD puede aplicarse a la matriz de covarianza para realizar el PCA u obtener los valores propios, de hecho, es mi método favorito para resolver los problemas propios

Resulta que la SVD es más estable que los procedimientos típicos de descomposición de valores propios, especialmente, para el aprendizaje automático. En el aprendizaje automático es fácil acabar con regresores muy colineales. La SVD funciona mejor en estos casos.

Aquí está el código de Python para demostrar el punto. He creado una matriz de datos altamente colineal, he obtenido su matriz de covarianza y he intentado obtener los valores propios de esta última. La SVD sigue funcionando, mientras que la descomposición eigen ordinaria falla en este caso.

import numpy as np
import math
from numpy import linalg as LA

np.random.seed(1)

# create the highly collinear series
T = 1000
X = np.random.rand(T,2)
eps = 1e-11
X[:,1] = X[:,0] + eps*X[:,1]

C = np.cov(np.transpose(X))
print('Cov: ',C)

U, s, V = LA.svd(C)
print('SVDs: ',s)

w, v = LA.eig(C)
print('eigen vals: ',w)

La salida:

Cov:  [[ 0.08311516  0.08311516]
 [ 0.08311516  0.08311516]]
SVDs:  [  1.66230312e-01   5.66687522e-18]
eigen vals:  [ 0.          0.16623031]

Actualización

Respondiendo al comentario de Federico Poloni, aquí está el código con la prueba de estabilidad de SVD vs Eig en 1000 muestras aleatorias de la misma matriz anterior. En muchos casos Eig muestra 0 valor propio pequeño, lo que llevaría a la singularidad de la matriz, y SVD no lo hace aquí. SVD es aproximadamente dos veces más preciso en la determinación del valor propio pequeño, que puede o no ser importante dependiendo de su problema.

import numpy as np
import math
from scipy.linalg import toeplitz
from numpy import linalg as LA

np.random.seed(1)

# create the highly collinear series
T = 100
p = 2
eps = 1e-8

m = 1000 # simulations
err = np.ones((m,2)) # accuracy of small eig value
for j in range(m):
    u = np.random.rand(T,p)
    X = np.ones(u.shape)
    X[:,0] = u[:,0]
    for i in range(1,p):
        X[:,i] = eps*u[:,i]+u[:,0]

    C = np.cov(np.transpose(X))

    U, s, V = LA.svd(C)

    w, v = LA.eig(C)

    # true eigen values
    te = eps**2/2 * np.var(u[:,1])*(1-np.corrcoef(u,rowvar=False)[0,1]**2)
    err[j,0] = s[p-1] - te
    err[j,1] = np.amin(w) - te

print('Cov: ',C)
print('SVDs: ',s)
print('eigen vals: ',w)
print('true small eigenvals: ',te)

acc = np.mean(np.abs(err),axis=0)    
print("small eigenval, accuracy SVD, Eig: ",acc[0]/te,acc[1]/te)

La salida:

Cov:  [[ 0.09189421  0.09189421]
 [ 0.09189421  0.09189421]]
SVDs:  [ 0.18378843  0.        ]
eigen vals:  [  1.38777878e-17   1.83788428e-01]
true small eigenvals:  4.02633695086e-18
small eigenval, accuracy SVD, Eig:  2.43114702041 3.31970128319

Aquí el código funciona. En lugar de generar la matriz de covarianza aleatoria para probar las rutinas, estoy generando la matriz de datos aleatorios con dos variables: $$x_1=u\\ x_2=u+\varepsilon v$$ donde $u,v$ - variables aleatorias uniformes independientes. Así, la matriz de covarianza es $$\begin{pmatrix} \sigma_1^2 & \sigma_1^2 + \varepsilon \rho \sigma_1 \sigma_2\\ \sigma_1^2 + \varepsilon \rho \sigma_1 \sigma_2 & \sigma_1^2 + 2 \varepsilon \rho \sigma_1 \sigma_2 + \varepsilon^2 \sigma_2^2\sigma^2\end{pmatrix}$$ donde $\sigma_1^2,\sigma_2^2,\rho$ - varianzas de los uniformes y el coeficiente de correlación entre ellos.

Su valor propio más pequeño: $$\lambda= \frac 1 2 \left(\sigma_2^2 \varepsilon^2 - \sqrt{\sigma_2^4 \varepsilon^4 + 4 \sigma_2^3 \rho \sigma_1 \varepsilon^3 + 8 \sigma_2^2 \rho^2 \sigma_1^2 \varepsilon^2 + 8 \sigma_2 \rho \sigma_1^3 \varepsilon + 4 \sigma_1^4} + 2 \sigma_2 \rho \sigma_1 \varepsilon + 2 \sigma_1^2\right)$$ El valor propio pequeño no se puede calcular simplemente enchufando el $\varepsilon$ en la fórmula debido a su limitada precisión, por lo que es necesario que Taylor la amplíe: $$\lambda\approx \sigma_2^2 \varepsilon^2 (1-\rho^2)/2$$

Corro $j=1,\dots,m$ simulaciones de las realizaciones de la matriz de datos, calcular los valores propios de la matriz de covarianza simulada $\hat\lambda_j$ y obtener los errores $e_j=\lambda-\hat\lambda_j$ .

Answer 3

Para los usuarios de Python, me gustaría señalar que para las matrices simétricas (como la matriz de covarianza), es mejor utilizar numpy.linalg.eigh en lugar de una función general numpy.linalg.eig función.

eigh es de 9 a 10 veces más rápido que eig en mi ordenador (independientemente del tamaño de la matriz) y tiene mejor precisión (según la prueba de precisión de @Aksakal).

No me convence la demostración del beneficio de la precisión de la SVD con pequeños valores propios. La prueba de @Aksakal es de 1 a 2 órdenes de magnitud más sensible al estado aleatorio que al algoritmo (intente trazar todos los errores en lugar de reducirlos a un máximo absoluto). Significa que pequeños errores en la matriz de covarianza tendrán un mayor efecto en la precisión que la elección de un algoritmo de eigendecomposición. Además, esto no está relacionado con la pregunta principal, que es sobre PCA. Los componentes más pequeños se ignoran en PCA.

Se puede hacer un argumento similar sobre la estabilidad numérica. Si tengo que utilizar el método de la matriz de covarianza para el ACP, lo descompondría con eigh en lugar de svd . Si falla (lo que aún no se ha demostrado aquí), probablemente merezca la pena replantearse el problema que se intenta resolver antes de empezar a buscar un algoritmo mejor.

Answer 4

Ya se han dado grandes respuestas a tus preguntas, así que no añadiré muchas cosas nuevas. Pero he intentado (i) basar mi respuesta en los conocimientos que pareces tener y (ii) ser lo más conciso posible. Así que puede que esta respuesta te resulte útil a ti, o a otras personas que se encuentren en una situación similar.

Explicación matemática (simple)

La SVD y la eigendecomposición están estrechamente relacionadas. Sea $X \in \mathbb{R}^{n \times p}$ sea una matriz de datos reales, por lo que se puede definir su matriz de covarianza $C \in \mathbb{R}^{p\times p}$ como \begin{equation} C = \frac{1}{n} X^T X. \end{equation}

1 | SVD de X

Como ha dicho correctamente, la aplicación de la SVD en $X$ descompone sus datos originales en \begin{equation*} X = U S V^T \end{equation*} con $U \in \mathbb{R}^{n \times n}$ , $V \in \mathbb{R}^{p \times p}$ siendo unitaria, conteniendo la (ortonormal) componentes principales y vectores propios respectivamente. La matriz diagonal $S \in \mathbb{R}^{n \times p}$ tiene el valores singulares $s$ .

2 | Eigendecomposición de C

Desde $C$ es hermitiana su eigendecomposición da como resultado vectores propios dada por la matriz unitaria $V$ con su correspondiente valores propios reales $\lambda$ como entradas de una matriz diagonal $\Lambda \in \mathbb{R}^{p \times p}$ : \begin{equation} CV = V\Lambda \end{equation} En este caso podemos calcular los componentes principales proyectando los vectores propios sobre los datos originales PCs $= X V^T$ . Obsérvese que estos PC están escalados por sus correspondientes valores propios y, por tanto, son correlacionado .

3 | SVD de C

Para responder a tus preguntas, recuerda que podemos factorizar $C$ - ya que es simétrica - a través de \begin{equation} C = V\Lambda V^T \end{equation} utilizando sus vectores y valores propios.

Obsérvese que esto también es cierto con sólo reordenar la ecuación de la sección (2). Por lo tanto, podemos calcular los vectores propios de $X$ aplicando la SVD a $C$ .

Con un poco más de esfuerzo podemos ahora establecer la relación entre los valores singulares y los valores propios. Utilizando la definición de la covarianza podemos escribir: \begin{align} C &= \frac{1}{n} ( U S V^T )^T ( U S V^T ) \\ &= \frac{1}{n} ( V S U^T U S V^T ) \\ &= \frac{1}{n} ( V S^2 V^T ) \end{align} La última ecuación se mantiene ya que $U$ es unitaria, es decir $U^T U = \mathbb{1}$ . Ahora simplemente comparando este resultado con el de arriba encontramos: \begin{equation} \frac{1}{n} ( V S^2 V^T ) = V \Lambda V^T \quad \Rightarrow \quad \lambda = \frac{s^2}{n} \end{equation}

---

Python, NumPy y algoritmos

Sólo un ejemplo básico para explorar el diferente comportamiento de numpy 's

• linalg.svd(X)
• linalg.svd(C)
• linalg.eig(C)
• linalg.eigh(C)

muestra que linalg.eig() revela un comportamiento inesperado (al menos para mí). Cálculo de la matriz $V^TV=\mathbb{1}$ para los cuatro casos, podemos hacernos una idea visual de la precisión respectiva. De la figura siguiente se desprende que linalg.eig() sólo proporciona una solución estable hasta la dimensión $d = \text{rank}(C) = \text{min}(n,p)$ .

---

# Create random data
n,p = [100,300]
X = np.random.randn(n,p)

# Covariane matrix
C = X.T @ X /n

# Create figure environment
fig = plt.figure(figsize=(14,5))
ax1 = fig.add_subplot(141)
ax2 = fig.add_subplot(142)
ax3 = fig.add_subplot(143)
ax4 = fig.add_subplot(144)

# 1. SVD on X
# ---------------------
U,s,VT = np.linalg.svd(X)
V = VT.T
ax1.imshow(V.T@V,cmap='Reds',vmin=0,vmax=1)

# 2. SVD on C
# ---------------------
V,eigenvalues,VT = np.linalg.svd(C)
ax2.imshow(V.T@V,cmap='Reds',vmin=0,vmax=1)
# 3. Eigendecomposition on C
# -> linalg.eig()
# ---------------------
eigenvalues,V = np.linalg.eig(C)
sortIdx = np.argsort(eigenvalues)[::-1]
V = V[:,sortIdx]
ax3.imshow((V.T@V).real,cmap='Reds',vmin=0,vmax=1)

# 4. Eigendecomposition on C
# -> linalg.eigh()
# ---------------------
eigenvalues,V = np.linalg.eigh(C)
sortIdx = np.argsort(eigenvalues)[::-1]
V = V[:,sortIdx]
ax4.imshow((V.T@V).real,cmap='Reds',vmin=0,vmax=1)

for a in [ax1,ax2,ax3,ax4]:
    a.set_xticks([])
    a.set_yticks([])
ax1.set_title('svd(X)')
ax2.set_title('svd(C)')
ax3.set_title('eig(C)')
ax4.set_title('eigh(C)')
fig.subplots_adjust(wspace=0,top=0.95)  
fig.suptitle('Eigendecomposition in NumPy')
plt.show()

Answer 5

Para responder a la última parte de su pregunta, "¿Por qué hacen la SVD de la matriz de covarianza y no de la matriz de datos?" Creo que es por razones de rendimiento y almacenamiento. Típicamente, $m$ será un número muy grande e incluso si $n$ es grande, esperaríamos que $m \gg n$ .

Calcular la matriz de covarianza y luego realizar la SVD en ella es mucho más rápido que calcular la SVD en la matriz de datos completa en estas condiciones, para el mismo resultado.

Incluso para valores bastante pequeños, las ganancias de rendimiento son factores de miles (milisegundos frente a segundos). Hice algunas pruebas en mi máquina para comparar usando Matlab:

Eso es sólo el tiempo de la CPU, pero las necesidades de almacenamiento son tanto o más importantes. Si intentas hacer la SVD en una matriz de un millón por mil en Matlab, dará un error por defecto, porque necesita un tamaño de matriz de trabajo de 7,4TB.