Sparse solvers for linear systems

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, ...).

Contents

    1 Iterative vs direct methods
    2 Sparse linear CG solver
    3 Sparse LSQR solver
    4 Downloads section

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:

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| 2 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:

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.

Download ALGLIB for C++ / C# / ...

ALGLIB Project offers you two editions of ALGLIB:

ALGLIB Free Edition:
delivered for free
offers full set of numerical functionality
extensive algorithmic optimizations
no low level optimizations
non-commercial license

ALGLIB Commercial Edition:
flexible pricing
offers full set of numerical functionality
extensive algorithmic optimizations
high performance (SMP, SIMD)
commercial license with support plan

Links to download sections for Free and Commercial editions can be found below:

ALGLIB 3.12.0 for C++

C++ library.
Delivered with sources.
Monolithic design.
Extreme portability.
Editions:   FREE   COMMERCIAL

ALGLIB 3.12.0 for C#

C# library with native kernels.
Delivered with sources.
VB.NET and IronPython wrappers.
Extreme portability.
Editions:   FREE   COMMERCIAL

ALGLIB 3.12.0 for Delphi

Delphi wrapper around C core.
Delivered as precompiled binary.
Compatible with FreePascal.
Editions:   FREE   COMMERCIAL

ALGLIB 3.12.0 for CPython

CPython wrapper around C core.
Delivered as precompiled binary.
Editions:   FREE   COMMERCIAL