Tutorial page 9 - Sparse Matrix

Table of contents

Manipulating and solving sparse problems involves various modules which are summarized below:

ModuleHeader fileContents
SparseCore
#include <Eigen/SparseCore>
SparseMatrix and SparseVector classes, matrix assembly, basic sparse linear algebra (including sparse triangular solvers)
SparseCholesky Direct sparse LLT and LDLT Cholesky factorization to solve sparse self-adjoint positive definite problems
IterativeLinearSolvers Iterative solvers to solve large general linear square problems (including self-adjoint positive definite problems)
#include <Eigen/Sparse>
Includes all the above modules

Sparse matrix representation

In many applications (e.g., finite element methods) it is common to deal with very large matrices where only a few coefficients are different from zero. In such cases, memory consumption can be reduced and performance increased by using a specialized representation storing only the nonzero coefficients. Such a matrix is called a sparse matrix.

The SparseMatrix class The class SparseMatrix is the main sparse matrix representation of Eigen's sparse module; it offers high performance and low memory usage. It implements a more versatile variant of the widely-used Compressed Column (or Row) Storage scheme. It consists of four compact arrays:

This storage scheme is better explained on an example. The following matrix

03000
2200017
75010
00000
001408

and one of its possible sparse, column major representation:

Values: 227_3514__1_178
InnerIndices: 12_024__2_14
OuterIndexPtrs:03581012
InnerNNZs: 22112

Currently the elements of a given inner vector are guaranteed to be always sorted by increasing inner indices. The "_" indicates available free space to quickly insert new elements. Assuming no reallocation is needed, the insertion of a random element is therefore in O(nnz_j) where nnz_j is the number of nonzeros of the respective inner vector. On the other hand, inserting elements with increasing inner indices in a given inner vector is much more efficient since this only requires to increase the respective InnerNNZs entry that is a O(1) operation.

The case where no empty space is available is a special case, and is refered as the compressed mode. It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS). Any SparseMatrix can be turned to this form by calling the SparseMatrix::makeCompressed() function. In this case, one can remark that the InnerNNZs array is redundant with OuterIndexPtrs because we the equality: InnerNNZs[j] = OuterIndexPtrs[j+1]-OuterIndexPtrs[j]. Therefore, in practice a call to SparseMatrix::makeCompressed() frees this buffer.

It is worth noting that most of our wrappers to external libraries requires compressed matrices as inputs.

The results of Eigen's operations always produces compressed sparse matrices. On the other hand, the insertion of a new element into a SparseMatrix converts this later to the uncompressed mode.

Here is the previous matrix represented in compressed mode:

Values: 22735141178
InnerIndices: 12024214
OuterIndexPtrs:024568

A SparseVector is a special case of a SparseMatrix where only the Values and InnerIndices arrays are stored. There is no notion of compressed/uncompressed mode for a SparseVector.

Matrix and vector properties

Here mat and vec represent any sparse-matrix and sparse-vector type, respectively.

Declarations:

SparseMatrix<std::complex<float> > mat(1000,2000); // declares a 1000x2000 col-major compressed sparse matrix of complex<float>
SparseMatrix<double,RowMajor> mat(1000,2000); // declares a 1000x2000 row-major compressed sparse matrix of double
SparseVector<std::complex<float> > vec(1000); // declares a column sparse vector of complex<float> of size 1000
SparseVector<double,RowMajor> vec(1000); // declares a row sparse vector of double of size 1000
Standard
dimensions
mat.rows()
mat.cols()
vec.size()
Sizes along the
inner/outer dimensions
mat.innerSize()
mat.outerSize()
Number of non
zero coefficients
mat.nonZeros()
vec.nonZeros()

Iterating over the nonzero coefficients

Iterating over the coefficients of a sparse matrix can be done only in the same order as the storage order. Here is an example:

SparseMatrixType mat(rows,cols);
for (int k=0; k<mat.outerSize(); ++k)
for (SparseMatrixType::InnerIterator it(mat,k); it; ++it)
{
it.value();
it.row(); // row index
it.col(); // col index (here it is equal to k)
it.index(); // inner index, here it is equal to it.row()
}
SparseVector<double> vec(size);
for (SparseVector<double>::InnerIterator it(vec); it; ++it)
{
it.value(); // == vec[ it.index() ]
it.index();
}

If the type of the sparse matrix or vector depends on a template parameter, then the typename keyword is required to indicate that InnerIterator denotes a type; see The template and typename keywords in C++ for details.

Filling a sparse matrix

Because of the special storage scheme of a SparseMatrix, special care has to be taken when adding new nonzero entries. For instance, the cost of inserting nnz non zeros in a a single purely random insertion into a SparseMatrix is O(nnz), where nnz is the current number of nonzero coefficients.

The simplest way to create a sparse matrix while guarantying good performance is to first build a list of so called triplets, and then convert it to a SparseMatrix.

Here is a typical usage example:

typedef Triplet<double> T;
std::vector<T> tripletList;
triplets.reserve(estimation_of_entries);
for(...)
{
// ...
tripletList.push_back(T(i,j,v_ij));
}
SparseMatrixType m(rows,cols);
m.setFromTriplets(tripletList.begin(), tripletList.end());
// m is ready to go!

The std::vector triplets might contain the elements in arbitrary order, and might even contain duplicated elements that will be summed up by setFromTriplets(). See the SparseMatrix::setFromTriplets() function and class Triplet for more details.

In some cases, however, slightly higher performance, and lower memory consumption can be reached by directly inserting the non zeros into the destination matrix. A typical scenario of this approach is illustrated bellow:

1: SparseMatrix<double> mat(rows,cols); // default is column major
2: mat.reserve(VectorXi::Constant(cols,6));
3: for each i,j such that v_ij != 0
4: mat.insert(i,j) = v_ij; // alternative: mat.coeffRef(i,j) += v_ij;
5: mat.makeCompressed(); // optional

Supported operators and functions

In the following sm denotes a sparse matrix, sv a sparse vector, dm a dense matrix, and dv a dense vector. In Eigen's sparse module we chose to expose only the subset of the dense matrix API which can be efficiently implemented. Moreover, not every combination is allowed; for instance, it is not possible to add two sparse matrices having two different storage orders. On the other hand, it is perfectly fine to evaluate a sparse matrix or expression to a matrix having a different storage order:

SparseMatrixType sm1, sm2, sm3;
sm3 = sm1.transpose() + sm2; // invalid, because transpose() changes the storage order
sm3 = SparseMatrixType(sm1.transpose()) + sm2; // correct, because evaluation reformats as column-major

Here are some examples of supported operations:

sm1 *= 0.5;
sm1 = sm2 * 0.5;
sm1 = sm2.transpose();
sm1 = sm2.adjoint();
sm4 = sm1 + sm2 + sm3; // only if sm1, sm2 and sm3 have the same storage order
sm3 = sm1 * sm2; // conservative sparse * sparse product preserving numerical zeros
sm3 = (sm1 * sm2).pruned(); // sparse * sparse product that removes numerical zeros (triggers a different algorithm)
sm3 = (sm1 * sm2).pruned(ref); // sparse * sparse product that removes elements much smaller than ref
sm3 = (sm1 * sm2).pruned(ref,epsilon); // sparse * sparse product that removes elements smaller than ref*epsilon
dv3 = sm1 * dv2;
dm3 = sm1 * dm2;
dm3 = dm2 * sm1;
sm3 = sm1.cwiseProduct(sm2); // only if sm1 and sm2 have the same storage order
dv2 = sm1.triangularView<Upper>().solve(dv2);

The product of a sparse symmetric matrix A with a dense matrix (or vector) d can be optimized by specifying the symmetry of A using selfadjointView:

res = A.selfadjointView<>() * d; // if all coefficients of A are stored
res = A.selfadjointView<Upper>() * d; // if only the upper part of A is stored
res = A.selfadjointView<Lower>() * d; // if only the lower part of A is stored

Solving linear problems

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:

ClassModuleSolver kindMatrix kindFeatures related to performance Dependencies,License

Notes

SimplicialLLt SparseCholesky Direct LLt factorizationSPDFill-in reducing built-in, LGPL SimplicialLDLt is often preferable
SimplicialLDLt SparseCholesky Direct LDLt factorizationSPDFill-in reducing built-in, LGPL Recommended for very sparse and not too large problems (e.g., 2D Poisson eq.)
ConjugateGradientIterativeLinearSolvers Classic iterative CGSPDPreconditionning built-in, LGPL Recommended for large symmetric problems (e.g., 3D Poisson eq.)
BiCGSTABIterativeLinearSolvers Iterative stabilized bi-conjugate gradientSquarePreconditionning built-in, LGPL

Might not always converge

CholmodDecompositionCholmodSupport Direct LLT factorizationSPDFill-in reducing, Leverage fast dense algebra Requires the SuiteSparse package, GPL
UmfPackLUUmfPackSupport Direct LU factorizationSquareFill-in reducing, Leverage fast dense algebra Requires the SuiteSparse package, GPL
SuperLUSuperLUSupport Direct LU factorizationSquareFill-in reducing, Leverage fast dense algebra Requires the SuperLU library, (BSD-like)

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()!=Succeeded) {
// decomposition failed
return;
}
x = solver.solve(b);
if(solver.info()!=Succeeded) {
// 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.:

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