如何用奇异值分解对numpy矩阵进行反演

Question

(在你告诉我之前，是的，我知道你应该永远不要反转矩阵。不幸的是，在我的计算中，我有一个我构造的矩阵，它必须以某种方式反转。

我有一个大矩阵M，它是病态的。numpy.linalg.cond(M)输出一个大小为e+22的值。基体M的形状为(1000,1000)。

当然，numpy.linalg.inv()将导致许多精度误差。所以，我用numpy.linalg.solve()来求矩阵的逆。

考虑矩阵逆A^{-1}由A * A^{-1} = Identity定义。numpy.linalg.solve()计算已确定的，即满秩线性矩阵方程ax = b的"精确"解x。

那么，我定义单位矩阵:

import numpy as np
iddmatrix = np.identity(100)

和解决:

inverse = np.linalg.solve(M, iddmatrix)

然而，由于我的矩阵如此大且如此病态，np.linalg.solve()不会给出"精确解"。我需要另一种方法来求矩阵的逆。

1. 用SVD实现这种逆的标准方法是什么?
2. 我怎样才能使这个病态矩阵....定义良好?

欢迎提出任何建议。谢谢!

Answer 1

由于SVD将矩阵A分解为U*S*V，其中S是对角的，U和V是正交的，因此它的逆是V'*inv(S)*U'，而对角矩阵的逆就是主对角上的数的逆。

>>> A=np.random.rand(1000,1000)
>>> u,s,v=np.linalg.svd(A)
>>> Ainv=np.dot(v.transpose(),np.dot(np.diag(s**-1),u.transpose()))

Answer 2

考虑取矩阵的SVD实际上意味着什么。这意味着对于某个矩阵M，我们可以将其表示为M=UDV*(这里让*表示转置，因为在堆栈溢出中我看不到这样做的好方法)。

if M=UDV*:
  then: M^-1 = (UDV*)^-1 = (V*^-1)(D^-1)(U^-1)

但是由于U的列是MM*的特征值，V的列是M*M的特征值，所以这些矩阵的逆是它们自己的转置(因为特征向量是正交的)。我们得到:M^-1 = V(D^-1)U*。求对角矩阵的逆就像求这些元素的乘法逆一样简单。

更好的排版:http://adrianboeing.blogspot.com/2010/05/inverting-matrix-svd-singular-value.html

Answer 3

点积的第一个参数应为v.transpose():

import numpy as np
from numpy.linalg import inv
def svdsolve(A):
    u, s, v = np.linalg.svd(A)
    Ainv = np.dot(v.transpose(), np.dot(np.diag(s**-1), u.transpose()))
    return Ainv
temp = np.random.rand(1000, 1000)
np.allclose(svdsolve(temp), inv(temp))
>>> True

np.linalg.solve将初始矩阵分解为A = USV，因此逆矩阵只是V' S^-1 U'