Solving Linear Systems of equations with SVD Decomposition

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!

Solution

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,rcond=1e-8):
    # compute svd of A
    U, s, Vh = svd(A)

    # filtering numerical zeroes
    m = (abs(s) / np.max(abs(s))) > rcond
    U, s, Vh = U[:,m], s[m], Vh[m, :]

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

Now, let's test it with a square matrix

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

to 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

In case you have to deal with a rectangular coefficient matrix, you may compare the resulting array with scipy.linalg.lstsq(A,b) which relies on least-squares solution.

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