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]]
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