Concatenate sparse matrix Eigen

2
votes

I have two sparse matrices in Eigen, and I would like to join them vertically into one. As an example the target of the code would be:

SparseMatrix<double> matrix1;
matrix1.resize(10, 10);
SparseMatrix<double> matrix2;
matrix2.resize(5, 10);

SparseMatrix<double> MATRIX_JOIN;
MATRIX_JOIN.resize(15, 10);
MATRIX_JOIN << matrix1, matrix2;

I found some solutions on a forum however, I wasn't able to implement it.

What's the proper way to join the matrices vertically?

Edit

My implementation:

SparseMatrix<double> L;
SparseMatrix<double> C;
// ... (Operations with the matrices)
SparseMatrix<double> EMATRIX;
EMATRIX.resize(L.rows() + C.rows(), L.cols());
EMATRIX.middleRows(0, L.rows()) = L;
EMATRIX.middleRows(L.rows(), C.rows()) = C;

I get an error of types, acording to the compiler the right hand side is an Eigen::Block and the left side is Eigen::SparseMatrix

3
c++sparse-matrixeigen
What does "I wasn't able to implement it." mean? - Code-Apprentice
@Code-Apprentice I wasn't able to implement the solution that I found in the forum. It was also meant to join (the solution in the forum) a sparse and a dense matrix - Javier
Repeating the same words doesn't explain what the words mean. Please show exactly what you did and what happens when you compile and run the code you wrote. - Code-Apprentice

3 Answers

5
votes

As far as I know, there is currently no built-in solution. You can be way more efficient than your solution by using the internal insertBack function:

SparseMatrix<double> M(L.rows() + C.rows(), L.cols());
M.reserve(L.nonZeros() + C.nonZeros());
for(Index c=0; c<L.cols(); ++c)
{
    M.startVec(c); // Important: Must be called once for each column before inserting!
    for(SparseMatrix<double>::InnerIterator itL(L, c); itL; ++itL)
         M.insertBack(itL.row(), c) = itL.value();
    for(SparseMatrix<double>::InnerIterator itC(C, c); itC; ++itC)
         M.insertBack(itC.row()+L.rows(), c) = itC.value();
}
M.finalize();
2
votes

Based on @Javier's answer.

  • Fixed the number of columns in the output matrix (just cols instead of cols+cols)
  • Fixed the lower matrice's triplet indices (upper.rows() + it.row() instead of just it.row())

using sparse_matrix_type = Eigen::SparseMatrix<T>;
using triplet_type = Eigen::Triplet<T, size_t>;

static sparse_matrix_type sparse_vstack(sparse_matrix_type const& upper, sparse_matrix_type const& lower) {
  assert(upper.cols() == lower.cols() && "vstack with mismatching number of columns");

  std::vector<triplet_type> triplets;
  triplets.reserve(upper.nonZeros() + lower.nonZeros());

  for (int k = 0; k < upper.outerSize(); ++k) {
    for (sparse_matrix_type::InnerIterator it(upper, k); it; ++it) {
      triplets.emplace_back(it.row(), it.col(), it.value());
    }
  }

  for (int k = 0; k < lower.outerSize(); ++k) {
    for (sparse_matrix_type::InnerIterator it(lower, k); it; ++it) {
      triplets.emplace_back(upper.rows() + it.row(), it.col(), it.value());
    }
  }

  sparse_matrix_type result(lower.rows() + upper.rows(), upper.cols());
  result.setFromTriplets(triplets.begin(), triplets.end());
  return result;
}

Unfortunately I coulnd't get @chtz's example to work with Eigen 3.3.4 due to a static assertion error THIS_SPARSE_BLOCK_SUBEXPRESSION_IS_READ_ONLY. It seems to be explicitly forbidden by Eigen (see https://eigen.tuxfamily.org/dox/SparseBlock_8h_source.html).

1
votes

I ended up doing the following:

MATRIX_JOIN.resize(matrix1.rows() + matrix2.rows(), matrix1.cols() + matrix2.cols());
MATRIX_JOIN.setZero();

// Fill MATRIX_JOIN with triples from the other matrices
std::vector<Triplet<double> > tripletList;
for (int k = 0; k < matrix1.outerSize(); ++k)
{
    for (SparseMatrix<double>::InnerIterator it(matrix1, k); it; ++it)
    {
        tripletList.push_back(Triplet<double>(it.row(), it.col(), it.value()));
    }
}
for (int k = 0; k < matrix2.outerSize(); ++k)
{
    for (SparseMatrix<double>::InnerIterator it(matrix2, k); it; ++it)
    {
        tripletList.push_back(Triplet<double>(it.row(), it.col(), it.value()));
    }
}
FINALMATRIX.setFromTriplets(tripletList.begin(), tripletList.end());

There can be a speedup by calling tripleList.reserve(X) with X being the expected amount of triplets to insert.