Fastest way to solve a sparse unsymmetric system many times

I have to solve a system $Ax^{(n)} = b^{(n)}$ many times, $A$ being a sparse (pentadiagonal in most part of its structure), unsymmetric, constant matrix.

Currently, I am performing the LU factorization with rows and columns permutations, i.e., $PAQ = LU$, and then I am using this factorization to compute $x^{(k)}$ in the following way:

$x^{(n)} = Q*(U \backslash L \backslash (P*b^{(n)}))$

I have read recently about iterative methods for sparse systems, but I am not very familiar with them. Could you please tell if there is any iterative method that could perform faster and if it is already implemented in Matlab?

Thanks, Manu

• Since you are solving the system many times with a constant $A$ it is very unlikely you will find an iterative method faster than direct. In matlab the most common iterative method for an unsymmetric system would be gmres but that is often preceded by an approximate LU factorization,ilu, to act as a preconditioner. Trying it is a very easy experiment. – Bill Greene Mar 21 '17 at 18:40

Generally speaking, for many right-hand side (RHS) problems, a direct solver is a more feasible solution for several reasons:

1. Major computations are performed during the factorization step (which is done only once for all RHS), and the solution find (for each RHS) is much cheaper.
2. Direct solvers do not suffer from poor conditioning of the matrix or, in other words, you don't have to look for a nice preconditioner (required for iterative solvers) ensuring $\mathcal O(1)$ iterations for every RHS (relatively easy) and every matrix (much harder) resulting from a numerical simulation.

NB. Sure, if the matrix is extremely ill-conditioned, the direct solver would also suffer, however, direct solvers extend "the dynamic range" of problems that can be solved without re-formulating a problem.

There are a couple of things that I should mention:

1. Usually, it is ill-advised to use an inverse based on LU-decomposition to find a solution. You might want to use a back-substitution LU subroutine that would be more efficient. From your formula, it is unclear if you are doing it in the right way.
2. Since your matrix is sparse and has a stable structure, I would look into sparse direct solvers as opposed to full direct solver. It should speed up the computations significantly.
3. In your case, I would consider doing partial pivoting ($PA=LU$) only, as opposed to full-pivoting ($PAQ=LU$). That might speed up the computations.
4. Connecting a high-performance LAPACK library to MATLAB is also an option.

Regarding sparse iterative solvers:

1. You can try different iterative solvers, including GMRES, BiCGStab, TFQMR, etc. Your task would be to provide a fast sparse matrix-vector product to them. Most of them are already included in MATLAB. And if you are using sparse matrix storage there, using them should not be too hard.
2. Repeating an item from above about advantages of a direct solver: you would have to implement a decent preconditioner for all those methods that would reduce the number of iterations they require to converge. One choice would be using incomplete LU (iLU) as a preconditioner. There are more sophisticated and easier options available.
3. Using iterative solvers is a bit tricky. You have to figure out how well is you system conditioned, what is a reasonable choice for a preconditioner, what is a good tolerance to use with an iterative solver (accuracy of the obtained solution), have a feel what kind of problems (in the subset of problems you are solving using numerical simulations) cause non- or slow-convergence of your iterative algorithm.
• The $PAQ=LU$ factorization does not indicate full rather than partial pivoting. In a sparse solver, the $Q$ permutation is defined to reduce fill-in and the $P$ permutation is used to preserve numerical stability. In MATLAB, the form of sparse LU with both $P$ and $Q$ returned indicates that UMFPACK is being used. Davis provides more details in this paper, pdfs.semanticscholar.org/918c/…, and also mentions that, in some cases, the older MATLAB sparse LU is faster than UMFPACK and how that can be invoked. – Bill Greene Mar 21 '17 at 19:36
• @BillGreene, good to know! Thanks for a nice reference. – Anton Menshov Mar 21 '17 at 19:42
• Thank you for your detailed answer. Currently, I am not inverting L nor U, but using back-substitutions through the backslash command in Matlab (which is already faster than a self-made back-substitution routine, despite it spends some time checking what type of matrices are L and U). I will have a look into direct solvers. – Manu Mar 21 '17 at 20:16
• For iterative solvers, GMRES is optimal (provably converges the fastest), whereas BiCGSTAB, TFQMR are only heuristics but are very cheap. So the sensible thing is to try GMRES first. If it works well, try the others. If not, give up, the others will only perform worse. – Richard Zhang Mar 25 '17 at 15:19
• @RichardZhang, true, only if GMRES restart parameter is equal to $N$ (no restart). Unfortunately, this may require a lot of memory. Thus other solvers worth a shot as well. – Anton Menshov Mar 25 '17 at 16:30