Sparse iterative linear solvers

Sparse linear systems are linear systems A·x=b with sparse matrix A. Sparsity of A allows us to use specialized algorithms which may achieve many-orders-of-magnitude speed-up over general purpose dense solvers.

ALGLIB package has well optimized implementations of sparse solvers in all programming languages it supports (C++, C#, Python, ...).

Iterative vs direct methods

All sparse linear solvers belong to one of two families: iterative (also called matrix-free) or direct algorithms.

Iterative solvers access system matrix A only by evaluating its products A·x (or A'·x with trial vectors x provided by the solver. Direct solvers work with explicit representation of A and perform some sparse triangular factorization (Cholesky or LU with pivoting), followed by solution of sparse triangular system.

Both approaches have benefits and drawbacks. From one side, matrix-free methods have modest memory requirements - typically only O(N) additional memory have to be allocated by the algorithm. It contrasts with almost unpredictable memory requirements of direct methods - sparse factorizations often produce so-called fill-in (elements which were zero in the original A, but become non-zero during triangular factorization). And the problem is that in most cases you can not predict amount of fill-in produced (it can be 2x, 10x, 100x increase - no guarantees here).

From the other side, direct methods have nice property of producing more precise answers. And they allow reuse of triangular factorization - say, if you want to solve a sequence of related systems.

Sparse linear CG solver

Sparse linear conjugate gradient algorithm is an iterative algorithm for solution of A·x=b with NxN sparse symmetric positive matrix A. This algorithm does not work for non-positive definite matrices - use LSQR (see below) for such systems. Being purely iterative method, this algorithm has modest - just O(N) - memory requirements (in addition to storing matrix A).

Sparse linear CG solver is provided by lincg subpackage. Typical use case is outlined below:

• first, you should create solver object with lincgcreate function call.
• stopping criteria (step size and/or iterations count) should be set with lincgsetcond call
• initial point can be provided with lincgsetstartingpoint (by default, x=0 is used)
• solution process is started with lincgsolvesparse call, and results can be obtained with lincgresults call

An example (source code) can be found in ALGLIB Reference Manual: lincg_d_1.

Note #1
You may also want to tweak algorithm settings with the help of the following functions: lincgsetprecdiag, lincgsetprecunit, lincgsetrestartfreq, lincgsetrupdatefreq

Sparse LSQR solver

LSQR solver is intended for the solution of sparse rectangular linear systems (in the least squares sense):

as you may see, A can be general rectangular matrix (and no positive definiteness requirement is placed upon it). Additional Tikhonov regularization term of the form λ|x|² can be easily set with this algorithm without explicit modification of A:

The wrong way to approach such tasks is by means of normal equations (A'·A)·x=A'·b. The problem with normal equations is that such reduction squares condition number of the system, with extremely slow convergence. In contrast to normal equations approach, LSQR algorithm does not square condition number (internally it uses variant of Lanczos bidiagonalization process), which results in much faster convergence.

LSQR is also an iterative matrix-free algorithm with O(max(N,M)) additional memory requirements. This solver is provided by linlsqr subpackage. Typical use case is outlined below:

• solver object is created with linlsqrcreate function call.
• stopping criteria (step size and/or iterations count) should be set with linlsqrsetcond call
• additional regularizing term can be set with linlsqrsetlambdai
• solution process is started with linlsqrsolvesparse call, and results can be obtained with linlsqrresults call

An example (source code) can be found in ALGLIB Reference Manual: linlsqr_d_1.

Note #2
You may also want to tweak algorithm settings with the help of the following functions: linlsqrsetprecdiag and linlsqrsetprecunit

This article is licensed for personal use only.