views:

1690

answers:

5

Hello, I'm trying to implement some basic linear algebra operations and one of these operations is the inversion of a triangular (upper and/or lower) matrix. Is there an easy and stable algorithm to do that?

Thank you.

+3  A: 

Yes, use backsubstitution. A standard algorithm to invert a matrix is to find its LU decomposition (decomposition into a lower-triangular and an upper-triangular matrix), use backsubsitution on the triangular pieces, and then combine the results to obtain the inverse of the original matrix.

Jason
I'm trying to get the inverse of a *triangular* matrix, not a square matrix. How back substitution should help me in obtaining the inverses of the triangulars?
tunnuz
By defintion, triangular matrices are square.
Jason
A: 

If you are talking about single precision reals, have a look at the source code for the LAPACK routines STRTRI and STRTI2.

Greg Rogers
A: 

Wow, that's practically half the contents of a numerical analysis course. The standard algorithms will do it, and there is a bunch of canned code here. The ultimate source for this and most other usual numerical analysis problems is Numerical Recipes.

Charlie Martin
Inverting a triangular matrix is not half the contents of a course in numerical analysis. Inverting a triangular matrix is trivial, and the naive algorithm is stable.
Jason
Gotta do row pivoting. Naive won't be stable.
duffymo
Pivoting (via, say, Gaussian elimination) puts a linear system into triangular form which is then solved with backsubstitution. As for stability see, for example, the book Accuracy and Stability of Numerical Algorithms by Nicholas Higham, page 140 of the second edition.
Jason
Uh, Jason, you might want to look up the word "hyperbole". It could come in handy.
Charlie Martin
A: 

Given a lower triangular matrix L, backsubstitution allows you to solve the system L x = b quickly for any right-hand side b.

To invert L, you can solve this system for right-hand sides e1=(1,0,...,0), e2=(0,1,...,0), ..., en=(0,0,...,1) and combine the resulting solution vectors into a single (necessarily lower-triangular) matrix.

If you are interested in a closed-form solution, the diagonal elements of the inverse are the inverses of the original diagonal elements, and the formula for the rest of the elements of the inverse gets more and more complicated as you move aways from the diagonal.

Fanfan
A: 

Don't invert it if you can. It's one of the basic commandments of numerical linear algebra.

It is much faster and numerically stabler to keep the matrix L itself in memory and compute

inv(L)b
with back-substitution whenever you need to do something else with inv(L).

Note that the customary algorithm for inverting it requires solving the systems

inv(L)[1 0 0 ...],
inv(L)[0 1 0 ....],
inv(L)[0 0 1 ....]
and so on, so you see it is much easier not to invert it at all.