SVD: The restored matrix doesn't match with the original matrix

Question

I'm trying a SVD with sympy on the matrix:

[3, 2, 2]
[2, 3, -2]

However, the restored matrix doesn't match with the original matrix.

Here's my code:

import sympy

A = sympy.Matrix([
  [3, 2, 2],
  [2, 3, -2]
])

A1 = A@A.transpose()
[U, D1] = A1.diagonalize(normalize=True)

A2 = A.transpose()@A
[V, D2] = A2.diagonalize(normalize=True)
V_T = V.transpose()

S = sympy.sqrt(D2).doit()
S = S.row_del(0)

U@S@V_T

The output is:

[2, 3, -2]
[3, 2, 2]

The 1st row and the 2nd row are swapped from the original one.
I know that U causes this result, which is:

Matrix([
[-sqrt(2)/2, sqrt(2)/2],
[ sqrt(2)/2, sqrt(2)/2]])

instead of:

[ sqrt(2)/2, sqrt(2)/2],
[-sqrt(2)/2, sqrt(2)/2]])

..., but I don't know how to fix this.

I added a parameter sort to diagonalize():

[U, D1] = A1.diagonalize(sort=True, normalize=True)

But, it makes no difference.

How would you solve this problem? Thank you in advance.

Answer 1

I was also tricked by my lecturer. I also thought that SVD was super simple and everything. But eigenvectors are not unique. Even when normalized, each column can be multiplied by a -1 to give the same answer.

The reason why your SVD does not work is that U and V are interlinked in that their column signs need to work together in some way. Doing separated diagonalizations ignores this inherent dependence and you are left with an incorrect answer most of the time.

See this Maths Exchange question and many others for similar explanations.

The following does what is is said in that Maths Exchange question.

import sympy

A = sympy.Matrix([
  [3, 2, 2],
  [2, 3, -2]
])
# A = U S V'

if A.shape[0] <= A.shape[1]:
    A1 = A * A.T
    U, S = A1.diagonalize(normalize=True)
    V_T = S**-1 * U.T * A
    print(U * S * V_T)
else:
    A2 = A.T * A
    V, S = A2.diagonalize(normalize=True)
    U = A * V * S**-1
    print(U * S * V.T)

I have added two separate cases if A is "landscape" or "portrait" just because of the way the maths plays out in each one. Hopefully it is correct.

NB: I highly advise you to use some kind of numerical computing library instead of creating your own function using a symbolic maths library. If you are not using symbols, the numerical computing library is likely to be way faster.