连接稀疏矩阵特征答案

【问题标题】：Concatenate sparse matrix Eigen连接稀疏矩阵特征
【发布时间】：2017-01-20 05:02:24
【问题描述】：

我在 Eigen 中有两个稀疏矩阵，我想将它们垂直合并为一个。例如，代码的目标是：

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;

我在forum 上找到了一些解决方案，但是我无法实施。

垂直连接矩阵的正确方法是什么？

编辑

我的实现：

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;

我收到类型错误，根据编译器，右侧是 Eigen::Block，左侧是 Eigen::SparseMatrix

【问题讨论】：

什么是“我无法实现它。”什么意思？
@Code-Apprentice 我无法实施我在论坛中找到的解决方案。它也意味着加入（论坛中的解决方案）一个稀疏和密集的矩阵
重复相同的词并不能解释这些词的含义。请准确说明您做了什么以及编译和运行您编写的代码时会发生什么。

标签： c++ sparse-matrix eigen

【解决方案1】：

据我所知，目前没有内置解决方案。通过使用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();

【讨论】：

这比使用三元组更有效吗？
@Javier 它应该在所有非常奇怪的情况下。在许多情况下，这可能并不重要。
@Pippo 感谢您纠正这个错误。我注意到另一个（我实际上从未测试过这段代码）

【解决方案2】：

基于@Javier 的回答。

修复了输出矩阵中的列数（只是cols 而不是cols+cols）
修复了下矩阵的三元组索引（upper.rows() + it.row() 而不仅仅是 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;
}

不幸的是，由于静态断言错误THIS_SPARSE_BLOCK_SUBEXPRESSION_IS_READ_ONLY，我无法让@chtz 的示例与 Eigen 3.3.4 一起使用。它似乎被 Eigen 明确禁止（参见https://eigen.tuxfamily.org/dox/SparseBlock_8h_source.html）。

【讨论】：

我认为你这里还有一个bug：当你创建result时，upper.rows() + upper.rows()应该是upper.rows() + lower.rows()我认为。
我的解决方案中没有任何类型的block。您是指@Javier 在他的问题中的尝试吗？这确实行不通。

【解决方案3】：

我最终做了以下事情：

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());

调用tripleList.reserve(X) 可以加快速度，其中 X 是要插入的三连音的预期数量。

【讨论】：