58
$\begingroup$

Sparse linear systems turn up with increasing frequency in applications. One has a lot of routines to choose from for solving these systems. At the highest level, there is a watershed between direct (e.g. sparse Gaussian elimination or Cholesky decomposition, with special ordering algorithms, and multifrontal methods) and iterative (e.g. GMRES, (bi-)conjugate gradient) methods.

How does one determine whether to use a direct or an iterative method? Having made that choice, how does one pick a particular algorithm? I already know about the exploitation of symmetry (e.g. use conjugate gradient for a sparse symmetric positive definite system), but are there any other considerations like this to be considered in picking a method?

Improve this question
$\endgroup$

4 Answers 4

Reset to default
42
$\begingroup$

The important thing when choosing iterative solvers is the spectrum of the operator, see this paper. However, there are so many negative results, see this paper where no iterative solver wins for all problems and this paper in which they prove they can get any convergence curve for GMRES for any spectrum. Thus, it seems impossible to predict the behavior of iterative solvers except in a few isolated cases, Therefore, your best option is to try them all, using a system like PETSc, which also has direct solvers.

Improve this answer
$\endgroup$
3
  • 3
    $\begingroup$ "Throw everything you can at it" was pretty much the advice I was accustomed to. :) The third paper you link to is something I haven't seen before; thanks for that! $\endgroup$
    – J. M.
    Nov 30, 2011 at 15:20
  • 4
    $\begingroup$ Matt has a great answer, but you have to take it within the context of the community he is coming from (large-scale scientific computing). You will find that for small problems (say, less than a hundred thousand unknowns), that direct solvers vastly outperform iterative methods if the problem is not strongly elliptic. I have not seen any good general papers in the literature that would steer you towards an initial starting strategy, which is a bit embarrassing to me. $\endgroup$
    – Aron Ahmadia
    Nov 30, 2011 at 15:42
  • 6
    $\begingroup$ Aron's estimate is good but heavily dependent on fill since sparse direct methods usually exhaust memory before they exhaust patience. $\endgroup$
    – Matt Knepley
    Nov 30, 2011 at 15:50
20
$\begingroup$

The choice between direct and iterative methods is dependent on goals and problem at hand.

For Direct methods, we can note:

  • The coefficient matrix of the linear system changes over the course of computation and may for sparse systems exhaust memory requirements and increase work effort due to fill-in
  • Must complete to give useful results
  • Factorization can be reused in subsequent steps if multiple right-hand sides are present
  • Can be used for solving linear systems only.
  • Seldom fails.

For Iterative methods, we can note:

  • The goal is to give a partial result only after a small number of iterations.
  • Solution effort should be less than direct methods for the same problem.
  • Economical with respect to storage (no fill-in)
  • Often easy to program.
  • A known approximate solution can be exploited.
  • Sometimes they are fast and sometimes they are not (sometimes even divergent).
  • For complex problems, iterative methods are considerably less robust compared to direct methods.

Guidelines for when to use direct or iterative methods?

  • Iterative methods when the coefficient matrix is sparse and direct methods cannot exploit sparsity efficiently (avoid creating fill-in).
  • Direct methods for multiple right-hand sides.
  • Iterative methods can be more efficient if accuracy is of less concern
  • Iterative methods for nonlinear systems of equations.
Improve this answer
$\endgroup$
1
  • 9
    $\begingroup$ I think that it is important to note that direct methods are not always better for multiple right-hand sides. Perhaps they are better for $O(n)$ right-hand sides, but if the iterative method is $O(n)$ while the direct method is $O(n^2)$, it is still advantageous to use the iterative solver for $O(1)$ right-hand sides. $\endgroup$
    – Jack Poulson
    Jan 8, 2012 at 18:08
8
$\begingroup$

I completely concur with the answers already given. I wanted to add that all iterative methods require some sort of initial guess. The quality of this initial guess can often affect the convergence rate of the method you choose. Methods like Jacobi, Gauss Seidel, and Successive Over Relaxation all work to iteratively "smoothen out" as much error as possible at each step (see this paper for details) . The first few steps reduce the high frequency error rather quickly, but the low frequency error takes many more itrations to smoothen out. This is what makes convergence slow for these methods. In cases such as this, we can accelerate convergence by solving low frequency error (e.g. solving the same problem on a coarser mesh) first, then solving the higher frequency error (e.g. on a finer mesh). If we apply this concept recursively by divide and conquer, we get what is called a Multi-grid method. Even if the linear system is not symmetric, there are alternative implementations of multi-grid method for any nonsingular sparse matrix system (e.g. algebraic multi-grid method) which can accelerate the convergence of the solver. Their scalability on parallel systems, however, is the subject of lots of research.

Improve this answer
$\endgroup$
2
  • 6
    $\begingroup$ This answer seems to give the impression that the effectiveness of multigrid comes from finding a good initial guess. In reality, the initial guess is a minor concern for linear problems and really only a concern for Full Multigrid. Multigrid works due to spectral separation. Note that making multigrid perform well for hard problems is a significant challenge. Multigrid works pretty well in parallel, it has been the key ingredient in several Gordon Bell prizes and a few open source packages run with high efficiency on today's largest machines. For GPU implementations, look at the CUSP library. $\endgroup$
    – Jed Brown
    Jan 14, 2012 at 15:42
  • $\begingroup$ Most times a random initial guess is good enough. In extracting eigenvalues using Lanczos algorithm, a random starting/restarting vector does help. Restarts do happen at times in the Lanczos Algorithm. $\endgroup$
    – AnilJ
    Apr 20, 2017 at 23:05
6
$\begingroup$

There is an important piece of information missing in your question: where did the matrix originate from. The structure of the problem you were trying to solve has great potential to suggest a solution method.

If your matrix originated from a partial differential equation with smooth coefficients, a geometric multigrid method will be hard to beat, in particular in three dimensions. If your problem is less regular, algebraic multigrid is a good method. Both usually combined with Krylov-space methods. Other efficient solvers can be derived from fast multipole methods or fast Fourier transform.

Improve this answer
$\endgroup$

Your Answer

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

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