Matrix inverse using Cholesky decomposition

I have to find a way to calculate the inverse of matrix A using Cholesky decomposition. I understand that using Cholesky we can re-write A^(-1) as A^(-1)=L^(-T) L^(-1) =U^(-1)U^(-T) and the problem is reduced to finding the inverse of the triangular matrix. Obviously I can't use inv(L), because the whole point of this task is to compare several inverse finding methods. However, I fail to figure out how to do this. Is there a built-in function in Matlab for computing an inverse of triangular matrices?

Answer 1

[John D'Errico]

If you have the Cholesky factor, then just use back substitution (twice), applied to an identity matrix. Actually, it will be a forward substitution, then a back subs.

A = L*L'

where L is a lower triangular matrix. So, the problem is to compute the inverse, Ainv here:

L*L'*Ainv = eye(n,n)

Think of it as first solving the problem

L*u = eye(n,n)

L is lower triangular. The right hand side is known, the identity. Apply forward substitution to each column of the RHS.

Then solve a second problem.

L'*Ainv = u

Again, u is now completely known. L' is a upper triangular matrix. Use back substitution to solve for the square matrix Ainv.

Easy, peasy. If you got as far as the Cholesky factor, then this part will be easy for you.

When you are done, VERIFY YOUYR RESULT! Test that

A*Ainv = eye(n,n)

Don't worry if there is some floating point trash in the off diagonal elements of the product if they are on the order of 1e-16. That is to be expected.