Computing the Inverse of a matrix, using the Cholesky decomposition

Question

I have to compute $CA^{-1}B$ and $CA^{-1}x$, where $A,B,C$ are conformable matrices and $x$ is a vector.

I've read that the a very computationally stable way to compute these inverses is by computing the Cholesky Decomposition.

However, I don't see how one can use that decomposition to compute Both terms I'm interested in...

Any help would be appreciated.

Answer 1

The Cholesky decomposition [the function dpotrf() in LAPACK] factors $\mathbf A = \mathbf L \mathbf L^{\mathrm T}$, or alternatively $\mathbf A^{-1} = \left(\mathbf L \mathbf L^\mathrm T \right)^{-1} = \mathbf L^{-\mathrm T}\mathbf L^{-1}$. If you insert the latter representation your other expressions you'll see how you can compute them efficiently:

• $\mathbf C \mathbf A^{-1} \mathbf B = \mathbf C \mathbf L^{-T}\mathbf L^{-1} \mathbf B$: if you allow yourself to overwrite $\mathbf B$, the quantity $\mathbf{\tilde B} = \mathbf L^{-T}\mathbf L^{-1} \mathbf B$ can be computed efficiently using two consecutive triangle\matrix backsolutions [the function dtrsm() in BLAS]. Then the sought quantity $\mathbf C \mathbf{\tilde B}$ can be computed using matrix-matrix multiplication [the function dgemm() in BLAS].

• $\mathbf C \mathbf A^{-1} \mathbf x = \mathbf C \mathbf L^{-T}\mathbf L^{-1} \mathbf x$: similar approach, overwrite $\mathbf{\tilde x} = \mathbf L^{-T}\mathbf L^{-1} \mathbf x$ using two triangle\vector backsolutions [dtrsv() in BLAS]. Then $\mathbf C \mathbf {\tilde x}$ follows from matrix-vector multiplication [dgemv() in BLAS].

• It is worth mentioning that a symmetric product involving $\mathbf A^{-1}$ leads to an especially efficient algorithm. For example, the quantity $\mathbf B^{\mathrm T} \mathbf A^{-1} \mathbf B = \mathbf B^{\mathrm T} \mathbf L^{-T}\mathbf L^{-1} \mathbf B $ can be found by introducing $\mathbf{\tilde B} = \mathbf L^{-1} \mathbf B$ [requires only one dstrm() call] and then using symmetric-rank-k multication to form $\mathbf {\tilde B}^{\mathrm T} \mathbf {\tilde B}$ [the function dsyrk() in BLAS]. This is a factor of 2 savings.

• A similar idea can be applied to computing the quadratic form $\mathbf x^{\mathrm T}\mathbf A^{-1} \mathbf x$ [one dtrsv() call, followed by ddot()]