I've long wondered this, but what is the algorithm(s) implemented in modern linear equation solvers for sparse systems?

The obvious answer I can think of is Gauss-Jordan with a bunch of tricks to make it computationally efficient and exploit sparsity, or LU decompositions (e.g. the ones in Eigen), but are there any different algorithms that are commonly used?

  • $\begingroup$ There are some newer matrix decompositions you don't generally see - Non-negative matrix factorizations and the interpolative decomposition are newer. The interpolative decomposition is good for sparse problems - amath.colorado.edu/faculty/martinss/Pubs/…. $\endgroup$ – Ryan Howe Feb 21 at 18:59

The following isn't meant to be exhaustive. It usually depends on the structure of the matrix because that impacts the way you choose it. In general there are sparse variants for many of the general matrix decompositions you see LU, QR, SVD and what not. There are also sparse Kyrlov methods such as the shifted Block Lanczos method. If it isn't symmetric then there is the Arnoldi algorithm which they use in GMRES (there is a sparse variant). There are a bunch of tricks they do use like

  • Preconditioning
  • Restarting
  • Shifting (in terms of eigenvalue methods).
  • The Hessenberg reduction which cuts the computations in half for Krylov methods.
  • Randomization has become very popular in the last 20 years.
  • Convergence speed up methods like Nesterov's method.

Usually the reasoning for choosing one over the other is a combination of the size of the matrix, the conditioning and any relevant structure. MATLAB has a detailed explanation of how the backslash operator works which you can see below. If you look on a website like the NETLIB it has the actual code for most of LAPACK (https://www.netlib.org/).

enter image description here


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.