# Computing SVD through Eigendecomposition of correlation matrix

I am following the excellent series on SVD by Steve Brunton from the University of Washington, on YouTube, but I have trouble interpreting his 4th video on the subject.

If I understand correctly, he mentions that one can compute the economy SVD decomposition $$X = \hat{U}\hat{\Sigma}V^T$$ with the following :

$$X^TX= V\hat{\Sigma}\hat{U}^T\hat{U}\hat{\Sigma}V^T = V\hat{\Sigma}^2V^T \implies X^TXV = V\hat{\Sigma}^2$$ $$XX^T= \hat{U}\hat{\Sigma}\hat{V}^T\hat{V}\hat{\Sigma}\hat{U}^T = \hat{U}\hat{\Sigma}^2\hat{U}^T \implies XX^T\hat{U} = \hat{U}\hat{\Sigma}^2$$

Instead of using U,S,VT = svd(X) in Python, I want to try decomposing the image in its SVD, and reconstructing the original using only $$r$$ values with this technique. I am trying to apply this on a picture of The Starry Night, loaded in a grayscale numpy array X, which I do this way :

r = 10
XT = X.transpose()
C_1 = XT @ X
C_2 = X @ XT
[Lambda_1, V_hat] = np.linalg.eig(C_1)
[Lambda_2, U_hat] = np.linalg.eig(C_2)
V_hat_T = V_hat.transpose()
X_tilde = U_hat[:,:r] @ np.diag(Lambda_1)[0:r,:r] @ V_hat_T[:r,:]


But I can't reproduce the image, whatever the value I attribute to $$r$$, whereas the product of U @ S @ VT , computed with the svd(X) function, perfectly reconstructs the original picture.

What am I doing wrong? It is certainly superfluous to mention that I am a beginner, and I probably made some big mistake.

• Please share more about your image tensors if possible. And try X_tilde = U_hat[:,:r] @ np.diag(Lambda_2)[0:r,:r] @ V_hat_T[:r,:] at the same time. Apr 24 at 23:14
• In this case I simply use the following code : A = Image.open('vangogh.jpg') X = np.mean(A, -1) # grayscale Running the X_tilde line with Lambda_1 or Lambda_2 does not produce any meaningful results... Apr 27 at 4:02

The problem is occurring for a few reasons:

• You must ensure the correct ordering of eigenvalues/eigenvectors.
• You should be taking the square root of the eigenvalues (Lambda_1) to get the singular values.
• You need to take special care to ensure you are using the right signs for the eigenvectors as these won't necessarily be correct by default. This is because we assume that the singular values are the positive roots of the eigenvalues, but the two sets of eigenvectors we retrieve don't necessarily respect this assumption when joined in the singular value decomposition.

Here is the working example:

r = 10
XT = X.transpose()
C_1 = XT @ X
C_2 = X @ XT
[Lambda_1, V_hat] = np.linalg.eig(C_1)
[Lambda_2, U_hat] = np.linalg.eig(C_2)

# The eigenvalues/eigenvectors returned from np.linalg.eig are not sorted so may differ between Lambda_1 and Lambda_2
# When performing economy SVD we also want to ensure the largest eigenvalues come first
i_1 = np.argsort(Lambda_1)[::-1]
i_2 = np.argsort(Lambda_2)[::-1]
V_hat = V_hat[:, i_1]
U_hat = U_hat[:, i_2]

# We must take the square root of eigenvalues to get the singular values
Lambda = np.sqrt(np.sort(Lambda_1)[::-1])

# X @ V_hat and U_hat @ np.diag(Lambda) should be equal. But the eigenvectors retrieved may not have the correct signs
# This code checks if  X @ V_hat and U_hat @ np.diag(Lambda) have different signs and then updates V_hat wherever the sign is incorrect
same_sign = np.sign((X @ V_hat)[0] * (U_hat @ np.diag(Lambda))[0])
V_hat = V_hat * same_sign.reshape(1, -1)

V_hat_T = V_hat.transpose()
X_tilde = U_hat[:,:r] @ np.diag(Lambda)[0:r,:r] @ V_hat_T[:r,:]

• Absolutely. I forgot to mention it, but I do check for ordering and square root, however I did not know about the different signs, thank you for the input. Nonetheless, I have trouble with the dimensions of matrices when running your code...Also, why aren't Lambda_1 and Lambda_2 the same? Apr 27 at 4:01
• To your first point: Are you running this with a square image? My guess is you are using a non-square X which may lead to the dimension issues as well as other issues - np.linalg.eig is not really designed for non-square matrices. To your second point: If you sort Lambda_1 and Lambda_2, they should be the same (or very close, there will be small differences due to numerical issues). But as returned from np.linalg.eig the eigenvalues/eigenvectors are not arranged in any particular order so you are not guaranteed to have the same ordering over eigenvalues in Lambda_1 and Lambda_2. Apr 27 at 5:23
• Good point! Actually, I tried first using a square image, but to keep track of operations, I used a rectanglar one to be sure I understand every step, and now I completely missed that eigendecomposition is indeed not designed for non-square matrices. And now it is working perfectly fine! Thank you so much for helping me learn. I will mark your post as an answer. Apr 30 at 21:10