Estoy estudiando PCA del curso Coursera de Andrew Ng y otros materiales. En el curso de Stanford NLP, la primera asignación de cs224n , y en el video de la conferencia de Andrew Ng , hacen una descomposición de valores singulares en lugar de la descomposición de vectores propios de la matriz de covarianza, y Ng incluso dice que SVD es numéricamente más estable que la descomposición propia.

Según tengo entendido, para PCA deberíamos hacer SVD de la matriz de (m,n)tamaño de datos , no de la matriz de (n,n)tamaño de covarianza . Y la descomposición del vector propio de la matriz de covarianza.

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

ameba ya dio una buena respuesta en los comentarios, pero si quieres una discusión formal, aquí va.

$A$$A=U\Sigma V^T$$V$$A^TA$$\Sigma$$\sigma_{ii}=\sqrt{\lambda_i(A^TA)}$

$\frac{1}{n-1}A^TA$$\lambda_i(\frac{1}{n-1}A^TA)$

$B$$\alpha \in \mathbb R$$v$$Bv=\lambda v$

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

$S=\frac{1}{n-1}A^TA$$S$$S^TS=\frac{1}{(n-1)^2}A^TAA^TA$

1. $(A^TA)^TA^TA=A^TAA^TA$$A^TA$
2. $\frac{1}{(n-1)^2}A^TAA^TA$$\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à!

Con respecto a la estabilidad numérica, uno necesitaría descubrir cuáles son los algoritmos empleados. Si estás preparado, creo que estas son las rutinas LAPACK utilizadas por numpy:

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

Actualización: en cuanto a la estabilidad, la implementación de SVD parece estar utilizando un enfoque de divide y vencerás, mientras que la descomposición propia usa un algoritmo QR simple. No puedo acceder a algunos documentos relevantes de SIAM de mi institución (culpar a los recortes de investigación) pero encontré 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 de división y conquista espectrales estables y eficientes para la descomposición simétrica del valor propio 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 (¡usan el mismo como numpy en uno de los experimentos!) es más estable que el algoritmo QR. Esto, junto con las afirmaciones de que los métodos de D&C son más estables, respalda la elección de Ng.

@amoeba tenía excelentes respuestas a las preguntas de PCA, incluida esta en relación con SVD a PCA. Respondiendo a su pregunta exacta, haré tres puntos:

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

Resulta que 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 terminar con regresores altamente colineales. SVD funciona mejor en estos casos.

Aquí está el código de Python para demostrar el punto. Creé una matriz de datos altamente colineal, obtuve su matriz de covarianza y traté de obtener los valores propios de esta última. SVD todavía funciona, mientras que la descomposición del eigen ordinario 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)

Salida:

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

Actualizar

En respuesta al comentario de Federico Poloni, aquí está el código con pruebas 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í. La SVD es aproximadamente dos veces más precisa en una determinación de valor de eigen 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)

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

$$x_1=u\\ x_2=u+\varepsilon v$$ $u,v$ $$\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}$$ $\sigma_1^2,\sigma_2^2,\rho$

$$\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)$$ $\varepsilon$ $$\lambda\approx \sigma_2^2 \varepsilon^2 (1-\rho^2)/2$$

$j=1,\dots,m$$\hat\lambda_j$$e_j=\lambda-\hat\lambda_j$

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

eighes 9-10 veces más rápido que eigen mi computadora (independientemente del tamaño de la matriz) y tiene una mejor precisión (según la prueba de precisión de @ Aksakal).

No estoy convencido con la demostración del beneficio de precisión de SVD con valores propios pequeños. La prueba de @ Aksakal es 1-2 órdenes de magnitud más sensibles 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 sobre la precisión que la elección de un algoritmo de descomposición propia. 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 usar el método de matriz de covarianza para PCA, lo descompondría con en eighlugar de svd. Si falla (que aún no se ha demostrado aquí), entonces probablemente valga la pena repensar el problema que está tratando de resolver antes de comenzar a buscar un algoritmo mejor.

$m$$n$$m \gg n$

Calcular la matriz de covarianza y luego realizar SVD en eso es mucho más rápido que calcular 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). Realicé algunas pruebas en mi máquina para comparar usando Matlab:

Eso es solo tiempo de CPU, pero las necesidades de almacenamiento son tan importantes, si no más. Si intenta SVD en una matriz de un millón por mil en Matlab, se producirá un error de forma predeterminada, ya que necesita un tamaño de matriz de trabajo de 7,4 TB.