# Computing the Inverse of a matrix, using the Cholesky decomposition

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.

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