问题描述

在我想要实现的算法中，我需要对一个矩阵进行重复操作，该矩阵由从另一个现有矩阵中提取的一组行组成。为清楚起见，我将“ A”称为原始矩阵，而“ B”称为“工作矩阵”（由的某些行组成A）。

算法进行如下：

r从中选择一行A
插入或删除这样的行B
执行一些操作B（如果它不为空）
重复

一些相关信息：

如果r必须插入B，我确信它不存在（没有创建重复项的风险）
如果r必须从中删除B，我确信它已经存在（删除不存在的行没有风险）
内部行的顺序B无关紧要

我的目标是使用 Eigen 有效地实现这部分算法。效率方面，我主要指的是执行时间。

我想出了 3 种不同的解决方案，我想知道您对它们的看法以及可以添加哪些改进。

解决方案1：调整大小并复制

这个想法是存储B在一个MatrixXd反复扩展/收缩的实例中。由于行的顺序无关紧要，我将新行堆叠在B. 当我需要从中删除一行时B，我只需向上复制较低的行。

class Algorithm1 {
public:
  bool show_b;

  Algorithm1(const Eigen::MatrixXd& M) : A(M), B(0,M.cols()), show_b(false) {
    rows.reserve(A.rows());
  }

  void run(const std::vector<std::pair<bool,int>>& ops) {
    for(const auto& op : ops) {
      if(rows.size() > 0) {
        // If B is not empty, do something with it (print is just an example)
        if(show_b)
          std::cout << "-------------------" << std::endl << "B is:" << std::endl << B << std::endl;
      }

      // check if we have to insert or remove the row
      if(op.first) {
        // Row 'op.second' should be added to B.
        B.conservativeResize(rows.size()+1, A.cols()); // add space for a column at the bottom
        B.row(rows.size()) = A.row(op.second); // copy the desired row from A to B
        rows.push_back(op.second); // keep track of which rows are inside B
      }
      else {
        // Row 'op.second' should be removed from B.
        auto it = std::find(rows.begin(), rows.end(), op.second); // find where the row is stored in B
        int r = std::distance(rows.begin(), it); // get the actual index
        for(int i=r; i<rows.size()-1; i++)
          B.row(i) = B.row(i+1); // copy each row to the one above
        B.conservativeResize(rows.size()-1, A.cols()); // remove the last row of B
        rows.erase(it); // remove the row also from this list
      }
    }
  }

private:
  Eigen::MatrixXd A; // Original matrix
  Eigen::MatrixXd B; // Working matrix
  std::vector<int> rows; // list of rows in B, such that A.row(rows[i]) equals B.row(i)
};

请注意，原则上要添加/删除的行是由算法在内部选择的。但是，为了简单起见，我手动将“操作列表”作为参数传递给算法。的每个元素std::vector都是 a pair，其第一个元素表明该行是必须添加 ( pair.firstis true) 还是删除 ( pair.firstis false)。要选择的行pair.second的索引在（索引是指的行A）中给出。

解决方案2：`Eigen::Map`

我能想到的另一个解决方案是将的数据存储B在“原始缓冲区”中，可以通过Eigen::Map. 为了能够从缓冲区中快速添加/删除整行，我必须对和都使用行优先A排序B。

class Algorithm2 {
public:
  // Row-major matrix
  typedef Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic,Eigen::RowMajor> RMatrixXd;

  bool show_b;

  Algorithm2(const Eigen::MatrixXd& M) : A(M), show_b(false) {
    rows.reserve(A.rows());
    data.reserve(A.size());
  }

  void run(const std::vector<std::pair<bool,int>>& ops) {
    for(const auto& op : ops) {
      if(rows.size() > 0) {
        // Create the sub-matrix B using a Map
        Eigen::Map<RMatrixXd> B(data.data(), rows.size(), A.cols());

        // Do something with it (print is just an example)
        if(show_b)
          std::cout << "-------------------" << std::endl << "B is:" << std::endl << B << std::endl;
      }

      // check if we have to insert or remove the row
      if(op.first) {
        // Row 'op.second' should be added to B.
        // append the row content inside the data
        data.insert(
          data.end(), // append at the end
          A.data() + op.second*A.cols(), // from the beginning of the row
          A.data() + (op.second+1)*A.cols() // until the beginning of the next one
        );
        rows.push_back(op.second); // keep track of which rows are inside B
      }
      else {
        // Row 'op.second' should be removed from B.
        auto it = std::find(rows.begin(), rows.end(), op.second); // find where the row is stored in B
        int r = std::distance(rows.begin(), it); // get the actual index
        data.erase(data.begin()+r*A.cols(), data.begin()+(r+1)*A.cols()); // remove the corresponding values from the data
        rows.erase(it); // remove the row also from this list
      }
    }
  }

private:
  RMatrixXd A; // The original matrix
  std::vector<int> rows; // list of rows in B, such that A.row(rows[i]) equals B.row(i)
  std::vector<double> data; // content of B (row-major ordering)
};

解决方案3：新的重载`operator()`

使用最新版本的 Eigen，可以直接通过选择子矩阵operator()。在实践中，矩阵B可以得到为A( rows-to-select , Eigen::all )。

class Algorithm3 {
public:
  bool show_b;

  Algorithm3(const Eigen::MatrixXd& M) : A(M), show_b(false) {
    rows.reserve(A.rows());
  }

  void run(const std::vector<std::pair<bool,int>>& ops) {
    for(const auto& op : ops) {
      if(rows.size() > 0) {
        // Create the sub-matrix B
        auto B = A(rows, Eigen::all);

        // Do something with it (print is just an example)
        if(show_b)
          std::cout << "-------------------" << std::endl << "B is:" << std::endl << B << std::endl;
      }

      // check if we have to insert or remove the row
      if(op.first) {
        // Row 'op.second' should be added to B.
        rows.push_back(op.second); // keep track of which rows are inside B
      }
      else {
        // Row 'op.second' should be removed from B.
        rows.erase(std::find(rows.begin(), rows.end(), op.second)); // remove the row from the list
      }
    }
  }

private:
  Eigen::MatrixXd A; // The original matrix
  std::vector<int> rows; // list of rows in B, such that A.row(rows[i]) equals B.row(i)
};

解决方案的快速比较

下面是一个记录执行时间的简单程序：

int main() {

  Eigen::MatrixXd A(3,3);
  A << 0,0,0,1,1,1,2,2,2;

  std::vector<std::pair<bool,int>> ops {
    // {add=true/remove=false, row-index} // which rows will now be inside B
    {true, 1}, // 1
    {true, 0}, // 0 1
    {true, 2}, // 0 1 2
    {false,0}, // 1 2
    {false,2}, // 1
    {true, 2}, // 1 2
    {false,1}, // 2
    {true, 0}, // 0 2
    {true, 1}, // 0 1 2
    {false,0}, // 1 2
    {false,1}, // 2
    {false,2}, // empty
  };

  // Instanciate
  Algorithm1 alg1(A);
  Algorithm2 alg2(A);
  Algorithm3 alg3(A);

  // warm-up
  for(int i=0; i<1000; i++) {
    alg1.run(ops);
    alg2.run(ops);
    alg3.run(ops);
  }

  // run the 3 algorithms
  auto start1 = std::chrono::high_resolution_clock::now();
  for(int i=0; i<1000; i++)
    alg1.run(ops);
  auto end1 = std::chrono::high_resolution_clock::now();

  auto start2 = std::chrono::high_resolution_clock::now();
  for(int i=0; i<1000; i++)
    alg2.run(ops);
  auto end2 = std::chrono::high_resolution_clock::now();

  auto start3 = std::chrono::high_resolution_clock::now();
  for(int i=0; i<1000; i++)
    alg3.run(ops);
  auto end3 = std::chrono::high_resolution_clock::now();

  // show the results
  std::cout << "Alg1 took " << std::chrono::duration_cast<std::chrono::nanoseconds>(end1-start1).count() << "ns" << std::endl;
  std::cout << "Alg2 took " << std::chrono::duration_cast<std::chrono::nanoseconds>(end2-start2).count() << "ns" << std::endl;
  std::cout << "Alg3 took " << std::chrono::duration_cast<std::chrono::nanoseconds>(end3-start3).count() << "ns" << std::endl;

  return 0;
}

在我的笔记本电脑上，第二个策略给出了最好的结果，其次是第三个（约 3.4 倍）和第一个（约 4.3 倍慢）。

你怎么看？与您的经验相比，结果是否有意义？你能提出更好的策略吗？

标签： c++eigen3

c++ - Eigen：在子矩阵上重复操作的有效方法

问题描述

解决方案1：调整大小并复制

解决方案2：`Eigen::Map`

解决方案3：新的重载`operator()`

解决方案的快速比较

解决方案

推荐阅读

c++ - Eigen：在子矩阵上重复操作的有效方法

问题描述

解决方案1：调整大小并复制

解决方案2：Eigen::Map

解决方案3：新的重载operator()

解决方案的快速比较

解决方案

推荐阅读

解决方案2：`Eigen::Map`

解决方案3：新的重载`operator()`