Python 中手动 SVD 实现失败

Question

我正在尝试自己在Python中实现奇异值分解（SVD）函数，使用A^T A和AA^T的特征值分解，但是重建的矩阵（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_AAT]))

# 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))

有人知道为什么这个重建只是偶尔与原始矩阵 A 匹配吗？

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

Answer 1

谢谢@Ghorban M. Tavakoly。现在我明白这是由于符号不确定性造成的。为了避免这种情况，我需要首先从 ATA（或 AAT）重建 U（或 Vt），然后使用重建的 U（或 Vt）重建另一个。这是更新的代码

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

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

# Sort eigenvaluess and eigenvectors
idx_ATA = eigenvalues_ATA.argsort()[::-1]
sorted_eigenvectors_ATA = eigenvectors_ATA[:,idx_ATA]
sorted_eigenvalues_ATA = eigenvalues_ATA[idx_ATA]

# Calculate singular values
sorted_singularvalues_ATA = np.sqrt(np.abs(sorted_eigenvalues_ATA))

# Construct diagonal matrix S and (1/S)^T
S = np.zeros_like(A)
np.fill_diagonal(S, sorted_singularvalues_ATA[:len(S)])
S_inverse_transpose = np.zeros_like(A).T
np.fill_diagonal(S_inverse_transpose, 1/sorted_singularvalues_ATA[:len(S)])

# Reconstruct U and then the original matrix as B
U = A @ sorted_eigenvectors_ATA.T @ S_inverse_transpose
B = U @ S @ sorted_eigenvectors_ATA
print(np.allclose(A,B))