Valores singulares de la matriz de datos y valores propios de la matriz de covarianza

Question

Tengo problemas para calcular SVD y PCA en Matlab. No sé si estoy cometiendo errores teóricos o errores de programación.

A partir de una matriz de datos $X$ PCA calcula los valores propios $\lambda_i$ de la matriz de covarianza $X^TX/(n-1)$ .

Por otro lado, la SVD de $X$ viene dada por $X=U\Sigma V^\top$ y así $$X^T X=V \Sigma^T U^T U\Sigma V^T=V \Sigma^2 V^T,$$ donde $\Sigma$ es una matriz diagonal de valores singulares con elementos $\sigma_i$ .

Así que tenemos que $$\lambda_i=\sigma_i^2 /(n-1).$$

Ahora pruebo con un ejemplo de Matlab:

    load hald;
    [u s v] = svd(ingredients);
    sigma = cov(ingredients);
    [a,b] = eig(sigma);
    disp('sigma')
    disp(diag(s)')
    disp('lambda')
    disp(diag(b)')

y aquí está la salida:

sigma
  211.3369   77.2356   28.4597   10.2667

lambda
    0.2372   12.4054   67.4964  517.7969

Los valores obtenidos no respetan la ecuación original. ¿Dónde está el error?

Answer 1

Aparte de decir lo obvio: eig ofrece los resultados en orden ascendente, mientras que svd en orden descendente; el svd los valores propios (y los vectores propios, obviamente) son distintos de los de eig descomposición porque su matriz ingredients no es simétrico para empezar. Parafraseando wikipedia un poco: " Cuando el $X$ es una matriz normal y/o semidefinida positiva, la descomposición $\ {X} = {U} {D} {U}^*$ es también una descomposición de valor singular ", no de otro modo. ( $U$ que son los vectores propios de $XX^\mathbf{T}$ )

Así, por ejemplo, si usted hizo algo como:

rng(0,'twister')        %just set the seed.
Q = random('normal', 0,1,5);
X =  Q' * Q;            %so X is PSD 
[U S V]=    svd(X);
[A,B]=      eig(X);

max( abs(diag(S)- fliplr(diag(B)')' ))
% ans =  7.1054e-15     % AKA equal to numerical precision.

usted encontraría que svd y eig le devuelven los mismos resultados. Mientras que antes exactamente porque matriz ingredients no era al menos PSD (o incluso cuadrado para el caso), bueno.. no obtuviste los mismos resultados. :)

Por decirlo de otra manera: $X= U\Sigma V^*$ prácticamente se traduce en: $X = \sum_1^r u_i s_i v_i^T$ ( $r$ siendo el rango de $X$ ). Lo que a su vez significa que se le permite (de forma bastante asombrosa) escribir $X v_i = \sigma_i u_i$ . Claro para volver a la eigen-descomposición $X u_i = \lambda_i u_i$ necesitas ante todo $u_i$ == $v_i$ . Algo que las matrices no normales no garantizan. Como nota final: Las pequeñas diferencias numéricas se deben a eig y svd con diferentes algoritmos trabajando en segundo plano; una variante de la Algoritmo QR de svd y un (normalmente) descomposición de Schur generalizada para eig .

Específicamente para tu problema lo que quieres es algo parecido a:

load hald;
[u s v]=svd(ingredients);
sigma=(ingredients' * ingredients); 
lambda =eig(sigma);     
max( abs(diag(s)- fliplr(sqrt(lambda)')' ))
% ans = 5.6843e-14

Como ves esto no tiene nada que ver con centrar tus datos para tener media $0$ en este punto la matriz ingredients no está centrado.

Ahora bien, si utilizas la matriz de covarianza (y no una simple matriz de producto interior como hice yo) tienes tendrá que centrarse sus datos. Digamos que ingredients2 es su muestra de media cero.

ingredients2 = ingredients - repmat(mean(ingredients), 13,1);

Entonces sí que se necesita esta normalización por $1/(n-1)$

[u s v] =svd(ingredients2 );        
sigma = cov(ingredients); % You don't care about centring here
lambda =eig(sigma);   

max( abs( diag(s)- fliplr(sqrt(lambda *12)')')) % n = 13 so multiply by n-1
% ans = 4.7962e-14

Así que sí, es el centrado ahora . Me equivoqué un poco al principio porque trabajaba con la noción de matrices PSD en lugar de matrices de covarianza. La respuesta antes de la edición estaba bien. Abordaba exactamente por qué tu descomposición de valores propios no se ajustaba a tu descomposición de valores singulares. Con la edición que mostrar por qué su descomposición de valor singular no se ajustaba a la eigen-descomposición. Es evidente que se puede ver el mismo problema de dos maneras diferentes :D