Linear solver recommendation(s) for small problems

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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.