NumPy：求解上三角矩阵的逆矩阵

Question

在 numpy/scipy 中，计算上三角矩阵的逆的规范方法是什么？

该矩阵以2D numpy数组的形式存储，并且具有零子对角线元素，结果也应存储为2D数组。 编辑 我目前找到的最佳方法是scipy.linalg.solve_triangular(A, np.identity(n))。就是这个方法吗？

- NPE

---

三角矩阵有多大？在我的机器上，对于大小不超过40x40的矩阵，直接使用numpy.linalg.inv比solve_triangular更快。 - amcnabb

@NPE 有更新吗？另外，您是否注意到调用上述 (*TRTRS) 存在任何问题？我的矩阵足够小，我可以为逆矩阵编写一个回代算法，但如果可能的话，我想避免这样做。 - Daniel

Answer 1

实际上并没有一个专门的逆转程序，本身。 scipy.linalg.solve是解决矩阵向量或矩阵矩阵方程的规范方式，它可以提供有关矩阵结构的明确信息，用于选择正确的例程（在这种情况下可能相当于BLAS3 dtrsm）。

LAPACK确实包括doptri来实现这个目的，而scipy.linalg确实公开了一个原始的C lapack接口。如果你真的需要逆矩阵，那么你可以尝试使用它。

- talonmies

---

你真的需要求逆矩阵吗？如果你确实需要求逆矩阵，LAPACK是最好的选择。否则，linalg.solve使用LAPACK来解决线性系统，做得相当不错。 - Ivan

“doptri”是否是正确的程序？如果我理解正确，“op”表示正交矩阵。由于发帖者有一个三角矩阵，不应该是“tp”或“tr”吗？我不是LAPACK专家——我的信息来自：http://en.wikipedia.org/wiki/LAPACK - amcnabb

好的，请问需要翻译成什么语言呢？ - dashesy

"DTRTRI" 确实是这个问题的正确答案，应该更加突出显示。 - Daniel

Answer 2

我认为dtrtri应该更加明显，因此我写了一个例子。

# Import.
import timeit
import numpy as np
from scipy.linalg.lapack import dtrtri

# Make a random upper triangular matrix.
rng = np.random.default_rng(12345)
n = 15
mat = np.triu(rng.random(size=(n, n)))
# The condition number is high, and grows quickly with n.
print('Condition number: ', np.linalg.cond(mat))

# Time the generic matrix inverse routine and the ad hoc one.
print('Time inv: ',    timeit.timeit(lambda: np.linalg.inv(mat), number=10000))
print('Time dtrtri: ', timeit.timeit(lambda: dtrtri(mat, lower=0), number=10000))

# Check the error.
inv_mat1 = np.linalg.inv(mat)
inv_mat2, _ = dtrtri(mat, lower=0)
print('Error inv: ',    np.max(np.abs(inv_mat1 @ mat - np.eye(n))))
print('Error dtrtri: ', np.max(np.abs(inv_mat2 @ mat - np.eye(n))))

至少对于这个简单的例子，我们得到：

Condition number:  227524.1404212523
Time inv:  0.1151930999999422
Time dtrtri:  0.03039009999974951
Error inv:  7.883022033401421e-12
Error dtrtri:  7.65486099651801e-13

这表明dtrtri()比inv()更快且更准确。在这种情况下，inv()和dtrtri()都计算了一个完全上三角的矩阵。然而，对于一个下三角矩阵来说，对角线上方的小元素会影响inv()的结果。

- Jommy