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

1 Iterative vs direct methods

2 Sparse linear CG solver

3 Sparse LSQR solver

4 Downloads section

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

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

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.*

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:

C++ library.

Delivered with sources.

Monolithic design.

Extreme portability.

Delivered with sources.

Monolithic design.

Extreme portability.

C# library with native kernels.

Delivered with sources.

VB.NET and IronPython wrappers.

Extreme portability.

Delivered with sources.

VB.NET and IronPython wrappers.

Extreme portability.

Delphi wrapper around C core.

Delivered as precompiled binary.

Compatible with FreePascal.

Delivered as precompiled binary.

Compatible with FreePascal.

CPython wrapper around C core.

Delivered as precompiled binary.

Delivered as precompiled binary.