Python手动实现SVD时矩阵重构偶发失败的原因排查

我尝试通过对A^T A和AA^T进行特征值分解，在Python中手动实现奇异值分解（SVD）函数，但重构矩阵B并非总能与原矩阵A匹配。以下是我的实现代码：

import numpy as np
# Generate a random matrix A
row, col = 3, 3
A = np.random.normal(size=row * col).reshape(row, col)

# Eigen decomposition of A^T*A and A*A^T
ATA = A.T @ A
AAT = A @ A.T
eigenvalues_ATA, eigenvectors_ATA = np.linalg.eig(ATA)
eigenvalues_AAT, eigenvectors_AAT = np.linalg.eig(AAT)

# Sort eigenvalues and eigenvectors
idx_ATA = eigenvalues_ATA.argsort()[::-1]
idx_AAT = eigenvalues_AAT.argsort()[::-1]
sorted_eigenvectors_ATA = eigenvectors_ATA[:, idx_ATA]
sorted_eigenvectors_AAT = eigenvectors_AAT[:, idx_AAT]

# Calculate singular values
sorted_singularvalues_ATA = np.sqrt(np.abs(eigenvalues_ATA[idx_ATA]))
sorted_singularvalues_AAT = np.sqrt(np.abs(eigenvalues_AAT[idx_ATA]))

# Construct diagonal matrix S
S = np.zeros_like(A)
np.fill_diagonal(S, sorted_singularvalues_ATA)

# Reconstruct matrix B
B = sorted_eigenvectors_AAT @ S @ sorted_eigenvectors_ATA.T
print(np.allclose(A, B))

同时提供了重构成功与失败的矩阵示例：

# Example it works
A_equal = [-1.59038869, -0.28431377,  0.36309318,
           0.07133563, -0.20420962,  1.82207923,
           0.84681193,  0.31419994, -0.93808105]

# Example it fails
A_not_equal = [ 1.61171729,  0.6436384,   0.47359656,
               -1.04121454,  0.17558459,  0.36595138,
                0.40957221,  0.20499528,  0.18525562]
A = np.array(A_not_equal).reshape(3,3)

# Expected output
[[ 1.61171729  0.6436384   0.47359656]
 [-1.04121454  0.17558459  0.36595138]
 [ 0.40957221  0.20499528  0.18525562]]

# Actual output
[[-1.61240387 -0.63872607 -0.47789066]
 [ 1.04102391 -0.17422069 -0.36714363]
 [-0.40734839 -0.22090623 -0.17134711]]

请问为何该重构仅偶尔匹配原矩阵A？

核心问题 #

1. 特征向量的符号不确定性：np.linalg.eig返回的特征向量符号是随机的——对于同一个特征值，v和-v都是合法的特征向量。当A^T A和AA^T的对应特征向量符号相反时，重构就会出现符号翻转，导致结果与原矩阵不匹配。
2. 特征向量的对应关系未对齐：单纯对A^T A和AA^T的特征向量分别排序，无法保证U（来自AA^T的特征向量）和V（来自A^T A的特征向量）的列满足SVD的核心关系：A @ V[:,i] = σ_i * U[:,i]。

修正方案 #

需要验证并调整U的列向量符号，确保每个U[:,i]与A @ V[:,i]方向一致：

import numpy as np

row, col = 3, 3
# 用失败示例测试
A_not_equal = [ 1.61171729,  0.6436384,   0.47359656,
               -1.04121454,  0.17558459,  0.36595138,
                0.40957221,  0.20499528,  0.18525562]
A = np.array(A_not_equal).reshape(row, col)

ATA = A.T @ A
AAT = A @ A.T

# 特征值分解
eigenvalues_ATA, V = np.linalg.eig(ATA)
eigenvalues_AAT, U = np.linalg.eig(AAT)

# 按奇异值降序排序（先处理V和奇异值）
idx_ATA = np.argsort(np.sqrt(np.abs(eigenvalues_ATA)))[::-1]
sorted_singular_values = np.sqrt(np.abs(eigenvalues_ATA[idx_ATA]))
V = V[:, idx_ATA]

# 对U排序，确保和V的奇异值顺序对应
idx_AAT = np.argsort(np.sqrt(np.abs(eigenvalues_AAT)))[::-1]
U = U[:, idx_AAT]

# 调整U的符号，保证A@V[:,i]和σ_i*U[:,i]方向一致
for i in range(len(sorted_singular_values)):
    expected = A @ V[:, i]
    actual = sorted_singular_values[i] * U[:, i]
    # 计算点积判断方向，点积为负则翻转符号
    if np.dot(expected, actual) < 0:
        U[:, i] = -U[:, i]

# 构建S矩阵
S = np.zeros_like(A)
np.fill_diagonal(S, sorted_singular_values)

# 重构矩阵
B = U @ S @ V.T
print(np.allclose(A, B))  # 输出True
print(B)

额外说明 #

• 对于非方阵，还需要处理U或V的维度，但核心逻辑一致：确保U的列与A@V[:,i]的方向匹配。
• 实际工程中优先使用np.linalg.svd，它已经处理了这些细节，稳定性更高。

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