I am interested in solving many linear systems $Ax = b$, where $A$ is symmetric positive definite and small (i.e. less than 25,000 rows) --- $b$ will be changing. We can assume that $A$ arises from the FEM discretization of the elasticity equation in 3D on an unstructured mesh for an arbitrarily-shaped domain.

I am looking for a method that uses less memory than a sparse solver and such that one solution takes less than twice the time for a forward-backward.

Which algorithm would you recommend under these assumptions?

I would be grateful if you could mention why it would be efficient for small matrices (i.e. less than 25,000 rows). For such matrix sizes, I believe that, even if the algorithm is $O(n)$, the constant in front of $n$ matters a lot.

  • 5
    $\begingroup$ At those sizes, nothing beats precomputing the LU factorization of A once and repeatedly applying forward and backward substitutions as b changes. Unless you are working on 10+ years old machines of course. $\endgroup$ Commented Mar 30, 2021 at 15:29
  • 2
    $\begingroup$ because of the symmetry, you could use LDL or cholesky decompositions to save a little memory and time, but yes, you will be hard pressed to beat direct sparse solvers. $\endgroup$ Commented Mar 30, 2021 at 15:34
  • 2
    $\begingroup$ Another important thing to do is reordering the rows/columns of $A$ to reduce fill-in during the Cholesky factorization. $\endgroup$ Commented Mar 30, 2021 at 15:37
  • $\begingroup$ Of course, a good re-ordering would be used. By "a sparse solver", I meant a code that does "LDLT" or Cholesky factorization and that can do forward-backward. So I am looking for something using less memory than a "LDLT" or Cholesky factorization. $\endgroup$
    – user7440
    Commented Mar 30, 2021 at 15:51
  • 1
    $\begingroup$ 25.000 rows is imo everything but small. For example, small is "7" -- just joking. To the point: do you consider a coefficient matrix that is constant, while only the rhs changes? Otherwise (and maybe even then) I can see no difference to the general problem of "solving SPD linear systems" ... $\endgroup$
    – davidhigh
    Commented Mar 31, 2021 at 5:10

3 Answers 3


For a single (maybe a few) $b$'s, I have found Conjugate Gradient can beat an $LL^T$ Cholesky factorization . To do this, I use MKL's Inspector-Executor sparse matrix-vector product mkl_sparse_d_mv with the following:

  1. a good re-ordering (COLAMD and METIS worked for me). This speeds up each matrix-vector product. To avoid a permutation operation each time, I re-ordered the nodes of my FEM mesh such that my $A$ and $b$ is "pre-re-ordered". This may help with your direct methods as well.
  2. prior to solving, I call mkl_sparse_set_mv_hint, which has a HUGE impact on performance. For this to work, you need a good estimate of the required iterations for convergence. The best estimate is to use whatever you needed for the previous solution (simple!). Therefore, the first solve will be a bit slower, but will inform subsequent solves. I suspect this must come at some cost, though I'm not sure what it is.
  3. Jacobi (diagonal) preconditioner. Simple, almost no cost.
  4. Easing up on the residual norm stopping criteria. Instead of requiring machine epsilon, I use something a little bigger, but not so much that I can even notice the difference.
  5. Using a decent initial guess can shave off an iteration or two. If you have access to a really good initial guess, then maybe its significant.

With all these tricks, a supernodal Cholesky factorization (using CHOLMOD) comes pretty close speed-wise. But, it has its uses.


For problems this small, sparse direct solvers are faster than most iterative solvers even if you include the cost of factorization. As a result, I don't believe that you will be able to find a preconditioner for an iterative solver that can work in less than twice the cost of the forward-backward solve. You either have to pay the memory cost of a sparse direct solver, or the runtime cost of an iterative method, but you can't have it both ways.

  • 1
    $\begingroup$ this concisely sums the post that I planned myself (which has the application of computational electromagnetics). $\endgroup$
    – Anton Menshov
    Commented Mar 30, 2021 at 23:38

You can use Krylov iterative solvers with preconditioners. Since you mentioned that your problem of interest is linear elasticity, you end up with symmetric positive definite matrices, assuming that you use standard techniques. For this case, you can use the Conjugate Gradient method with the ILU preconditioner.

These solvers are available in many languages. If you use MATLAB, then you can call in-built functions in MATLAB. For Python, you can use SciPy. In the case of C++, you can use the Eigen library.

  • $\begingroup$ Any reason why it would be cheaper memory-wise? $\endgroup$
    – user7440
    Commented Mar 30, 2021 at 14:17
  • $\begingroup$ The ILU preconditioner allows you to control the amount of storage used in the incomplete LU factorization- in the most extreme case you have 0 fill-in. $\endgroup$ Commented Mar 30, 2021 at 15:34
  • $\begingroup$ My concern is the storage of the directions for the Krylov solver. $\endgroup$
    – user7440
    Commented Mar 30, 2021 at 15:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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