[FIXED] I implemented an SVD, but sometimes it is not possible to restore the original matrix

Issue

Here is my code

import numpy as np


def svd(A):
    eigvals_left, eigvecs_left = np.linalg.eig(A @ A.T)

    eigvals_right, eigvecs_right = np.linalg.eig(A.T @ A)

    sigma = np.sqrt(np.abs(eigvals_right))

    num = min(A.shape)

    sorted_indices = np.argsort(-sigma)
    sigma = sigma[sorted_indices[:num]]

    U = eigvecs_left[:, np.argsort(-eigvals_left)[:num]]
    V = eigvecs_right[:, np.argsort(-eigvals_right)[:num]]

    return U, sigma, V.T


h = np.random.rand(2, 3)

print(h)

j, k, l = svd(h)
x, y, z = np.linalg.svd(h, compute_uv=True, full_matrices=False)

print('---------------')
print(j @ np.diag(k) @ l)
print('---------------')
print(x @ np.diag(y) @ z)
print('---------------')

print(l @ z.T)

I am confused.In most cases, the absolute values of l and z are similar, but there is always a sign difference between 1 and z.I have sorted eigvecs correctly. The matrix obtained by multiplying the transpose of l and z is also something similar to a diagonal matrix, but the sign is sometimes positive and sometimes negative.

output1:

[[0.53221807 0.59786549 0.09906266]
 [0.4031512  0.38389025 0.16900433]]
---------------
[[-0.55253742 -0.55534658 -0.19184685]
 [-0.37481911 -0.44317609 -0.03963159]]
---------------
[[0.53221807 0.59786549 0.09906266]
 [0.4031512  0.38389025 0.16900433]]
---------------
[[ 1.00000000e+00  7.40766293e-17]
 [-4.02870278e-17  1.00000000e+00]]

output2:

[[0.23044744 0.86795113 0.64293873]
 [0.80871952 0.17031212 0.33708637]]
---------------
[[0.23044744 0.86795113 0.64293873]
 [0.80871952 0.17031212 0.33708637]]
---------------
[[0.23044744 0.86795113 0.64293873]
 [0.80871952 0.17031212 0.33708637]]
---------------
[[-1.00000000e+00 -2.56188682e-16]
 [-4.53777597e-16  1.00000000e+00]]

Solution

The sign (and more generally the phase) of the eigenvectors is arbitrary, two different algorithms may output different signs. In the case of an SVD, the signs of the singular vectors is a gauge freedom in the SVD definition.

If you have an SVD M = U @ S @ V, where S is diagonal and U and V are unitaries, you any sign change on U columns and V rows also gives an SVD.

n = 10
M = np.random.random((n, n))
U, s, V = np.linalg.svd(M)
S = np.diag(s)
print(np.allclose(M, U @ S @ V))
random_signs = 1 - 2 * (np.random.random(n) > 0.5)
D = np.diag(random_signs)
U2 = U @ D
V2 = D @ V
print(np.allclose(M, U2 @ S @ V2))

In the absence of degeneracy in the singular values, this is the only degree of freedom in an SVD. Some algorithm have been proposed to fix this gauge.

For more details on SVD gauge, see e.g. https://www.sciencedirect.com/science/article/pii/037704278890386X?via%3Dihub

Answered By - Olivier Gauthé