-1
\$\begingroup\$

I use the following code for a manual single value decomposition using numpy. Depending on the array I choose it sometimes works out well and I can verify the svd and sometimes it does not work out straight aways and requires a sign flip.

Code

import numpy as np

array = np.array([[-3,3,6],
                  [3,8,7]])  # here you can choose any matrix values. Some matrices work, some don't.

# left singular vectors
AAT = np.matmul(array, array.T)
eAAT_values, eAAT_vectors = np.linalg.eig(AAT)
idx = eAAT_values.argsort()[::-1]      #sorting: largest singular values first
eAAT_values = eAAT_values[idx]
eAAT_vectors = eAAT_vectors[:,idx]
SL = eAAT_vectors
# Some arrays require an additional line:
# SL[:,0] = -SL[:,0]


# right singular vectors
ATA = np.matmul(array.T,array)
eATA_values, eATA_vectors = np.linalg.eig(ATA)
idx = eATA_values.argsort()[::-1]     #sorting: largest singular values first
eATA_values = eATA_values[idx]
eATA_vectors = eATA_vectors[:,idx]
SR = eATA_vectors.T


# singular values
S0 = np.zeros(np.shape(array))
np.fill_diagonal(S0, np.sqrt(eATA_values), wrap=True)

# verifying. Expected proof == array
Proof = np.matmul(SL,np.matmul(S0,SR))     # works out with some arrays with others it does not

# Python linalg.svd works consistently
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))    # always works out

The array in the code above works well. Contrary examples where the svd does not work out without an additional sign flip are presented here.

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

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

How can I stabilize the outcome or determine when a signflip of the first eigenvector from eAAT_vectors is necessary?

Another working array without an additional sign flip:

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

\$\endgroup\$

1 Answer 1

1
\$\begingroup\$

The sign flip arises because the order of the eigenvectors returned by np.linalg.eig may not always match the expected order. To fix this you should:

Compute the singular value decomposition (SVD) using your code. Check the consistency of the result by comparing the sign of the first element of SL[:, 0] with the sign of the first element of U[:, 0] from Python's np.linalg.svd.

Here's how you can go about it on your script:

import numpy as np

def manual_svd(array):
    # left singular vectors
    AAT = np.matmul(array, array.T)
    eAAT_values, eAAT_vectors = np.linalg.eig(AAT)
    idx = eAAT_values.argsort()[::-1]  # sorting: largest singular values first
    eAAT_values = eAAT_values[idx]
    eAAT_vectors = eAAT_vectors[:, idx]
    SL = eAAT_vectors

    # Check the sign of the first element
    if np.sign(SL[0, 0]) != np.sign(array[0, 0]):
        SL[:, 0] = -SL[:, 0]  # Apply sign flip if necessary

    # right singular vectors
    ATA = np.matmul(array.T, array)
    eATA_values, eATA_vectors = np.linalg.eig(ATA)
    idx = eATA_values.argsort()[::-1]  # sorting: largest singular values first
    eATA_values = eATA_values[idx]
    eATA_vectors = eATA_vectors[:, idx]
    SR = eATA_vectors.T

    # singular values
    S0 = np.zeros(np.shape(array))
    np.fill_diagonal(S0, np.sqrt(eATA_values), wrap=True)

    # Verify with the sign-corrected SL
    Proof = np.matmul(SL, np.matmul(S0, SR))

    return SL, S0, SR, Proof

# Test your code with different arrays
array1 = np.array([[4, 5, 9], [3, 2, 6]])
array2 = np.array([[3, -1, 4], [1, 5, 9]])

SL1, S0_1, SR1, Proof1 = manual_svd(array1)
SL2, S0_2, SR2, Proof2 = manual_svd(array2)

print("Proof1:")
print(Proof1)
print("Proof2:")
print(Proof2)

This should hopefully help stabilize the outcome and automatically apply the sign flip when necessary.

\$\endgroup\$

Not the answer you're looking for? Browse other questions tagged or ask your own question.