In numpy/scipy, what's the canonical way to compute the inverse of an upper triangular matrix?
The matrix is stored as 2D numpy array with zero sub-diagonal elements, and the result should also be stored as a 2D array.
edit The best I've found so far is scipy.linalg.solve_triangular(A, np.identity(n)). Is that it?
There really isn't an inversion routine, per se. scipy.linalg.solve is the canonical way of solving a matrix-vector or matrix-matrix equation, and it can be given explicit information about the structure of the matrix which it will use to choose the correct routine (probably the equivalent of BLAS3 dtrsm in this case).
LAPACK does include doptri for this purpose, and scipy.linalg does expose a raw C lapack interface. If the inverse matrix is really what you want, then you could try using that.
doptri be the right routine? If I understand correctly, the op indicates an orthogonal matrix. Since the poster has a triangular matrix, shouldn't it be tp or tr? I'm not a LAPACK expert--my information comes from: en.wikipedia.org/wiki/LAPACK
dtrtri is perhaps the right one as in scipy.linalg.lapack.clapack.dtrtri
DTRTRI is really the correct answer to this question and should be more visible.I agree that dtrtri should be more visible, so I wrote an example.
# 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))))
At least for this simple example we get:
Condition number: 227524.1404212523
Time inv: 0.1151930999999422
Time dtrtri: 0.03039009999974951
Error inv: 7.883022033401421e-12
Error dtrtri: 7.65486099651801e-13
Which shows that dtrtri() is both faster and accurate than inv(). In this case, both inv() and dtrtri() compute a matrix that is exactly upper triangular. However, this is not the case for a lower triangular matrix, where small entries above the diagonal pollute the result of inv().
numpy.linalg.invis faster thansolve_triangularfor matrices up to about 40x40.