使用 SVD 分解求解线性方程组

Question

我想编写一个使用 SVD 分解来求解方程组 ax=b 的函数，其中 a 是方阵，b 是值向量。scipy 函数 scipy.linalg.svd() 应该将 a 转换为矩阵 UW V。对于 U 和 V，我可以简单地进行转置来找到它们的逆矩阵。但对于 W，该函数为我提供了一个一维值数组，我需要将其记下矩阵的对角线，然后在该值上输入 1。

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

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