Solving Linear Systems of equations with SVD Decomposition

Question

I want to write a function that uses SVD decomposition to solve a system of equations ax=b, where a is a square matrix and b is a vector of values. The scipy function scipy.linalg.svd() should turn a into the matrices U W V. For U and V, I can simply take the transpose of to find their inverse. But for W the function gives me a 1-D array of values that I need to put down the diagonal of a matrix and then enter one over the value.

def solveSVD(a,b):

    U,s,V=sp.svd(a,compute_uv=True)

    Ui=np.transpose(a)
    Vi=np.transpose(V)

    W=np.diag(s)

    Wi=np.empty(np.shape(W)[0],np.shape(W)[1])
    for i in range(np.shape(Wi)[0]):
        if W[i,i]!=0:
            Wi[i,i]=1/W[i,i]
    
    ai=np.matmul(Ui,np.matmul(Wi,Vi))
    x=np.matmul(ai,b)

    return(x)

However, I get a "TypeError: data type not understood" error. I think part of the issue is that

W=np.diag(s)

is not producing a square diagonal matrix.

This is my first time working with this library so apologies if I've done something very stupid, but I cannot work out why this line hasn't worked. Thanks all!

empty function takes a tuple, just use W.shape, and moreover your function does not solve linear system correctly — Yacola
@Yacola Thanks for the help with tuple, can you expand on why my function won't solve LSs correctly? — Eli Rees
Sorry, but I am not sure how were you willing to solve this with your implementation... Is Ui=np.transpose(a) or more like Ui=np.transpose(U), also the for loop could just be replace by Wi=np.diag(1/s)? Check the comments in my code to help you understand what the function is doing. — Yacola

Yacola Yacola · Accepted Answer · 2019-12-11T19:38:05

To be short, using singular value decomposition let you replace your initial problem which is A x = b by U diag(s) Vh x = b. Using a bit of algebra on the latter, give you the following 3 steps function which is really easy to read :

import numpy as np
from scipy.linalg import svd

def solve_svd(A,b):
    # compute svd of A
    U,s,Vh = svd(A)

    # U diag(s) Vh x = b <=> diag(s) Vh x = U.T b = c
    c = np.dot(U.T,b)
    # diag(s) Vh x = c <=> Vh x = diag(1/s) c = w (trivial inversion of a diagonal matrix)
    w = np.dot(np.diag(1/s),c)
    # Vh x = w <=> x = Vh.H w (where .H stands for hermitian = conjugate transpose)
    x = np.dot(Vh.conj().T,w)
    return x

Now, let's test it with

A = np.random.random((100,100))
b = np.random.random((100,1))

and compare it with LU decomposition of np.linalg.solve function

x_svd = solve_svd(A,b)
x_lu = np.linalg.solve(A,b)

which gives

np.allclose(x_lu,x_svd)
>>> True

Please, feel free to ask more explanations in comments if needed. Hope this helps.

Solving Linear Systems of equations with SVD Decomposition

1 Answers