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.
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.
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.
If you are talking about single precision reals, have a look at the source code for the LAPACK routines STRTRI and STRTI2.
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.
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.
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.