5
$\begingroup$

I am having some problem in computing SVD and PCA in Matlab. I do not know if I am doing theoretical mistakes or programming mistakes.

Starting with a data matrix $X$, PCA computes the eigenvalues $\lambda_i$ of the covariance matrix $X^TX/(n-1)$.

On the other side SVD of $X$ is given by $X=U\Sigma V^\top$, and so $$X^T X=V \Sigma^T U^T U\Sigma V^T=V \Sigma^2 V^T,$$ where $\Sigma$ is diagonal matrix of singular values with elements $\sigma_i$.

So we have that $$\lambda_i=\sigma_i^2 /(n-1).$$

Now I try a Matlab example:

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

and here is the output:

sigma
  211.3369   77.2356   28.4597   10.2667

lambda
    0.2372   12.4054   67.4964  517.7969

The obtained values do not respect the original equation. Where is the mistake?

$\endgroup$
  • 2
    $\begingroup$ Is $\mathbf X$ centered? $\endgroup$ – cbeleites Oct 8 '13 at 12:13
  • $\begingroup$ @Donbeo, welcome. Please know that it is a good practice to show the data whenever possible. I hope you don't expect people here to have your ingredients dataset right at hand? $\endgroup$ – ttnphns Oct 8 '13 at 13:16
  • 3
    $\begingroup$ @ttnphns: Totally agreed; just given that the OP asks about Matlab and that the hald dataset is part of Matlab's "Sample Data Sets" mathworks.co.uk/help/stats/_bq9uxn4.html it is an honest mistake by a new user. It is available here: qsar.org/resource/datasets/hald.htm For practical purposes ingredients is a $13\times4$ full rank matrix of positive integer values. $\endgroup$ – usεr11852 Oct 8 '13 at 13:41
  • 1
    $\begingroup$ Donbeo, a short version of the @usεr11852's answer is that in your code the $X$ matrix (ingredients) was not centered, and that's why the computations did not match. The important thing to understand is that covariance matrix is given by $X^\top X/(n-1)$ only if $X$ is centered. $\endgroup$ – amoeba Jan 19 '15 at 14:04
  • 2
    $\begingroup$ thanks a lot @amoeba. I did this question so long time ago. Indeed the problem was that the data matrix was not centred. $\endgroup$ – Donbeo Jan 19 '15 at 15:15
6
$\begingroup$

Aside stating the obvious: eig gives the results in ascending order while svd in descending one; the svd eigenvalues (and eigenvectors obviously) are dissimilar to those of eig decomposition because your matrix ingredients is not symmetric to start with. To paraphrase wikipedia a bit: "When the $X$ is a normal and/or a positive semi-definite matrix, the decomposition $\ {X} = {U} {D} {U}^*$ is also a singular value decomposition", not otherwise. ($U$ being the eigenvectors of $XX^\mathbf{T}$)

So example if you did something like:

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.

you would find that svd and eig do give you back the same results. While before exactly because matrix ingredients was not at least PSD (or even square for that matter), well.. you didn't get the same results. :)

Just to state it in another way: $X= U\Sigma V^*$ practically translates into: $X = \sum_1^r u_i s_i v_i^T$ ($r$ being the rank of $X$). Which itself means that you are (pretty awesomely) allowed to write $X v_i = \sigma_i u_i$. Clear to get back to the eigen-decomposition $X u_i = \lambda_i u_i$ you need first all $u_i$ == $v_i$. Something that non-normal matrices do not guarantee. As final note: The small numerical differences are due to eig and svd having different algorithms working in the background; a variant of the QR algorithm of svd and a (usually) generalized Schur decomposition for eig.

Specific to your problem what you want is something akin to:

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

As you see this is nothing to do with centring you data to have mean $0$ at this point; the matrix ingredients is not centered.

Now if you use the covariance matrix (and not a simple inner product matrix as I did) you will have to centre your data. Let's say that ingredients2 is your zero-meaned sample.

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

Then indeed you need this normalization by $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

So yeah, it the centring now. I was a bit misleading originally because I worked with the notion of PSD matrices rather than covariance matrices. The answer before the editing was fine. It addressed exactly why your eigen-decomposition did not fit your singular value decomposition. With the editting I show why your singular value decomposition did not fit the eigen-decomposition. Clearly one can view the same problem in two different ways. :D

$\endgroup$
  • $\begingroup$ But when we talk about PCA are we looking at the value of $\sigma_i$ or at the value of $\lambda_i$ ? $\endgroup$ – Donbeo Oct 8 '13 at 14:31
  • $\begingroup$ @Donbeo PCA uses $\lambda$. You are getting the results for : $X v= \lambda {v}$ with PCA (eig). (The actual symbol used makes little difference) $\endgroup$ – usεr11852 Oct 8 '13 at 14:38
  • $\begingroup$ Yes sure of course. I was wondering in case we decide to do variable selection removing the smaller eigenvalue we should refer to the value computed using the function eig and not to the one coputed using the function svd. Is that right? $\endgroup$ – Donbeo Oct 8 '13 at 14:46
  • $\begingroup$ @Donbeo: Variable selection has nothing to do with this. If you mean reducing your sample's dimensions then using PCA would be the standard way to do it. (so yeah exclude the dimensions related with the $k$-smallest $\lambda$'s) PCA and SVD give comparable results; they are not totally the same thing though except in special cases (as explained above). :) $\endgroup$ – usεr11852 Oct 8 '13 at 14:53
  • $\begingroup$ Glad I could help. $\endgroup$ – usεr11852 Oct 8 '13 at 15:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.