Assemble block sparse matrix in Eigen

Question

Lets say I have a small sparse matrix B. I want to build a bigger sparse matrix like

BtB      0
0     (BtB)^-1

I want to know if Eigen provides some functionality to assemble something like this. I have been searching and found nothing. One option that I can use is to compute the operations, extract the triplets and assemble the matrix based on the triplets. Is there any easier way?

Answer 1

If you want to explicitly build a sparse block matrix (there are legitimate reasons!), this is what I use (only works for column major -- easily adapted to row major). Here, sparse_stack_v means to "stack" two matrices vertically. sparse_stack_his likewise used to "stack" two sparse matrices horizontally. sparse_stack_h_inplace is a more efficient operation as you can reuse most of the sparsity structure for the left matrix (and overwrite it).

void sparse_stack_v(
    const SparseMatrix<Scalar, ColMajor, StorageIndex>& top,
    const SparseMatrix<Scalar, ColMajor, StorageIndex>& bottom,
    SparseMatrix<Scalar, ColMajor, StorageIndex>& stacked)
{
    assert(top.cols() == bottom.cols());
    stacked.resize(top.rows() + bottom.rows(), top.cols());
    stacked.resizeNonZeros(top.nonZeros() + bottom.nonZeros());

    StorageIndex i = 0;

    for (StorageIndex col = 0; col < top.cols(); col++)
    {
        stacked.outerIndexPtr()[col] = i;

        for (StorageIndex j = top.outerIndexPtr()[col]; j < top.outerIndexPtr()[col + 1]; j++, i++)
        {
            stacked.innerIndexPtr()[i] = top.innerIndexPtr()[j];
            stacked.valuePtr()[i] = top.valuePtr()[j];
        }

        for (StorageIndex j = bottom.outerIndexPtr()[col]; j < bottom.outerIndexPtr()[col + 1]; j++, i++)
        {
            stacked.innerIndexPtr()[i] = (StorageIndex)top.rows() + bottom.innerIndexPtr()[j];
            stacked.valuePtr()[i] = bottom.valuePtr()[j];
        }
    }
    stacked.outerIndexPtr()[top.cols()] = i;
    assert(stacked.isCompressed());
}

template<typename Scalar, typename StorageIndex>
void sparse_stack_h(
    const SparseMatrix<Scalar, ColMajor, StorageIndex>& left,
    const SparseMatrix<Scalar, ColMajor, StorageIndex>& right,
    SparseMatrix<Scalar, ColMajor, StorageIndex>& stacked)
{
    assert(left.rows() == right.rows());

    stacked.resize(left.rows(), left.cols() + right.cols());
    stacked.resizeNonZeros(left.nonZeros() + right.nonZeros());

    std::copy(left.innerIndexPtr(), left.innerIndexPtr() + left.nonZeros(), stacked.innerIndexPtr());
    std::copy(right.innerIndexPtr(), right.innerIndexPtr() + right.nonZeros(), stacked.innerIndexPtr() + left.nonZeros());

    std::copy(left.valuePtr(), left.valuePtr() + left.nonZeros(), stacked.valuePtr());
    std::copy(right.valuePtr(), right.valuePtr() + right.nonZeros(), stacked.valuePtr() + left.nonZeros());

    std::copy(left.outerIndexPtr(), left.outerIndexPtr() + left.cols(), stacked.outerIndexPtr());//dont need the last entry of A.outerIndexPtr() -- total length is AB.cols() + 1 = A.cols() + B.cols() + 1
    std::transform(right.outerIndexPtr(), right.outerIndexPtr() + right.cols() + 1, stacked.outerIndexPtr() + left.cols(), [&](StorageIndex i) { return i + left.nonZeros(); });

    assert(stacked.isCompressed());
}

template<typename Scalar, typename StorageIndex>
void sparse_stack_h_inplace(
    SparseMatrix<Scalar, ColMajor, StorageIndex>& left,
    const SparseMatrix<Scalar, ColMajor, StorageIndex>& right)
{
    assert(left.rows() == right.rows());

    const StorageIndex leftcol = (StorageIndex)left.cols();
    const StorageIndex leftnz = (StorageIndex)left.nonZeros();

    left.conservativeResize(left.rows(), left.cols() + right.cols());
    left.resizeNonZeros(left.nonZeros() + right.nonZeros());

    std::copy(right.innerIndexPtr(), right.innerIndexPtr() + right.nonZeros(), left.innerIndexPtr() + leftnz);
    std::copy(right.valuePtr(), right.valuePtr() + right.nonZeros(), left.valuePtr() + leftnz);
    std::transform(right.outerIndexPtr(), right.outerIndexPtr() + right.cols() + 1, left.outerIndexPtr() + leftcol, [&](StorageIndex i) { return i + leftnz; });

    assert(left.isCompressed());
}

If you want to 'pad' a sparse matrix with zeros, just stack the matrices with an empty sparse matrix of an appropriate size.

Answer 2

I know this is years late, but in case anyone else finds their way here, I have a solution. My particular implementation comes from the perspective of using Rcpp and RcppEigen, but I believe this could be written using the std library list object class. Since I'm working with precision matrices, I do assume that the input matrices are square, but again, it would not take too much thought to convert this code to allow for arbitrary dimensions from each of the constituent matrices. This will create a block-diagonal sparse matrix efficiently:

#include <Rcpp.h>
#include <RcppEigen.h>

Eigen::SparseMatrix<double> sparseBdiag(Rcpp::List B_list)
{
  int K = B_list.length();
   Eigen::VectorXi B_cols(K);
  for(int k=0;k<K;k++) {
    Eigen::SparseMatrix<double> Bk = B_list(k);
    B_cols[k] = Bk.cols();
  }
  int sumCols = B_cols.sum();
  Eigen::SparseMatrix<double> A(sumCols,sumCols);
  int startCol=0, stopCol=0, Bk_cols;
  for(int k=0;k<K;k++) {
    Eigen::SparseMatrix<double> Bk = B_list(k);
    Bk_cols = Bk.cols();
    stopCol = startCol + Bk_cols;
    for(int j=startCol;j<stopCol;j++){
      A.startVec(j);
      for(Eigen::SparseMatrix<double,0,int>::InnerIterator it(Bk,j-startCol); it; ++it) {
        A.insertBack(it.row()+startCol,j) = it.value();
      }
    }
    startCol = stopCol;
  }
  A.finalize();
  return A;
}