12
votes

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?

1
How big is the triangular matrix? On my machine, the straightforward numpy.linalg.inv is faster than solve_triangular for matrices up to about 40x40. - amcnabb
@NPE Any update to this? Also have you noted any issues with calling the above (*TRTRS)? My matrix is small enough I can just write a back substitution out for the inverse, but would like to avoid if possible. - Daniel

1 Answers

7
votes

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.