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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.