自定义SVD实现无法还原原矩阵，与numpy官方实现存在符号差异

我自己实现了一个SVD算法，但有时候无法还原原始矩阵。以下是我的实现代码：

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)

我发现一个奇怪的现象：自己实现得到的l（即V^T）和numpy自带SVD返回的z，绝对值基本一致，但始终存在符号差异。虽然已经对特征向量做了排序，l和z的转置乘积是近似对角矩阵，但对角线上的符号时而正时而负，导致偶尔还原出的矩阵和原始矩阵符号完全相反。

输出示例1 #

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

输出示例2 #

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

问题根源 #

特征向量的符号是不唯一的：np.linalg.eig返回的特征向量，其方向（正负）是随机的。你的实现中，U的列是基于A@A.T的特征向量排序，V的列是基于A.T@A的特征向量排序，两者的符号没有同步校准。当U的某列和对应的V的列符号相反时，U @ diag(sigma) @ V.T就会出现符号翻转，导致还原矩阵错误。

解决方法 #

需要根据原矩阵A来校准U和V的符号，确保对于每个奇异值sigma[i]，满足U[:,i].T @ A @ V[:,i] = sigma[i]（因为奇异值都是非负的）。修改后的代码如下：

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_sigma = np.argsort(-sigma)[:num]
    sigma = sigma[sorted_indices_sigma]

    # 对U和V的特征向量按奇异值对应的顺序排序（A@A.T和A.T@A的非零特征值相同）
    sorted_indices_left = np.argsort(-eigvals_left)[:num]
    U = eigvecs_left[:, sorted_indices_left]
    sorted_indices_right = np.argsort(-eigvals_right)[:num]
    V = eigvecs_right[:, sorted_indices_right]

    # 校准U和V的符号
    for i in range(num):
        val = U[:, i].T @ A @ V[:, i]
        if val < 0:
            U[:, i] = -U[:, i]
            V[:, i] = -V[:, i]

    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)

这样修改后，U和V的符号会同步，还原出的矩阵就能和原始矩阵一致，符号差异问题也会解决。

内容的提问来源于stack exchange，提问作者chengyongru