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 nonzero coefficients. Such a matrix is called a sparse matrix.
The SparseMatrix class is the main sparse matrix representation of Eigen's sparse module; it offers high performance, low memory usage, and compatibility with most sparse linear algebra packages. These advantages come at the cost of some loss of flexibility, particularly during the assembly of the sparse matrix; consequently, a variant called DynamicSparseMatrix is offered which is tailored for low-level sparse matrix assembly. Both of them can be either row major or column major:
As with dense vectors, the size of a sparse vector denotes its dimension and not the number of nonzero coefficients. At the time of allocation, both sparse matrices and sparse vectors do not have any nonzero coefficients---they correspond to the "all zeros" matrix or vector.
In order to get the most out of Eigen's sparse objects, it is important to have a rough idea of the way they are represented internally. The SparseMatrix class implements the widely-used Compressed Column (or Row) Storage scheme. It consists of three compact arrays: one for the coefficient values, and two for the indices of the nonzero entries. However, the indices are \em not stored as a direct column, row list; instead, the beginning of each column (or row) is encoded as a pointer index. For instance, let \c m be a column-major sparse matrix. Then its nonzero coefficients are sequentially stored in memory in column-major order (\em values). A second array of integers stores the respective row index of each coefficient (\em inner \em indices). Finally, a third array of integers, having the same length as the number of columns, stores the index in the previous arrays of the first element of each column (\em outer \em indices).
As you might guess, here the storage order is even more important than with dense matrices. We will therefore often make a clear difference between the \em inner and \em outer dimensions. For instance, it is efficient to loop over the coefficients of an \em inner \em vector (e.g., a column of a column-major matrix), but completely inefficient to do the same for an \em outer \em vector (e.g., a row of a column-major matrix).
Since all nonzero coefficients of such a matrix are sequentially stored in memory, inserting a new nonzero near the "beginning" of the matrix can be extremely costly. As described below (\ref TutorialSparseFilling), one strategy is to fill nonzero coefficients in order. In cases where this is not possible, Eigen's sparse module also provides a DynamicSparseMatrix class which allows efficient random insertion. DynamicSparseMatrix is essentially implemented as an array of SparseVector, where the values and inner-indices arrays have been split into multiple small and resizable arrays. Assuming the number of nonzeros per inner vector is relatively small, this modification allows for very fast random insertion at the cost of a slight memory overhead (due to extra memory preallocated by each inner vector to avoid an expensive memory reallocation at every insertion) and a loss of compatibility with other sparse libraries used by some of our high-level solvers. Once complete, a DynamicSparseMatrix can be converted to a SparseMatrix to permit usage of these sparse libraries.
To summarize, it is recommended to use SparseMatrix whenever possible, and reserve the use of DynamicSparseMatrix to assemble a sparse matrix in cases when a SparseMatrix is not flexible enough. The respective pros/cons of both representations are summarized in the following table:
Because of the special storage scheme of a SparseMatrix, adding new nonzero entries can have consequences for performance. For instance, the cost of a purely random insertion into a SparseMatrix is O(nnz), where nnz is the current number of nonzero coefficients. In order to cover all use cases with best efficiency, Eigen provides various mechanisms, from the easiest but slowest, to the fastest but most restrictive.
If you don't have any prior knowledge about the order your matrix will be filled, then the best choice is to use a DynamicSparseMatrix. With a DynamicSparseMatrix, you can add or modify any coefficients at any time using the coeffRef(row,col) method. Here is an example:
- \code m.insert(i,j) = value; \endcode which assumes the coefficient of coordinate (i,j) does not already exist (otherwise this is a programming error and your program will stop).
- \code m.insertBack(i,j) = value; \endcode which, in addition to the requirements of insert(), also assumes that the coefficient of coordinate (i,j) will be inserted at the end of the target inner-vector. More precisely, if the matrix m is column major, then the row index of the last non zero coefficient of the j-th column must be smaller than i.
The SparseMatrix class also supports random insertion via the insert() method. However, it should only be used when the inserted coefficient is nearly the last one of the compact storage array. In practice, this means it should be used only to perform random (or sorted) insertion into the current inner-vector while filling the inner-vectors in increasing order. Moreover, with a SparseMatrix an insertion session must be closed by a call to finalize() before any use of the matrix. Here is an example for a column major matrix:
Finally, the fastest way to fill a SparseMatrix object is to insert the elements in purely increasing order (increasing inner index per outer index, and increasing outer index) using the insertBack() function:
In the following \em sm denotes a sparse matrix, \em sv a sparse vector, \em dm a dense matrix, and \em 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:
The product of a sparse \em symmetric matrix A with a dense matrix (or vector) d can be optimized by specifying the symmetry of A using selfadjointView:
\warning Those modules are currently considered to be unsupported because 1) they are not documented, and 2) their API is likely to change in the future.
