Home | History | Annotate | Download | only in doc
      1 namespace Eigen {
      2 
      3 /** \eigenManualPage TutorialSparse Sparse matrix manipulations
      4 
      5 \eigenAutoToc
      6 
      7 Manipulating and solving sparse problems involves various modules which are summarized below:
      8 
      9 <table class="manual">
     10 <tr><th>Module</th><th>Header file</th><th>Contents</th></tr>
     11 <tr><td>\link SparseCore_Module SparseCore \endlink</td><td>\code#include <Eigen/SparseCore>\endcode</td><td>SparseMatrix and SparseVector classes, matrix assembly, basic sparse linear algebra (including sparse triangular solvers)</td></tr>
     12 <tr><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>\code#include <Eigen/SparseCholesky>\endcode</td><td>Direct sparse LLT and LDLT Cholesky factorization to solve sparse self-adjoint positive definite problems</td></tr>
     13 <tr><td>\link SparseLU_Module SparseLU \endlink</td><td>\code #include<Eigen/SparseLU> \endcode</td>
     14 <td>%Sparse LU factorization to solve general square sparse systems</td></tr>
     15 <tr><td>\link SparseQR_Module SparseQR \endlink</td><td>\code #include<Eigen/SparseQR>\endcode </td><td>%Sparse QR factorization for solving sparse linear least-squares problems</td></tr>
     16 <tr><td>\link IterativeLinearSolvers_Module IterativeLinearSolvers \endlink</td><td>\code#include <Eigen/IterativeLinearSolvers>\endcode</td><td>Iterative solvers to solve large general linear square problems (including self-adjoint positive definite problems)</td></tr>
     17 <tr><td>\link Sparse_Module Sparse \endlink</td><td>\code#include <Eigen/Sparse>\endcode</td><td>Includes all the above modules</td></tr>
     18 </table>
     19 
     20 \section TutorialSparseIntro Sparse matrix format
     21 
     22 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.
     23 
     24 \b The \b %SparseMatrix \b class
     25 
     26 The class SparseMatrix is the main sparse matrix representation of Eigen's sparse module; it offers high performance and low memory usage.
     27 It implements a more versatile variant of the widely-used Compressed Column (or Row) Storage scheme.
     28 It consists of four compact arrays:
     29  - \c Values: stores the coefficient values of the non-zeros.
     30  - \c InnerIndices: stores the row (resp. column) indices of the non-zeros.
     31  - \c OuterStarts: stores for each column (resp. row) the index of the first non-zero in the previous two arrays.
     32  - \c InnerNNZs: stores the number of non-zeros of each column (resp. row).
     33 The word \c inner refers to an \em inner \em vector that is a column for a column-major matrix, or a row for a row-major matrix.
     34 The word \c outer refers to the other direction.
     35 
     36 This storage scheme is better explained on an example. The following matrix
     37 <table class="manual">
     38 <tr><td> 0</td><td>3</td><td> 0</td><td>0</td><td> 0</td></tr>
     39 <tr><td>22</td><td>0</td><td> 0</td><td>0</td><td>17</td></tr>
     40 <tr><td> 7</td><td>5</td><td> 0</td><td>1</td><td> 0</td></tr>
     41 <tr><td> 0</td><td>0</td><td> 0</td><td>0</td><td> 0</td></tr>
     42 <tr><td> 0</td><td>0</td><td>14</td><td>0</td><td> 8</td></tr>
     43 </table>
     44 
     45 and one of its possible sparse, \b column \b major representation:
     46 <table class="manual">
     47 <tr><td>Values:</td>        <td>22</td><td>7</td><td>_</td><td>3</td><td>5</td><td>14</td><td>_</td><td>_</td><td>1</td><td>_</td><td>17</td><td>8</td></tr>
     48 <tr><td>InnerIndices:</td>  <td> 1</td><td>2</td><td>_</td><td>0</td><td>2</td><td> 4</td><td>_</td><td>_</td><td>2</td><td>_</td><td> 1</td><td>4</td></tr>
     49 </table>
     50 <table class="manual">
     51 <tr><td>OuterStarts:</td><td>0</td><td>3</td><td>5</td><td>8</td><td>10</td><td>\em 12 </td></tr>
     52 <tr><td>InnerNNZs:</td>    <td>2</td><td>2</td><td>1</td><td>1</td><td> 2</td><td></td></tr>
     53 </table>
     54 
     55 Currently the elements of a given inner vector are guaranteed to be always sorted by increasing inner indices.
     56 The \c "_" indicates available free space to quickly insert new elements.
     57 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.
     58 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 \c InnerNNZs entry that is a O(1) operation.
     59 
     60 The case where no empty space is available is a special case, and is refered as the \em compressed mode.
     61 It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS).
     62 Any SparseMatrix can be turned to this form by calling the SparseMatrix::makeCompressed() function.
     63 In this case, one can remark that the \c InnerNNZs array is redundant with \c OuterStarts because we the equality: \c InnerNNZs[j] = \c OuterStarts[j+1]-\c OuterStarts[j].
     64 Therefore, in practice a call to SparseMatrix::makeCompressed() frees this buffer.
     65 
     66 It is worth noting that most of our wrappers to external libraries requires compressed matrices as inputs.
     67 
     68 The results of %Eigen's operations always produces \b compressed sparse matrices.
     69 On the other hand, the insertion of a new element into a SparseMatrix converts this later to the \b uncompressed mode.
     70 
     71 Here is the previous matrix represented in compressed mode:
     72 <table class="manual">
     73 <tr><td>Values:</td>        <td>22</td><td>7</td><td>3</td><td>5</td><td>14</td><td>1</td><td>17</td><td>8</td></tr>
     74 <tr><td>InnerIndices:</td>  <td> 1</td><td>2</td><td>0</td><td>2</td><td> 4</td><td>2</td><td> 1</td><td>4</td></tr>
     75 </table>
     76 <table class="manual">
     77 <tr><td>OuterStarts:</td><td>0</td><td>2</td><td>4</td><td>5</td><td>6</td><td>\em 8 </td></tr>
     78 </table>
     79 
     80 A SparseVector is a special case of a SparseMatrix where only the \c Values and \c InnerIndices arrays are stored.
     81 There is no notion of compressed/uncompressed mode for a SparseVector.
     82 
     83 
     84 \section TutorialSparseExample First example
     85 
     86 Before describing each individual class, let's start with the following typical example: solving the Laplace equation \f$ \Delta u = 0 \f$ on a regular 2D grid using a finite difference scheme and Dirichlet boundary conditions.
     87 Such problem can be mathematically expressed as a linear problem of the form \f$ Ax=b \f$ where \f$ x \f$ is the vector of \c m unknowns (in our case, the values of the pixels), \f$ b \f$ is the right hand side vector resulting from the boundary conditions, and \f$ A \f$ is an \f$ m \times m \f$ matrix containing only a few non-zero elements resulting from the discretization of the Laplacian operator.
     88 
     89 <table class="manual">
     90 <tr><td>
     91 \include Tutorial_sparse_example.cpp
     92 </td>
     93 <td>
     94 \image html Tutorial_sparse_example.jpeg
     95 </td></tr></table>
     96 
     97 In this example, we start by defining a column-major sparse matrix type of double \c SparseMatrix<double>, and a triplet list of the same scalar type \c  Triplet<double>. A triplet is a simple object representing a non-zero entry as the triplet: \c row index, \c column index, \c value.
     98 
     99 In the main function, we declare a list \c coefficients of triplets (as a std vector) and the right hand side vector \f$ b \f$ which are filled by the \a buildProblem function.
    100 The raw and flat list of non-zero entries is then converted to a true SparseMatrix object \c A.
    101 Note that the elements of the list do not have to be sorted, and possible duplicate entries will be summed up.
    102 
    103 The last step consists of effectively solving the assembled problem.
    104 Since the resulting matrix \c A is symmetric by construction, we can perform a direct Cholesky factorization via the SimplicialLDLT class which behaves like its LDLT counterpart for dense objects.
    105 
    106 The resulting vector \c x contains the pixel values as a 1D array which is saved to a jpeg file shown on the right of the code above.
    107 
    108 Describing the \a buildProblem and \a save functions is out of the scope of this tutorial. They are given \ref TutorialSparse_example_details "here" for the curious and reproducibility purpose.
    109 
    110 
    111 
    112 
    113 \section TutorialSparseSparseMatrix The SparseMatrix class
    114 
    115 \b %Matrix \b and \b vector \b properties \n
    116 
    117 The SparseMatrix and SparseVector classes take three template arguments:
    118  * the scalar type (e.g., double)
    119  * the storage order (ColMajor or RowMajor, the default is ColMajor)
    120  * the inner index type (default is \c int).
    121 
    122 As for dense Matrix objects, constructors takes the size of the object.
    123 Here are some examples:
    124 
    125 \code
    126 SparseMatrix<std::complex<float> > mat(1000,2000);         // declares a 1000x2000 column-major compressed sparse matrix of complex<float>
    127 SparseMatrix<double,RowMajor> mat(1000,2000);              // declares a 1000x2000 row-major compressed sparse matrix of double
    128 SparseVector<std::complex<float> > vec(1000);              // declares a column sparse vector of complex<float> of size 1000
    129 SparseVector<double,RowMajor> vec(1000);                   // declares a row sparse vector of double of size 1000
    130 \endcode
    131 
    132 In the rest of the tutorial, \c mat and \c vec represent any sparse-matrix and sparse-vector objects, respectively.
    133 
    134 The dimensions of a matrix can be queried using the following functions:
    135 <table class="manual">
    136 <tr><td>Standard \n dimensions</td><td>\code
    137 mat.rows()
    138 mat.cols()\endcode</td>
    139 <td>\code
    140 vec.size() \endcode</td>
    141 </tr>
    142 <tr><td>Sizes along the \n inner/outer dimensions</td><td>\code
    143 mat.innerSize()
    144 mat.outerSize()\endcode</td>
    145 <td></td>
    146 </tr>
    147 <tr><td>Number of non \n zero coefficients</td><td>\code
    148 mat.nonZeros() \endcode</td>
    149 <td>\code
    150 vec.nonZeros() \endcode</td></tr>
    151 </table>
    152 
    153 
    154 \b Iterating \b over \b the \b nonzero \b coefficients \n
    155 
    156 Random access to the elements of a sparse object can be done through the \c coeffRef(i,j) function.
    157 However, this function involves a quite expensive binary search.
    158 In most cases, one only wants to iterate over the non-zeros elements. This is achieved by a standard loop over the outer dimension, and then by iterating over the non-zeros of the current inner vector via an InnerIterator. Thus, the non-zero entries have to be visited in the same order than the storage order.
    159 Here is an example:
    160 <table class="manual">
    161 <tr><td>
    162 \code
    163 SparseMatrix<double> mat(rows,cols);
    164 for (int k=0; k<mat.outerSize(); ++k)
    165   for (SparseMatrix<double>::InnerIterator it(mat,k); it; ++it)
    166   {
    167     it.value();
    168     it.row();   // row index
    169     it.col();   // col index (here it is equal to k)
    170     it.index(); // inner index, here it is equal to it.row()
    171   }
    172 \endcode
    173 </td><td>
    174 \code
    175 SparseVector<double> vec(size);
    176 for (SparseVector<double>::InnerIterator it(vec); it; ++it)
    177 {
    178   it.value(); // == vec[ it.index() ]
    179   it.index();
    180 }
    181 \endcode
    182 </td></tr>
    183 </table>
    184 For a writable expression, the referenced value can be modified using the valueRef() function.
    185 If the type of the sparse matrix or vector depends on a template parameter, then the \c typename keyword is
    186 required to indicate that \c InnerIterator denotes a type; see \ref TopicTemplateKeyword for details.
    187 
    188 
    189 \section TutorialSparseFilling Filling a sparse matrix
    190 
    191 Because of the special storage scheme of a SparseMatrix, special care has to be taken when adding new nonzero entries.
    192 For instance, the cost of a single purely random insertion into a SparseMatrix is \c O(nnz), where \c nnz is the current number of non-zero coefficients.
    193 
    194 The simplest way to create a sparse matrix while guaranteeing good performance is thus to first build a list of so-called \em triplets, and then convert it to a SparseMatrix.
    195 
    196 Here is a typical usage example:
    197 \code
    198 typedef Eigen::Triplet<double> T;
    199 std::vector<T> tripletList;
    200 tripletList.reserve(estimation_of_entries);
    201 for(...)
    202 {
    203   // ...
    204   tripletList.push_back(T(i,j,v_ij));
    205 }
    206 SparseMatrixType mat(rows,cols);
    207 mat.setFromTriplets(tripletList.begin(), tripletList.end());
    208 // mat is ready to go!
    209 \endcode
    210 The \c std::vector of triplets might contain the elements in arbitrary order, and might even contain duplicated elements that will be summed up by setFromTriplets().
    211 See the SparseMatrix::setFromTriplets() function and class Triplet for more details.
    212 
    213 
    214 In some cases, however, slightly higher performance, and lower memory consumption can be reached by directly inserting the non-zeros into the destination matrix.
    215 A typical scenario of this approach is illustrated bellow:
    216 \code
    217 1: SparseMatrix<double> mat(rows,cols);         // default is column major
    218 2: mat.reserve(VectorXi::Constant(cols,6));
    219 3: for each i,j such that v_ij != 0
    220 4:   mat.insert(i,j) = v_ij;                    // alternative: mat.coeffRef(i,j) += v_ij;
    221 5: mat.makeCompressed();                        // optional
    222 \endcode
    223 
    224 - The key ingredient here is the line 2 where we reserve room for 6 non-zeros per column. In many cases, the number of non-zeros per column or row can easily be known in advance. If it varies significantly for each inner vector, then it is possible to specify a reserve size for each inner vector by providing a vector object with an operator[](int j) returning the reserve size of the \c j-th inner vector (e.g., via a VectorXi or std::vector<int>). If only a rought estimate of the number of nonzeros per inner-vector can be obtained, it is highly recommended to overestimate it rather than the opposite. If this line is omitted, then the first insertion of a new element will reserve room for 2 elements per inner vector.
    225 - The line 4 performs a sorted insertion. In this example, the ideal case is when the \c j-th column is not full and contains non-zeros whose inner-indices are smaller than \c i. In this case, this operation boils down to trivial O(1) operation.
    226 - When calling insert(i,j) the element \c i \c ,j must not already exists, otherwise use the coeffRef(i,j) method that will allow to, e.g., accumulate values. This method first performs a binary search and finally calls insert(i,j) if the element does not already exist. It is more flexible than insert() but also more costly.
    227 - The line 5 suppresses the remaining empty space and transforms the matrix into a compressed column storage.
    228 
    229 
    230 
    231 \section TutorialSparseFeatureSet Supported operators and functions
    232 
    233 Because of their special storage format, sparse matrices cannot offer the same level of flexibility than dense matrices.
    234 In Eigen's sparse module we chose to expose only the subset of the dense matrix API which can be efficiently implemented.
    235 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.
    236 
    237 \subsection TutorialSparse_BasicOps Basic operations
    238 
    239 %Sparse expressions support most of the unary and binary coefficient wise operations:
    240 \code
    241 sm1.real()   sm1.imag()   -sm1                    0.5*sm1
    242 sm1+sm2      sm1-sm2      sm1.cwiseProduct(sm2)
    243 \endcode
    244 However, <strong>a strong restriction is that the storage orders must match</strong>. For instance, in the following example:
    245 \code
    246 sm4 = sm1 + sm2 + sm3;
    247 \endcode
    248 sm1, sm2, and sm3 must all be row-major or all column-major.
    249 On the other hand, there is no restriction on the target matrix sm4.
    250 For instance, this means that for computing \f$ A^T + A \f$, the matrix \f$ A^T \f$ must be evaluated into a temporary matrix of compatible storage order:
    251 \code
    252 SparseMatrix<double> A, B;
    253 B = SparseMatrix<double>(A.transpose()) + A;
    254 \endcode
    255 
    256 Binary coefficient wise operators can also mix sparse and dense expressions:
    257 \code
    258 sm2 = sm1.cwiseProduct(dm1);
    259 dm2 = sm1 + dm1;
    260 dm2 = dm1 - sm1;
    261 \endcode
    262 Performance-wise, the adding/subtracting sparse and dense matrices is better performed in two steps. For instance, instead of doing <tt>dm2 = sm1 + dm1</tt>, better write:
    263 \code
    264 dm2 = dm1;
    265 dm2 += sm1;
    266 \endcode
    267 This version has the advantage to fully exploit the higher performance of dense storage (no indirection, SIMD, etc.), and to pay the cost of slow sparse evaluation on the few non-zeros of the sparse matrix only.
    268 
    269 
    270 %Sparse expressions also support transposition:
    271 \code
    272 sm1 = sm2.transpose();
    273 sm1 = sm2.adjoint();
    274 \endcode
    275 However, there is no transposeInPlace() method.
    276 
    277 
    278 \subsection TutorialSparse_Products Matrix products
    279 
    280 %Eigen supports various kind of sparse matrix products which are summarize below:
    281   - \b sparse-dense:
    282     \code
    283 dv2 = sm1 * dv1;
    284 dm2 = dm1 * sm1.adjoint();
    285 dm2 = 2. * sm1 * dm1;
    286     \endcode
    287   - \b symmetric \b sparse-dense. The product of a sparse symmetric matrix with a dense matrix (or vector) can also be optimized by specifying the symmetry with selfadjointView():
    288     \code
    289 dm2 = sm1.selfadjointView<>() * dm1;        // if all coefficients of A are stored
    290 dm2 = A.selfadjointView<Upper>() * dm1;     // if only the upper part of A is stored
    291 dm2 = A.selfadjointView<Lower>() * dm1;     // if only the lower part of A is stored
    292     \endcode
    293   - \b sparse-sparse. For sparse-sparse products, two different algorithms are available. The default one is conservative and preserve the explicit zeros that might appear:
    294     \code
    295 sm3 = sm1 * sm2;
    296 sm3 = 4 * sm1.adjoint() * sm2;
    297     \endcode
    298     The second algorithm prunes on the fly the explicit zeros, or the values smaller than a given threshold. It is enabled and controlled through the prune() functions:
    299     \code
    300 sm3 = (sm1 * sm2).pruned();                  // removes numerical zeros
    301 sm3 = (sm1 * sm2).pruned(ref);               // removes elements much smaller than ref
    302 sm3 = (sm1 * sm2).pruned(ref,epsilon);       // removes elements smaller than ref*epsilon
    303     \endcode
    304 
    305   - \b permutations. Finally, permutations can be applied to sparse matrices too:
    306     \code
    307 PermutationMatrix<Dynamic,Dynamic> P = ...;
    308 sm2 = P * sm1;
    309 sm2 = sm1 * P.inverse();
    310 sm2 = sm1.transpose() * P;
    311     \endcode
    312 
    313 
    314 \subsection TutorialSparse_SubMatrices Block operations
    315 
    316 Regarding read-access, sparse matrices expose the same API than for dense matrices to access to sub-matrices such as blocks, columns, and rows. See \ref TutorialBlockOperations for a detailed introduction.
    317 However, for performance reasons, writing to a sub-sparse-matrix is much more limited, and currently only contiguous sets of columns (resp. rows) of a column-major (resp. row-major) SparseMatrix are writable. Moreover, this information has to be known at compile-time, leaving out methods such as <tt>block(...)</tt> and <tt>corner*(...)</tt>. The available API for write-access to a SparseMatrix are summarized below:
    318 \code
    319 SparseMatrix<double,ColMajor> sm1;
    320 sm1.col(j) = ...;
    321 sm1.leftCols(ncols) = ...;
    322 sm1.middleCols(j,ncols) = ...;
    323 sm1.rightCols(ncols) = ...;
    324 
    325 SparseMatrix<double,RowMajor> sm2;
    326 sm2.row(i) = ...;
    327 sm2.topRows(nrows) = ...;
    328 sm2.middleRows(i,nrows) = ...;
    329 sm2.bottomRows(nrows) = ...;
    330 \endcode
    331 
    332 In addition, sparse matrices expose the SparseMatrixBase::innerVector() and SparseMatrixBase::innerVectors() methods, which are aliases to the col/middleCols methods for a column-major storage, and to the row/middleRows methods for a row-major storage.
    333 
    334 \subsection TutorialSparse_TriangularSelfadjoint Triangular and selfadjoint views
    335 
    336 Just as with dense matrices, the triangularView() function can be used to address a triangular part of the matrix, and perform triangular solves with a dense right hand side:
    337 \code
    338 dm2 = sm1.triangularView<Lower>(dm1);
    339 dv2 = sm1.transpose().triangularView<Upper>(dv1);
    340 \endcode
    341 
    342 The selfadjointView() function permits various operations:
    343  - optimized sparse-dense matrix products:
    344     \code
    345 dm2 = sm1.selfadjointView<>() * dm1;        // if all coefficients of A are stored
    346 dm2 = A.selfadjointView<Upper>() * dm1;     // if only the upper part of A is stored
    347 dm2 = A.selfadjointView<Lower>() * dm1;     // if only the lower part of A is stored
    348     \endcode
    349  - copy of triangular parts:
    350     \code
    351 sm2 = sm1.selfadjointView<Upper>();                               // makes a full selfadjoint matrix from the upper triangular part
    352 sm2.selfadjointView<Lower>() = sm1.selfadjointView<Upper>();      // copies the upper triangular part to the lower triangular part
    353     \endcode
    354  - application of symmetric permutations:
    355  \code
    356 PermutationMatrix<Dynamic,Dynamic> P = ...;
    357 sm2 = A.selfadjointView<Upper>().twistedBy(P);                                // compute P S P' from the upper triangular part of A, and make it a full matrix
    358 sm2.selfadjointView<Lower>() = A.selfadjointView<Lower>().twistedBy(P);       // compute P S P' from the lower triangular part of A, and then only compute the lower part
    359  \endcode
    360 
    361 Please, refer to the \link SparseQuickRefPage Quick Reference \endlink  guide for the list of supported operations. The list of linear solvers available is \link TopicSparseSystems here. \endlink
    362 
    363 */
    364 
    365 }
    366