numpy：反转上三角矩阵

Question

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

矩阵存储为具有零个下对角线元素的 2D numpy 数组，并且结果也应存储为 2D 数组。

编辑到目前为止我发现的最好的是scipy.linalg.solve_triangle(A, np.identity(n))。是这样吗？

Answer 1

本身确实不存在反转例程。 scipy.linalg.solve 是求解矩阵向量或矩阵矩阵方程的规范方法，并且可以给出有关矩阵结构的明确信息，它将用于选择正确的例程（在本例中可能相当于 BLAS3 dtrsm）。

LAPACK 确实包含用于此目的的 doptri，并且 scipy.linalg 确实公开了原始 C lapack 接口。如果逆矩阵确实是您想要的，那么您可以尝试使用它。

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() 的结果。