Why is my SVD calculation different than numpy's SVD calculation of this matrix?

I'm trying to manually compute the SVD of the matrix A defined below but I am having some problems. Computing it manually and with the svd method in numpy yields two different results.

Computed manually below:

import numpy as np
A = np.array([[3,2,2], [2,3,-2]])
V = np.linalg.eig(A.T @ A)[1]
U = np.linalg.eig(A @ A.T)[1]
S = np.c_[np.diag(np.sqrt(np.linalg.eig(A @ A.T)[0])), [0,0]]
print(A)
print(U @ S @ V.T)

And computed via numpy's svd method:

X,Y,Z = np.linalg.svd(A)
Y = np.c_[np.diag(Y), [0,0]]
print(A)
print(X @ Y @ Z)

When these two codes are run. The manual calculation doesn't equal the svd method. Why is there a discrepancy between these two calculations?

Solution

Look at the eigenvalues returned by np.linalg.eig(A.T @ A):

In [57]: evals, evecs = np.linalg.eig(A.T @ A)

In [58]: evals
Out[58]: array([2.50000000e+01, 3.61082692e-15, 9.00000000e+00])

So (ignoring the normal floating point imprecision), it computed [25, 0, 9]. The eigenvectors associated with those eigenvalues are in the columns of evecs, in the same order. But your construction of S doesn't match that order; here's your S:

In [60]: S
Out[60]: 
array([[5., 0., 0.],
       [0., 3., 0.]])

When you compute U @ S @ V.T, the values in S @ V.T are not correctly aligned.

As a quick fix, you can rerun your code with S set explicitly as follows:

S = np.array([[5, 0, 0],
              [0, 0, 3]])

With that change, your code outputs

[[ 3  2  2]
 [ 2  3 -2]]
[[-3. -2. -2.]
 [-2. -3.  2.]]

That's better, but why are the signs wrong? Now the problem is that you have independently computed U and V. Eigenvectors are not unique; they are the basis of an eigenspace, and such a basis is not unique. If the eigenvalue is simple, and if the vector is normalized to have length one (which numpy.linalg.eig does), there is still a choice of the sign to be made. That is, if v is an eigenvector, then so is -v. The choices made by eig when computing U and V won't necessarily result in restoring the sign of A when U @ S @ V.T is computed.

It turns out that you can get the result that you expect by simply reversing all the signs in either U or V. Here is a modified version of your script that generates the output that you expected:

import numpy as np

A = np.array([[3,  2,  2],
              [2,  3, -2]])

U = np.linalg.eig(A @ A.T)[1]
V = -np.linalg.eig(A.T @ A)[1]
#S = np.c_[np.diag(np.sqrt(np.linalg.eig(A @ A.T)[0])), [0,0]]
S = np.array([[5, 0, 0],
              [0, 0, 3]])

print(A)
print(U @ S @ V.T)

Output:

[[ 3  2  2]
 [ 2  3 -2]]
[[ 3.  2.  2.]
 [ 2.  3. -2.]]