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! Continue reading...
