Linear System Solvers¶

Solving a system of linear equations is a key operation in many scientific applications. More specifically, for a matrix \(\mathbf{A} \in \mathbb{C}^{m \times n}\) and a vector \(\vec{b} \in \mathbb{C}^{m}\), we want to find a solution vector \(\vec{x} \in \mathbb{C}^{n}\) that satisfies

\[\mathbf{A} \vec{x} = \vec{b}\]

If \(m < n\) or \(m > n\), then the system do not have an exact solution, and thus, we are typically looking for the vector \(\vec{x}\) that minimize the norm \(\|\vec{b} - \mathbf{A} \vec{x} \|\) (i.e., the least-square solution). If the matrix \(\mathbf{A}\) is dense and relatively small, then we can use routines from the LAPACK to solve the system with a high degree of accuracy. However, if the matrix is sparse and/or very large, then we need to look for alternative numerical methods.

Sparse Direct Solvers¶

Sparse direct solvers first decompose a sparse matrix \(\mathbf{A}\) into an upper triangular matrix \(\mathbf{U}\) and a lower triangular matrix \(\mathbf{L}\) and then use the backward/forward substitution to find the solution vector \(\vec{x}\). Since most of the elements of \(\mathbf{A}\) are zero, these solvers use supernodal/multifrontal techniques as well as column/row permutation to improve performance and avoid fill-ins. The most computational intensive part is the LU factorization, and thus, it is recommended to compute once the factorization and reuse the factors \(\mathbf{L}\) and \(\mathbf{U}\) to solve multiple linear systems \(\mathbf{A} \vec{x} = \vec{b}\). For example,

using namespace HighFM; CSRMatrix<double> A; Vector<double> x, b; //========== Initialize A, b ==============//

SuperLU<double> solver;

// LU factorization solver.factorize(A);

// Solve 10 systems with the same LU factorization for (int i = 0; i < 10; ++i){

x = solver.solve(b); b = x;

}

AMGCL¶

The iterative solvers are provided by AMGCL. AMGCL is a header-only C++ library for solving large sparse linear systems with algebraic multigrid (AMG) method. AMG is one of the most effective iterative methods for solution of equation systems arising, for example, from discretizing PDEs on unstructured grids. The method can be used as a black-box solver for various computational problems, since it does not require any information about the underlying geometry. AMG is often used not as a standalone solver but as a preconditioner within an iterative solver (e.g. Conjugate Gradients, BiCGStab, or GMRES). See the documentation of AMGCL here.

Usage:

using namespace HighFM;
CSRMatrix<double> A;
Vector<double> x, b;
//========== Initialize A, b ==============//

// Create a preconditioner based on the AMG, Smoothed Aggregation and SPAI0 methods
AMGSA<double, SPAI0> P(A);

// Create a solver using the GMRES algorithm
GMRES<double> gmres(A.rows());

// Solve the system
auto [steps, error] = solve_linsys(P, A, x, b, gmres);

AMGCL

Least-Square¶

LSMR is a variant of the conjugate-gradient method for solving sparse linear equations and sparse least-squares problems:

\[\begin{split}\begin{align*} \text{Solve } & \mathbf{A} \vec{x} = \vec{b} \\ \text{or minimize } & \|\mathbf{A} \vec{x} - \vec{b}\|^2 \\ \text{or minimize } & \|\mathbf{A} \vec{x} - \vec{b}\|^2 + \lambda^2 \|\vec{x}\|^2 \end{align*}\end{split}\]

where the matrix \(\mathbf{A}\) may be square or rectangular (over-determined or under-determined), and may have any rank. The scalar \(\lambda\) is a damping parameter. If \(\lambda > 0\), the solution is “regularized” in the sense that a unique solution always exists, and \(\vec{x}\) is bounded. The method is based on the Golub-Kahan bidiagonalization process. It is algebraically equivalent to applying MINRES to the normal equation \(\mathbf{A}^* \mathbf{A} x = \mathbf{A}^* b\) but has better numerical properties, especially if \(\mathbf{A}\) is ill-conditioned.

LSMR