CODEWITHSUNDEEP
pythonnumpymathsvd
LaLeLo
LaLeLo

Reputation: 157

Why does my own simple Implemantion of the svd algorithm in python not work?

I implemented a very simple singular value decompositon in Python (don't worry i'm aware of the numpy version), it is more for fun. It is based on the computation of the eigenvectors of A^TA and AA^T (see for example here). This is in theory very simple to implement in Python:

   def calc_svd(A: np.array):
    ATA = A.T@A
    AAT = [email protected]
    #compute the eigenvectors and sort them in ascending size of the eigenvalues
    eigenvalues_ATA, eigenvectors_ATA = np.linalg.eig(ATA)
    ATA_sorted = np.flip(np.argsort(eigenvalues_ATA))
    eigenvectors_ATA = eigenvectors_ATA[:, ATA_sorted]

    eigenvalues_AAT, eigenvectors_AAT = np.linalg.eig(AAT)
    AAT_sorted = np.flip(np.argsort(eigenvalues_AAT))
    eigenvalues_AAT = eigenvalues_AAT[AAT_sorted]
    eigenvectors_AAT = eigenvectors_AAT[:, AAT_sorted]
    sing_values = np.sqrt(eigenvalues_AAT)

    return U, sing_values, V.T

As far as I see this this should return the correct matrices. However there is a problem with this:

test = np.array([[1.0, 1.0, 0.0], [1.0, 3.0, 1.0], [2.0, -1.0, 1.0]])
u_np, svd_np, v_np = np.linalg.svd(test2)
u_own, svd_own, v_own = calc_svd(test2)
print([email protected](svd_np)@v_np)
print([email protected](svd_own)@v_own)

results in the following output:

[[ 1.00000000e+00  1.00000000e+00  5.76517542e-16]
 [ 1.00000000e+00  3.00000000e+00  1.00000000e+00]
 [ 2.00000000e+00 -1.00000000e+00  1.00000000e+00]]
[[-0.53412934 -0.91106401 -0.94056803]
 [-1.17456455 -3.03332485 -0.64756347]
 [-2.08209125  0.98432856 -0.83426215]]

The problem seems to be that some eigenvectors of U and V in my implemention are the negative of the one of numpy:

print(u_np)
print(u_own)
#results in
[[-0.35881632  0.13270783 -0.92392612]
 [-0.9317945  -0.10910581  0.34620071]
 [-0.05486216  0.98513174  0.16280537]]
[[ 0.35881632  0.13270783 -0.92392612]
 [ 0.9317945  -0.10910581  0.34620071]
 [ 0.05486216  0.98513174  0.16280537]]
#and
print(v_np)
print(v_own)
#in
[[-0.39543728 -0.87521453 -0.27861961]
 [ 0.80500544 -0.47631006  0.35368767]
 [-0.44226191 -0.08442901  0.89290321]]
[[-0.39543728 -0.87521453 -0.27861961]
 [-0.80500544  0.47631006 -0.35368767]
 [-0.44226191 -0.08442901  0.89290321]]

My problem is that I don't understand why this is not working. From a mathematical point of view it shouldn't make a difference whether to take the eigenvector v or -v since they are both eigenvectors to the same eigenvalue. So can maybe someone tell my error in my reasoning.

Thank you very much in advance :)

Upvotes: 0

Views: 517

Answers (1)

Max Kasper
Max Kasper

Reputation: 38

The approach is correct in itself, but you have a problem because the eigenvectors of U and V are not uniquely determined with respect to multiplication by -1.

It is important that the identity

A * v_i = σ_i * u_i

holds for the columns v_i, u_i and the singular values u_i.

Upvotes: 1

Related Questions