SVD Implementation Using QR and Upper Triangular Form

Question

I have to implement svd functino for my term project. I try to use QR method to find eigenvalues and it seems ok but when I try to find eigenvectors using this eigenvalues there is a mistake about the signs of the vectors

Here the part of the Hesenberg form and the functions to find eigenvalues:

def houseeig(A):
    """
    Reduce A to upper Hessenberg form using Householder reflections.
    
    Arguments:
    A -- Input square matrix of size (n x n)
    
    Returns:
    A -- Matrix in upper Hessenberg form
    """
    A = A.astype(np.float64)
    n = A.shape[0]
    for k in range(n-2):
        z = A[k+1:n, k]
        e1 = np.zeros(n-k-1)
        e1[0] = 1
        u = z + np.sign(z[0]) * np.linalg.norm(z) * e1
        u = u / np.linalg.norm(u)
        u = u.reshape(-1,1)
        # print("u",u)
        # Q = np.eye(n-k-1) - 2 * np.outer(u, u.T)
        # print("Q:",Q)
        A[k+1:n, k:n] -= 2* u @ (np.dot(u.T, A[k+1:n, k:n]))
        A[0:n, k+1:n] -= 2 *  (np.dot(A[0:n, k+1:n], u)) @ u.T
        # print("A:",A)
        # print("-------------------------")
    return A

def qrshift(A, tol):
    """
    Find one eigenvalue lam of A in upper Hessenberg form and return iteration count.
    Also improve A for future.

    Arguments:
    A -- Input matrix in upper Hessenberg form
    tol -- Tolerance for convergence

    Returns:
    lam -- Eigenvalue of A
    iter -- Iteration count
    A -- Improved matrix for future
    """
    m = A.shape[0]
    lam = A[m-1, m-1]
    iter = 0
    I = np.eye(m)
    if m == 1:
        return lam, iter, A
    while iter < 100:  # max number of iterations
        if abs(A[m-1, m-2]) < tol:
            return lam, iter, A  # check convergence
        iter += 1
        Q, R = np.linalg.qr(A - lam * I)  # compute the QR decomposition
        A = np.dot(R, Q) + lam * I  # find the next iterate
        lam = A[m-1, m-1]  # next shift
    return lam, iter, A

def qreig(A, tol):
    """
    Find all real eigenvalues lambda of A and return iteration counters.
    First stage: bring to upper Hessenberg form.
    Second stage: deflation loop.
    
    Arguments:
    A -- Input square matrix of size (n x n)
    tol -- Tolerance for convergence
    
    Returns:
    lambda -- List of eigenvalues
    itn -- List of iteration counts
    """
    A = houseeig(A)  # First stage: bring to upper Hessenberg form
    n = A.shape[0]
    eigenvalues = []
    iterations = []
    for j in range(n, 0, -1):
        # find jth eigenvalue
        eigenval, it_count, A = qrshift(A[:j, :j], tol)
        eigenvalues.append(eigenval)
        iterations.append(it_count)
    eigenvalues = eigenvalues[::-1]
    iterations = iterations[::-1]
    return eigenvalues, iterations

And after findigng the eigenvalues I use A - λI = 0 to find eigenvector for specific eigenvalue. First I make A - λI this matrix upper triangular

def upper_triangular(A):
    # Number of rows and columns in the matrix
    A = A.astype(np.float64)
    rows, cols = A.shape
    
    # Iterate over each row
    for i in range(rows):
        # Find the pivot row index (non-zero entry) in the current column
        pivot_index = np.argmax(np.abs(A[i:, i])) + i
        
        # Swap the current row with the pivot row
        A[[i, pivot_index]] = A[[pivot_index, i]]
        
        # Make other elements in the current column zero
        for k in range(i + 1, rows):
            factor = A[k, i] / A[i, i]
            A[k, i:] -= factor * A[i, i:]
    
    return A

Normally there can be some free variable when find the eigenvector but since I cannot give a leter to represent it I assume the x_last is free variable and value equal to 1

def sum_of_xs(A,X,row):
    return np.dot(A[row],X)

def eigenvectors_find(A,eig_vals):
    eigenvectors = np.zeros(A.shape)
    for i in range(len(eig_vals)):
        eigenvectors[:,i] = eigenvector_find(A,eig_vals[i]).flatten()
    return eigenvectors

def eigenvector_find(A,eigen):

    A_ = upper_triangular(A - eigen *np.eye(A.shape[0]).copy())
    

    X = np.zeros((A_.shape[0],1))
    X[-1] = 1
    B = np.zeros(A_.shape[0])

    for i in range(A_.shape[0]-2,-1,-1):
        sum = sum_of_xs(A_,X,i)
        X[i] = -1*sum / A_[i,i]


    X = X / np.linalg.norm(X)
    return X

After all this operations the svd function is

def my_svd(A):
    eigenvals, itr = qreig(A.T@A,1e-12)
    singular_vals = np.sqrt(np.array(eigenvals))
    #find U using A@A.T
    eig_val_u , itr_u = qreig(A@A.T,1e-12)
    U = eigenvectors_find(A@A.T,sorted(eig_val_u,reverse=True))

    #find V using  A.T@A
    eig_vals_v ,itr_v = qreig(A.T@A,1e-12)
    V = eigenvectors_find(A.T@A,sorted(eig_vals_v,reverse=True))

    return U,sorted(singular_vals,reverse=True),V.T

for matrix A

A = np.array([[5,-1,-5,4],
             [-1,3,-2,-3],
             [-3,-2,6,3],
             [1,-3,3,1]])

My function returns

U,sing,V = my_svd(A)
print(U)
print(sing)
print(V)
print("-----------------------------")
print(U @ np.diag(sing) @ V)

[[-0.60936572  0.69094756 -0.27691892  0.27309485]
 [-0.20895515 -0.55114869 -0.12892175  0.79746602]
 [ 0.7110522   0.28235673 -0.57577775  0.28837377]
 [ 0.28180833  0.37296283  0.75840298  0.4542111 ]]
[10.05640977818836, 7.790324714015723, 2.4707266446903646, 0.2738121081517082]
[[-0.46629183 -0.22722116  0.85283687  0.06009774]
 [ 0.45335455 -0.51705105  0.05912314  0.72362443]
 [-0.4978549   0.49924574 -0.18738638  0.68394483]
 [ 0.57374391  0.65710155  0.48379755  0.0705366 ]]
-----------------------------
[[ 5.68125565 -1.68315884 -4.74358387  3.06410227]
 [-0.68283615  2.68195011 -1.88062348 -3.43571444]
 [-1.58351375 -3.42044343  6.53314777  1.05405463]
 [-0.86576745 -1.12902023  2.29774839  3.56316045]]

and numpy svd returns

u,s,v = np.linalg.svd(A)
print(u)
print(s)
print(v)
print("-----------------------------")
print(u @ np.diag(s) @ v)

[[-0.60936572  0.69094756 -0.27691892  0.27309485]
 [-0.20895515 -0.55114869 -0.12892175  0.79746602]
 [ 0.7110522   0.28235673 -0.57577775  0.28837377]
 [ 0.28180833  0.37296283  0.75840298  0.4542111 ]]
[10.05640978  7.79032471  2.47072664  0.27381211]
[[-0.46629183 -0.22722116  0.85283687  0.06009774]
 [ 0.45335455 -0.51705105  0.05912314  0.72362443]
 [ 0.4978549  -0.49924574  0.18738638 -0.68394483]
 [ 0.57374391  0.65710155  0.48379755  0.0705366 ]]
-----------------------------
[[ 5. -1. -5.  4.]
 [-1.  3. -2. -3.]
 [-3. -2.  6.  3.]
 [ 1. -3.  3.  1.]]