为啥 numpy.linalg.solve() 提供比 numpy.linalg.inv() 更精确的矩阵求逆？

Posted 2023-02-23

【中文标题】为啥 numpy.linalg.solve() 提供比 numpy.linalg.inv() 更精确的矩阵求逆？

【英文标题】：Why does numpy.linalg.solve() offer more precise matrix inversions than numpy.linalg.inv()?为什么 numpy.linalg.solve() 提供比 numpy.linalg.inv() 更精确的矩阵求逆？

【发布时间】：2015-09-24 05:08:13

【问题描述】：

我不太明白为什么numpy.linalg.solve() 给出了更准确的答案，而numpy.linalg.inv() 则有些崩溃，给出了（我相信的）估计值。

举一个具体的例子，我正在求解方程C^-1 * d，其中C 表示一个矩阵，d 是一个向量数组。为便于讨论，C 的尺寸为形状(1000,1000)，d 的尺寸为形状(1,1000)。

numpy.linalg.solve(A, b) 为 x 求解方程 A*x=b，即 x = A^-1 * b. 因此，我可以通过以下方式求解这个方程

(1)

inverse = numpy.linalg.inv(C)
result = inverse * d

或 (2)

numpy.linalg.solve(C, d)

方法 (2) 给出了更精确的结果。为什么是这样？

究竟发生了什么使得一个“工作得更好”比另一个？

【问题讨论】：

最好拿起一本教科书。 Golub 和 van Loan 将是一个不错的选择。

只是一个猜测（可能来自 Gilbert Strang Books 的某处）：我认为在调用底层 LAPACK 例程期间从未真正计算过逆。如果您尝试 pinv() 而不是 inv() 怎么办？看pinv()的引用

在大多数情况下，线性方程解法比反演和乘法更快，产生的舍入误差更小。从根本上说，这些是对“浮动”数字表示的操作，而不是实数本身。

知道pinv与solve的比较会很有趣

你几乎从不想计算一个矩阵的逆......它效率低下，而且误差通常高于其他方法。

【参考方案1】：

np.linalg.solve(A, b)不计算A的逆。相反，它调用gesv LAPACK routines 之一，它首先使用LU 分解分解A，然后使用前向和后向替换求解x（参见here）。

np.linalg.inv 使用相同的方法来计算 A 的逆，通过在 A·A 中求解 A^{-1 -1} = I 其中 I 是身份*。分解步骤与上面完全相同，但需要更多的浮点运算来求解A^-1（一个n×n矩阵) 而不是 x (一个 n-long 向量)。此外，如果您想通过身份 A^-1·b = x 获得 x，那么额外的矩阵乘法将导致更多的浮动点操作，因此性能更慢，数值错误更多。

不需要计算A^-1的中间步骤——直接得到x更快更准确。 p>

---

* inv 的相关源位是here。不幸的是，由于它是模板化的 C，因此有点难以理解。需要注意的重要一点是，单位矩阵作为参数 B 传递给 LAPACK 求解器。