用SVD分解法解线性方程组

Question

我想写一个函数，它使用奇异值分解来求解一个方程组ax=b，其中a是一个方阵，b是一个值向量。scipy函数scipy.linalg.svd()应将a转换为矩阵U、W、V。对于U和V，我只需转置，即可求出它们的逆矩阵。但是对于W，该函数给出了一个一维值的数组，我需要用它来表示矩阵的对角线，然后在该值上输入一个值。

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)

然而，我得到了一个"TypeError:数据类型未被理解“的错误。我认为问题的一部分是

W=np.diag(s) 

并不是生成一个正方形对角矩阵。

这是我第一次使用这个库，所以如果我做了一些非常愚蠢的事情，我很抱歉，但我不明白为什么这行不起作用。谢谢大家！

Answer 1

简而言之，使用奇异值分解可以让你用U diag(s) Vh x = b替换你的初始问题A x = b。在后者上使用一点代数，给出以下3个步骤的函数，非常容易阅读：

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

现在，让我们用以下命令进行测试

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

并与np.linalg.solve函数的LU分解进行了比较

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

这给了我们

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

如果需要，请随时在评论中询问更多解释。希望这能有所帮助。