CODEWITHSUNDEEP
pythonsvd
Jihyun
Jihyun

Reputation: 1105

Manual SVD Implementation in Python failed

I am trying to implement a function for Singular Value Decomposition (SVD) in Python by myself, using the eigenvalue decomposition of A^T A and AA^T, but the reconstructed matrix (B) does not always match the original matrix (A). Here is my code:

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

Could anyone know why this reconstruction only occasionally match the original matrix 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]]

Upvotes: 5

Views: 152

Answers (1)

Jihyun
Jihyun

Reputation: 1105

Thanks to @Ghorban M. Tavakoly, now I understand it's due to sign indeterminacy. To avoid it, I need to reconstruct U (or Vt) first from ATA (or AAT) and then the other one using the reconstructed U (or Vt). Here is the updated code:

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

Upvotes: 2

Related Questions