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.

  • $\begingroup$ Do you happen to need that for Schur complement of a symmetric saddlepoint matrix? In that case this survey is a must-read. $\endgroup$ – Henrik Schumacher Jan 17 '19 at 11:38
  • $\begingroup$ Any extra information about $x$ that we can use? Or totally arbitrary? $\endgroup$ – MikeY Jan 17 '19 at 17:20
  • $\begingroup$ @MikeY it's arbitrary... $\endgroup$ – 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]; 

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


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.

| improve this answer | |
  • $\begingroup$ Henrik, in your answer, it seems you didn't end up explaining how one can enforce a Cholesky factorization with linear solve... $\endgroup$ – An old man in the sea. Jan 17 '19 at 11:57
  • $\begingroup$ Hah! You're right (of course). Please see my edit. $\endgroup$ – Henrik Schumacher Jan 17 '19 at 12:06
  • $\begingroup$ 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... $\endgroup$ – An old man in the sea. Jan 17 '19 at 12:51
  • $\begingroup$ What kind of matrix is A? E.g., is is exact (i.e., integers) or does it have floating point numbers? $\endgroup$ – Sjoerd Smit Jan 17 '19 at 14:42
  • $\begingroup$ @SjoerdSmit Floating point numbers. $\endgroup$ – An old man in the sea. Jan 17 '19 at 17:51

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.