# Fastest, and numerically stable way to compute $CA^{-1}B$ and $CA^{-1}x$?

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

I'm not sure if it helps, but $$A$$ is symmetric and positive definite. The $$A$$ matrix will be of dim approximately $$100 \times 100$$ and not sparse...

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

How would we, in Mathematica, compute the expressions above in the fastest way possible, using a numerically stable method? It can be through Cholesky or any other method.

• Do you happen to need that for Schur complement of a symmetric saddlepoint matrix? In that case this survey is a must-read. – Henrik Schumacher Jan 17 '19 at 11:38
• Any extra information about $x$ that we can use? Or totally arbitrary? – MikeY Jan 17 '19 at 17:20
• @MikeY it's arbitrary... – An old man in the sea. Jan 17 '19 at 17:51

In general, you may use the following for computing $$C \, A^{-1} \, x$$:

sol = LinearSolve[a];
c.sol[x]


Here, sol is LinearSolveFunction object which stored a suitable factorization; this factorization can also be reused. Per default, it is an LU-factorization.

In case that a is symmetric positive-definite or Hermitian positive-definite, you can enforce to use a Cholesky factorizaion as follows:

sol = LinearSolve[a, Method -> "Cholesky"];


You can also apply sol to matrices, so $$C \, A^{-1} \, B$$ can be computed with

c.sol[b]


Of course, sol[b] will take much longer than sol[x] because sol has to be applied to each column of b. So if you can, avoid it.

If matrix $$A$$ is not too large but still ill-conditioned, you may also consider to use c.PseudoInverse[a].b. PseudoInverse employs (I guess) singular value decomposition (with truncation of singular values up to some tolerance) and thus can deal also quite well with very ill-conditioned matrices. PseudoInverse is rather expensive and will generate a dense matrix in general. So it is not well-suited for sparse matrices. If you know a priorily that the singular values of a decay rapidly, you can employ a truncated SVD (using only a limited number of the largest singular values); with the Arnoldi method, this might also be doable for certain sparse matrices.

rank = 10;
{U, Σ, V} = SingularValueDecomposition[a, rank, Method -> "Arnoldi"];
c.(U\[Transpose].(1/Diagonal[Σ] V.b))


Unfortunately, there is no way to use perconditioners in Mathematica's Arnoldi methods which really limits their usefulness for large sparse matrices.

• Henrik, in your answer, it seems you didn't end up explaining how one can enforce a Cholesky factorization with linear solve... – An old man in the sea. Jan 17 '19 at 11:57
• Hah! You're right (of course). Please see my edit. – Henrik Schumacher Jan 17 '19 at 12:06
• Henry, thanks for it. However, I was also looking for a faster way to compute those expressions than just using Inverse... I've used AbsoluteTiming, and the linearsolve method is a bit worse... – An old man in the sea. Jan 17 '19 at 12:51
• What kind of matrix is A? E.g., is is exact (i.e., integers) or does it have floating point numbers? – Sjoerd Smit Jan 17 '19 at 14:42
• @SjoerdSmit Floating point numbers. – An old man in the sea. Jan 17 '19 at 17:51