手动实现的Numpy SVD仅适配部分矩阵，如何稳定其效果？

问题背景 #

我用Numpy手动实现了奇异值分解（SVD），代码如下：

import numpy as np

array = np.array([[-3,3,6],
                  [3,8,7]])  # 可替换为任意矩阵，部分矩阵能正常运行，部分不行

# 左奇异向量计算
AAT = np.matmul(array, array.T)
eAAT_values, eAAT_vectors = np.linalg.eig(AAT)
idx = eAAT_values.argsort()[::-1]      # 按奇异值从大到小排序
eAAT_values = eAAT_values[idx]
eAAT_vectors = eAAT_vectors[:,idx]
SL = eAAT_vectors
# 部分矩阵需要添加这行才能得到正确结果：
# SL[:,0] = -SL[:,0]

# 右奇异向量计算
ATA = np.matmul(array.T,array)
eATA_values, eATA_vectors = np.linalg.eig(ATA)
idx = eATA_values.argsort()[::-1]     # 按奇异值从大到小排序
eATA_values = eATA_values[idx]
eATA_vectors = eATA_vectors[:,idx]
SR = eATA_vectors.T

# 奇异值矩阵
S0 = np.zeros(np.shape(array))
np.fill_diagonal(S0, np.sqrt(eATA_values), wrap=True)

# 验证重构结果
Proof = np.matmul(SL,np.matmul(S0,SR))     # 部分矩阵能匹配原矩阵，部分不行

# numpy自带SVD的对比验证
U, S, Vt = np.linalg.svd(array) 
Sm = np.zeros(np.shape(array))
np.fill_diagonal(Sm, S, wrap=True)
proof2 = np.matmul(U, np.matmul(Sm,Vt))    # 始终能正确重构

这个实现存在一个问题：部分矩阵可以直接得到正确的SVD重构结果，但有些矩阵必须对eAAT_vectors的首个特征向量做符号翻转才能生效。

• 无需翻转即可正常运行的矩阵示例：

  array = np.array([[-3,3,6], [3,8,7]])
  array = np.array([[7,-12,2], [-4,3,1]])
  

• 直接运行无法得到正确结果的矩阵示例：

  array_not_working1 = np.array([[4,5,9], [3,2,6]])
  array_not_working2 = np.array([[3,-1,4], [1,5,9]])
  

问题根源 #

SVD的解本身存在符号歧义性：对于任意合法的分解 A = UΣVᵀ，如果把U的某一列乘以-1，同时把Vᵀ对应的行也乘以-1，结果依然满足A = UΣVᵀ。

手动实现中，U（即代码里的SL）来自AAᵀ的特征向量，Vᵀ（即SR）来自AᵀA的特征向量，但numpy的eig函数返回特征向量时，符号没有统一约定——特征向量的正负都是合法解，这就导致两组特征向量的符号可能不匹配，最终重构出的矩阵和原矩阵不符。

稳定实现的两种方法 #

方法1：通过原矩阵推导匹配的奇异向量（最可靠） #

既然A = UΣVᵀ，那么可以通过其中一组奇异向量直接推导另一组，保证符号完全匹配：

• 已知V（右奇异向量）和奇异值σ，则U的第i列可通过 U[:,i] = (A @ V[:,i]) / σ_i 计算
• 反之，已知U和σ，则V的第i列可通过 V[:,i] = (Aᵀ @ U[:,i]) / σ_i 计算

这种方法完全避开了特征向量符号不统一的问题，实现代码如下：

import numpy as np

array = np.array([[4,5,9], [3,2,6]])  # 用之前不工作的矩阵测试

# 先计算右奇异向量与奇异值
ATA = np.matmul(array.T, array)
eATA_values, eATA_vectors = np.linalg.eig(ATA)
# 按奇异值从大到小排序
idx = eATA_values.argsort()[::-1]
eATA_values = eATA_values[idx]
eATA_vectors = eATA_vectors[:, idx]
SR = eATA_vectors.T  # Vᵀ
S = np.sqrt(eATA_values)
# 构建奇异值矩阵
S0 = np.zeros(np.shape(array))
np.fill_diagonal(S0, S, wrap=True)

# 通过V推导U，保证符号完全匹配
SL = np.zeros((array.shape[0], len(S)))
for i in range(len(S)):
    if S[i] > 1e-10:  # 忽略数值误差导致的极小奇异值
        SL[:, i] = (array @ eATA_vectors[:, i]) / S[i]

# 验证重构结果
Proof = np.matmul(SL, np.matmul(S0, SR))
print("重构误差：", np.linalg.norm(Proof - array))
print("重构矩阵：\n", Proof)
print("原矩阵：\n", array)

方法2：符号校验与修正 #

如果坚持用特征向量分别计算U和V，可以通过原矩阵校验两组向量的符号是否匹配，不匹配时翻转符号：
对每个非零奇异值对应的索引i，检查A @ V[:,i]和σ_i * U[:,i]的符号是否一致——若点积为负，说明符号相反，翻转U[:,i]（或V[i,:]）的符号即可。

实现代码如下：

import numpy as np

array = np.array([[4,5,9], [3,2,6]])

# 原有左奇异向量计算
AAT = np.matmul(array, array.T)
eAAT_values, eAAT_vectors = np.linalg.eig(AAT)
idx = eAAT_values.argsort()[::-1]
eAAT_values = eAAT_values[idx]
eAAT_vectors = eAAT_vectors[:, idx]
SL = eAAT_vectors.copy()

# 原有右奇异向量与奇异值计算
ATA = np.matmul(array.T, array)
eATA_values, eATA_vectors = np.linalg.eig(ATA)
idx = eATA_values.argsort()[::-1]
eATA_values = eATA_values[idx]
eATA_vectors = eATA_vectors[:, idx]
SR = eATA_vectors.T
S = np.sqrt(eATA_values)
S0 = np.zeros(np.shape(array))
np.fill_diagonal(S0, S, wrap=True)

# 校验并修正符号
for i in range(len(S)):
    if S[i] < 1e-10:
        continue
    # 计算理论上的U列向量
    expected_u = (array @ eATA_vectors[:, i]) / S[i]
    # 点积为负说明符号相反
    if np.dot(SL[:, i], expected_u) < 0:
        SL[:, i] = -SL[:, i]

# 验证重构结果
Proof = np.matmul(SL, np.matmul(S0, SR))
print("重构误差：", np.linalg.norm(Proof - array))

内容的提问来源于stack exchange，提问作者LF-137