python - 用 SVD 分解求解线性方程组

Question

我想编写一个函数，它使用 SVD 分解来求解方程组 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

简而言之，使用奇异值分解可以让您将初始问题 A x = b 替换为 U diag(s) Vh 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

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