$treeview $search $mathjax
Eigen  3.2.5
$projectbrief
$projectbrief
$searchbox

Solving Sparse Linear Systems
[Sparse linear algebra]

In Eigen, there are several methods available to solve linear systems when the coefficient matrix is sparse. Because of the special representation of this class of matrices, special care should be taken in order to get a good performance. See Sparse matrix manipulations for a detailed introduction about sparse matrices in Eigen. This page lists the sparse solvers available in Eigen. The main steps that are common to all these linear solvers are introduced as well. Depending on the properties of the matrix, the desired accuracy, the end-user is able to tune those steps in order to improve the performance of its code. Note that it is not required to know deeply what's hiding behind these steps: the last section presents a benchmark routine that can be easily used to get an insight on the performance of all the available solvers.

Sparse solvers

Eigen currently provides a limited set of built-in solvers, as well as wrappers to external solver libraries. They are summarized in the following table:

Class

Module

Solver kind

Matrix kind

Features related to performance

Dependencies,License

Notes

SimplicialLLT

SparseCholesky

Direct LLt factorization

SPD

Fill-in reducing

built-in, LGPL

SimplicialLDLT is often preferable

SimplicialLDLT

SparseCholesky

Direct LDLt factorization

SPD

Fill-in reducing

built-in, LGPL

Recommended for very sparse and not too large problems (e.g., 2D Poisson eq.)

ConjugateGradient

IterativeLinearSolvers

Classic iterative CG

SPD

Preconditionning

built-in, MPL2

Recommended for large symmetric problems (e.g., 3D Poisson eq.)

BiCGSTAB

IterativeLinearSolvers

Iterative stabilized bi-conjugate gradient

Square

Preconditionning

built-in, MPL2

To speedup the convergence, try it with the IncompleteLUT preconditioner.

SparseLU

SparseLU

LU factorization

Square

Fill-in reducing, Leverage fast dense algebra

built-in, MPL2

optimized for small and large problems with irregular patterns

SparseQR

SparseQR

QR factorization

Any, rectangular

Fill-in reducing

built-in, MPL2

recommended for least-square problems, has a basic rank-revealing feature

Wrappers to external solvers

PastixLLT
PastixLDLT
PastixLU

PaStiXSupport

Direct LLt, LDLt, LU factorizations

SPD
SPD
Square

Fill-in reducing, Leverage fast dense algebra, Multithreading

Requires the PaStiX package, CeCILL-C

optimized for tough problems and symmetric patterns

CholmodSupernodalLLT

CholmodSupport

Direct LLt factorization

SPD

Fill-in reducing, Leverage fast dense algebra

Requires the SuiteSparse package, GPL

UmfPackLU

UmfPackSupport

Direct LU factorization

Square

Fill-in reducing, Leverage fast dense algebra

Requires the SuiteSparse package, GPL

SuperLU

SuperLUSupport

Direct LU factorization

Square

Fill-in reducing, Leverage fast dense algebra

Requires the SuperLU library, (BSD-like)

SPQR

SPQRSupport

QR factorization

Any, rectangular

fill-in reducing, multithreaded, fast dense algebra

requires the SuiteSparse package, GPL

recommended for linear least-squares problems, has a rank-revealing feature

Here SPD means symmetric positive definite.

All these solvers follow the same general concept. Here is a typical and general example:

#include <Eigen/RequiredModuleName>
// ...
SparseMatrix<double> A;
// fill A
VectorXd b, x;
// fill b
// solve Ax = b
SolverClassName<SparseMatrix<double> > solver;
solver.compute(A);
if(solver.info()!=Success) {
  // decomposition failed
  return;
}
x = solver.solve(b);
if(solver.info()!=Success) {
  // solving failed
  return;
}
// solve for another right hand side:
x1 = solver.solve(b1);

For SPD solvers, a second optional template argument allows to specify which triangular part have to be used, e.g.:

#include <Eigen/IterativeLinearSolvers>

ConjugateGradient<SparseMatrix<double>, Eigen::Upper> solver;
x = solver.compute(A).solve(b);

In the above example, only the upper triangular part of the input matrix A is considered for solving. The opposite triangle might either be empty or contain arbitrary values.

In the case where multiple problems with the same sparsity pattern have to be solved, then the "compute" step can be decomposed as follow:

SolverClassName<SparseMatrix<double> > solver;
solver.analyzePattern(A);   // for this step the numerical values of A are not used
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);
...
A = ...;                    // modify the values of the nonzeros of A, the nonzeros pattern must stay unchanged
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);
...

The compute() method is equivalent to calling both analyzePattern() and factorize().

Finally, each solver provides some specific features, such as determinant, access to the factors, controls of the iterations, and so on. More details are availble in the documentations of the respective classes.

The Compute Step

In the compute() function, the matrix is generally factorized: LLT for self-adjoint matrices, LDLT for general hermitian matrices, LU for non hermitian matrices and QR for rectangular matrices. These are the results of using direct solvers. For this class of solvers precisely, the compute step is further subdivided into analyzePattern() and factorize().

The goal of analyzePattern() is to reorder the nonzero elements of the matrix, such that the factorization step creates less fill-in. This step exploits only the structure of the matrix. Hence, the results of this step can be used for other linear systems where the matrix has the same structure. Note however that sometimes, some external solvers (like SuperLU) require that the values of the matrix are set in this step, for instance to equilibrate the rows and columns of the matrix. In this situation, the results of this step should not be used with other matrices.

Eigen provides a limited set of methods to reorder the matrix in this step, either built-in (COLAMD, AMD) or external (METIS). These methods are set in template parameter list of the solver :

DirectSolverClassName<SparseMatrix<double>, OrderingMethod<IndexType> > solver;

See the OrderingMethods module for the list of available methods and the associated options.

In factorize(), the factors of the coefficient matrix are computed. This step should be called each time the values of the matrix change. However, the structural pattern of the matrix should not change between multiple calls.

For iterative solvers, the compute step is used to eventually setup a preconditioner. For instance, with the ILUT preconditioner, the incomplete factors L and U are computed in this step. Remember that, basically, the goal of the preconditioner is to speedup the convergence of an iterative method by solving a modified linear system where the coefficient matrix has more clustered eigenvalues. For real problems, an iterative solver should always be used with a preconditioner. In Eigen, a preconditioner is selected by simply adding it as a template parameter to the iterative solver object.

IterativeSolverClassName<SparseMatrix<double>, PreconditionerName<SparseMatrix<double> > solver; 

The member function preconditioner() returns a read-write reference to the preconditioner to directly interact with it. See the Iterative solvers module and the documentation of each class for the list of available methods.

The Solve step

The solve() function computes the solution of the linear systems with one or many right hand sides.

X = solver.solve(B);

Here, B can be a vector or a matrix where the columns form the different right hand sides. The solve() function can be called several times as well, for instance when all the right hand sides are not available at once.

x1 = solver.solve(b1);
// Get the second right hand side b2
x2 = solver.solve(b2); 
//  ...

For direct methods, the solution are computed at the machine precision. Sometimes, the solution need not be too accurate. In this case, the iterative methods are more suitable and the desired accuracy can be set before the solve step using setTolerance(). For all the available functions, please, refer to the documentation of the Iterative solvers module .

BenchmarkRoutine

Most of the time, all you need is to know how much time it will take to qolve your system, and hopefully, what is the most suitable solver. In Eigen, we provide a benchmark routine that can be used for this purpose. It is very easy to use. In the build directory, navigate to bench/spbench and compile the routine by typing make spbenchsolver. Run it with --help option to get the list of all available options. Basically, the matrices to test should be in MatrixMarket Coordinate format, and the routine returns the statistics from all available solvers in Eigen.

The following table gives an example of XML statistics from several Eigen built-in and external solvers.

Matrix

N

NNZ

UMFPACK

SUPERLU

PASTIX LU

BiCGSTAB

BiCGSTAB+ILUT

GMRES+ILUT

LDLT

CHOLMOD LDLT

PASTIX LDLT

LLT

CHOLMOD SP LLT

CHOLMOD LLT

PASTIX LLT

CG

vector_graphics

12855

72069

Compute Time

0.0254549

0.0215677

0.0701827

0.000153388

0.0140107

0.0153709

0.0101601

0.00930502

0.0649689

Solve Time

0.00337835

0.000951826

0.00484373

0.0374886

0.0046445

0.00847754

0.000541813

0.000293696

0.00485376

Total Time

0.0288333

0.0225195

0.0750265

0.037642

0.0186552

0.0238484

0.0107019

0.00959871

0.0698227

Error(Iter)

1.299e-16

2.04207e-16

4.83393e-15

3.94856e-11 (80)

1.03861e-12 (3)

5.81088e-14 (6)

1.97578e-16

1.83927e-16

4.24115e-15

poisson_SPD

19788

308232

Compute Time

0.425026

1.82378

0.617367

0.000478921

1.34001

1.33471

0.796419

0.857573

0.473007

0.814826

0.184719

0.861555

0.470559

0.000458188

Solve Time

0.0280053

0.0194402

0.0268747

0.249437

0.0548444

0.0926991

0.00850204

0.0053171

0.0258932

0.00874603

0.00578155

0.00530361

0.0248942

0.239093

Total Time

0.453031

1.84322

0.644241

0.249916

1.39486

1.42741

0.804921

0.862891

0.4989

0.823572

0.190501

0.866859

0.495453

0.239551

Error(Iter)

4.67146e-16

1.068e-15

1.3397e-15

6.29233e-11 (201)

3.68527e-11 (6)

3.3168e-15 (16)

1.86376e-15

1.31518e-16

1.42593e-15

3.45361e-15

3.14575e-16

2.21723e-15

7.21058e-16

9.06435e-12 (261)

sherman2

1080

23094

Compute Time

0.00631754

0.015052

0.0247514

-

0.0214425

0.0217988

Solve Time

0.000478424

0.000337998

0.0010291

-

0.00243152

0.00246152

Total Time

0.00679597

0.01539

0.0257805

-

0.023874

0.0242603

Error(Iter)

1.83099e-15

8.19351e-15

2.625e-14

1.3678e+69 (1080)

4.1911e-12 (7)

5.0299e-13 (12)

bcsstk01_SPD

48

400

Compute Time

0.000169079

0.00010789

0.000572538

1.425e-06

9.1612e-05

8.3985e-05

5.6489e-05

7.0913e-05

0.000468251

5.7389e-05

8.0212e-05

5.8394e-05

0.000463017

1.333e-06

Solve Time

1.2288e-05

1.1124e-05

0.000286387

8.5896e-05

1.6381e-05

1.6984e-05

3.095e-06

4.115e-06

0.000325438

3.504e-06

7.369e-06

3.454e-06

0.000294095

6.0516e-05

Total Time

0.000181367

0.000119014

0.000858925

8.7321e-05

0.000107993

0.000100969

5.9584e-05

7.5028e-05

0.000793689

6.0893e-05

8.7581e-05

6.1848e-05

0.000757112

6.1849e-05

Error(Iter)

1.03474e-16

2.23046e-16

2.01273e-16

4.87455e-07 (48)

1.03553e-16 (2)

3.55965e-16 (2)

2.48189e-16

1.88808e-16

1.97976e-16

2.37248e-16

1.82701e-16

2.71474e-16

2.11322e-16

3.547e-09 (48)

sherman1

1000

3750

Compute Time

0.00228805

0.00209231

0.00528268

9.846e-06

0.00163522

0.00162155

0.000789259

0.000804495

0.00438269

Solve Time

0.000213788

9.7983e-05

0.000938831

0.00629835

0.000361764

0.00078794

4.3989e-05

2.5331e-05

0.000917166

Total Time

0.00250184

0.00219029

0.00622151

0.0063082

0.00199698

0.00240949

0.000833248

0.000829826

0.00529986

Error(Iter)

1.16839e-16

2.25968e-16

2.59116e-16

3.76779e-11 (248)

4.13343e-11 (4)

2.22347e-14 (10)

2.05861e-16

1.83555e-16

1.02917e-15

young1c

841

4089

Compute Time

0.00235843

0.00217228

0.00568075

1.2735e-05

0.00264866

0.00258236

Solve Time

0.000329599

0.000168634

0.00080118

0.0534738

0.00187193

0.00450211

Total Time

0.00268803

0.00234091

0.00648193

0.0534865

0.00452059

0.00708447

Error(Iter)

1.27029e-16

2.81321e-16

5.0492e-15

8.0507e-11 (706)

3.00447e-12 (8)

1.46532e-12 (16)

mhd1280b

1280

22778

Compute Time

0.00234898

0.00207079

0.00570918

2.5976e-05

0.00302563

0.00298036

0.00144525

0.000919922

0.00426444

Solve Time

0.00103392

0.000211911

0.00105

0.0110432

0.000628287

0.00392089

0.000138303

6.2446e-05

0.00097564

Total Time

0.0033829

0.0022827

0.00675918

0.0110692

0.00365392

0.00690124

0.00158355

0.000982368

0.00524008

Error(Iter)

1.32953e-16

3.08646e-16

6.734e-16

8.83132e-11 (40)

1.51153e-16 (1)

6.08556e-16 (8)

1.89264e-16

1.97477e-16

6.68126e-09

crashbasis

160000

1750416

Compute Time

3.2019

5.7892

15.7573

0.00383515

3.1006

3.09921

Solve Time

0.261915

0.106225

0.402141

1.49089

0.24888

0.443673

Total Time

3.46381

5.89542

16.1594

1.49473

3.34948

3.54288

Error(Iter)

1.76348e-16

4.58395e-16

1.67982e-14

8.64144e-11 (61)

8.5996e-12 (2)

6.04042e-14 (5)