Extract a block from a sparse matrix as another sparse matric

Question

How to extract a block from a Eigen::SparseMatrix<double>. It seems there aren't the methods I used for the dense ones.

‘class Eigen::SparseMatrix<double>’ has no member named ‘topLeftCorner’
‘class Eigen::SparseMatrix<double>’ has no member named ‘block’

There is a way to extract a block as a Eigen::SparseMatrix<double> ?

Answer 1

I made this function to extract blocks from a Eigen::SparseMatrix<double,ColMaior>

typedef Triplet<double> Tri;
SparseMatrix<double> sparseBlock(SparseMatrix<double,ColMajor> M,
        int ibegin, int jbegin, int icount, int jcount){
        //only for ColMajor Sparse Matrix
    assert(ibegin+icount <= M.rows());
    assert(jbegin+jcount <= M.cols());
    int Mj,Mi,i,j,currOuterIndex,nextOuterIndex;
    vector<Tri> tripletList;
    tripletList.reserve(M.nonZeros());

    for(j=0; j<jcount; j++){
        Mj=j+jbegin;
        currOuterIndex = M.outerIndexPtr()[Mj];
        nextOuterIndex = M.outerIndexPtr()[Mj+1];

        for(int a = currOuterIndex; a<nextOuterIndex; a++){
            Mi=M.innerIndexPtr()[a];

            if(Mi < ibegin) continue;
            if(Mi >= ibegin + icount) break;

            i=Mi-ibegin;    
            tripletList.push_back(Tri(i,j,M.valuePtr()[a]));
        }
    }
    SparseMatrix<double> matS(icount,jcount);
    matS.setFromTriplets(tripletList.begin(), tripletList.end());
    return matS;
}

And these if the sub-matrix is in one of the four corners:

SparseMatrix<double> sparseTopLeftBlock(SparseMatrix<double> M,
        int icount, int jcount){
    return sparseBlock(M,0,0,icount,jcount);
}
SparseMatrix<double> sparseTopRightBlock(SparseMatrix<double> M,
        int icount, int jcount){
    return sparseBlock(M,0,M.cols()-jcount,icount,jcount);
}
SparseMatrix<double> sparseBottomLeftBlock(SparseMatrix<double> M,
        int icount, int jcount){
    return sparseBlock(M,M.rows()-icount,0,icount,jcount);
}
SparseMatrix<double> sparseBottomRightBlock(SparseMatrix<double> M,
        int icount, int jcount){
    return sparseBlock(M,M.rows()-icount,M.cols()-jcount,icount,jcount);
}

Answer 2

This is now supported in Eigen 3.2.2 Docs (though maybe earlier versions support it too).

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Sparse>

using namespace Eigen;

int main()
{
  MatrixXd silly(6, 3);

  silly << 0, 1, 2,
           0, 3, 0,
           2, 0, 0,
           3, 2, 1,
           0, 1, 0,
           2, 0, 0;



  SparseMatrix<double, RowMajor> sparse_silly = silly.sparseView();

  std::cout <<"Whole Matrix" << std::endl;
  std::cout << sparse_silly << std::endl;

  std::cout << "block of matrix" << std::endl;
  std::cout << sparse_silly.block(1,1,3,2) << std::endl;

  return 0;
}

Answer 3

There is very sparse support (sorry, no pun intended) for submatrices in sparse matrices. Effectively you can only access a continuous set of rows for row-major, and columns for column major. The reason for that is not that the matrices could be empty, but rather that the indexing scheme is somewhat more complicated than with dense matrices. With dense matrices you only need an additional stride number in order to support sub-matrix support.