如何在特征稀疏矩阵中执行逐元素加法

Question

我有以下代码，其中 AMAT 当前是一个密集矩阵。然而，大多数元素都是零，因此本质上它是一个稀疏矩阵。我了解 Eigen 稀疏矩阵不支持块操作。想知道如果我将 AMAT 替换为稀疏矩阵，我该如何重写这段代码。 BMAT 是一个 9x9 的密集矩阵，BMAT 的每个 3x3 块都被添加到 AMAT 中的特定块中。 BMAT 在此循环外计算。

    for(j=0;j<5000;j++) {

      id1 = ids(0,j);
      id2 = ids(1,j);
      id3 = ids(2,j);

      AMAT.block<3,3>(id1*3,id1*3) = AMAT.block<3,3>(id1*3,id1*3) + BMAT.block<3,3>(0,0);
      AMAT.block<3,3>(id1*3,id2*3) = AMAT.block<3,3>(id1*3,id2*3) + BMAT.block<3,3>(0,3);
      AMAT.block<3,3>(id1*3,id3*3) = AMAT.block<3,3>(id1*3,id3*3) + BMAT.block<3,3>(0,6);

      AMAT.block<3,3>(id2*3,id1*3) = AMAT.block<3,3>(id2*3,id1*3) + BMAT.block<3,3>(3,0);
      AMAT.block<3,3>(id2*3,id2*3) = AMAT.block<3,3>(id2*3,id2*3) + BMAT.block<3,3>(3,3);
      AMAT.block<3,3>(id2*3,id3*3) = AMAT.block<3,3>(id2*3,id3*3) + BMAT.block<3,3>(3,6);

      AMAT.block<3,3>(id3*3,id1*3) = AMAT.block<3,3>(id3*3,id1*3) + BMAT.block<3,3>(6,0);
      AMAT.block<3,3>(id3*3,id2*3) = AMAT.block<3,3>(id3*3,id2*3) + BMAT.block<3,3>(6,3);
      AMAT.block<3,3>(id3*3,id3*3) = AMAT.block<3,3>(id3*3,id3*3) + BMAT.block<3,3>(6,6);

   }

Answer 1

据我了解 https://eigen.tuxfamily.org/dox/group__TutorialBlockOperations.html 中的本教程，块操作是可能的。但是你需要在编译时知道列和行。

Answer 2

这可行（未经测试，我不知道您的矩阵的实际类型）。这个想法是编写一个自定义迭代器，它提供 AMAT 的每个条目的索引和值，并将其传递给 setFromTriplets （重复的条目将被加在一起）。这将在您的索引列表中迭代两次，但遗憾的是不会利用 AMAT 的块结构。但它会在 O(nnz) 时间内执行。

#include <Eigen/SparseCore>

struct AMAT_constructor {
  struct AMAT_iterator {
    bool operator==(AMAT_iterator const& other) const {
      return j == other.j && k == other.k;
    }

    bool operator!=(AMAT_iterator const& other) const {
      return !(*this == other);
    }

    Eigen::Index operator-(AMAT_iterator const& other) const {
      return (j - other.j) * 81 + k - other.k;
    }

    AMAT_iterator const* operator->() const { return this; }
    AMAT_iterator const& operator*() const { return *this; }
    float value() const { return BMAT(k); }
    Eigen::Index row() const { return ids((k / 3) % 3, j) * 3 + k % 3; }

    Eigen::Index col() const { return ids(k / 27, j) * 3 + (k / 9) % 3; }

    AMAT_iterator& operator++() {
      if (++k == 81) {
        k = 0;
        ++j;
      }
      return *this;
    }

    Eigen::Index j, k;
    Eigen::Matrix3Xi const& ids;
    Eigen::Matrix<float, 9, 9> const& BMAT;
  };

  Eigen::Matrix3Xi const& ids;
  Eigen::Matrix<float, 9, 9> const& BMAT;

  AMAT_iterator begin() const { return AMAT_iterator{0, 0, ids, BMAT}; }
  AMAT_iterator end() const { return AMAT_iterator{ids.cols(), 0, ids, BMAT}; }
};

// use it like this:
Eigen::SparseMatrix<float> foo(Eigen::Matrix3Xi const& ids,
                               Eigen::Matrix<float, 9, 9> const& BMAT,
                               Eigen::Index sizeA) {

  Eigen::SparseMatrix<float> AMAT(sizeA, sizeA);
  AMAT_constructor Ac{ids, BMAT};
  AMAT.setFromTriplets(Ac.begin(), Ac.end());

  return AMAT;
}

如何在特征稀疏矩阵中执行逐元素加法

How to perform an elementwise addition in an Eigen sparse matrix

c++

sparse-matrix

eigen