16

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?

NPE
  • 486,780
  • 108
  • 951
  • 1,012
  • 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 Jun 14 '13 at 18:19
  • @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 Jul 05 '16 at 14:59

2 Answers2

8

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.

talonmies
  • 70,661
  • 34
  • 192
  • 269
  • 1
    Do you really need the inverted matrix? The LAPACK route is the best one if you really needed it inverted. Otherwise linalg.solve does a pretty decent job using LAPACK to solve a linear system. – Ivan May 18 '11 at 13:43
  • 4
    Would `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: http://en.wikipedia.org/wiki/LAPACK – amcnabb Jun 14 '13 at 18:08
  • 7
    Great pointer but `dtrtri` is perhaps the right one as in `scipy.linalg.lapack.clapack.dtrtri` – dashesy Jul 27 '13 at 23:28
  • 5
    @dashesy `DTRTRI` is really the correct answer to this question and should be more visible. – Daniel Jul 05 '16 at 15:06
3

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().

Jommy
  • 1,020
  • 1
  • 7
  • 14