如何在 C++ 中创建交换矩阵（使用 Eigen 或 Armadillo？）

Question

我想知道如何使用 Eigen 或 Armadillo 库在 C++ 中创建交换矩阵（参见 https://en.wikipedia.org/wiki/Commutation_matrix#MATLAB）。维基页面上有一些 MATLAB 代码：

function P = com_mat(m, n)

% determine permutation applied by K
A = reshape(1:m*n, m, n);
v = reshape(A', 1, []);

% apply this permutation to the rows (i.e. to each column) of identity matrix
P = eye(m*n);
P = P(v,:);

我想知道是否有人在 C++ 中有一个函数来执行此操作或可以将其转换为 C++ 代码？

谢谢

Answer 1

从 Eigen-3.4 开始，您几乎可以将它逐步转换为 Eigen 代码，它引入了

reshaped

方法。

  Eigen::MatrixXf A(3, 2);
  A <<
    1.f, 4.f,
    2.f, 5.f,
    3.f, 6.f;
  std::cout << A << "\n\n";
  /*
   * output:
   * 1 4
   * 2 5
   * 3 6
   */
  using PermutationMatrixXi =
    Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic, int>;
  PermutationMatrixXi p(A.size());
  p.indices() = Eigen::VectorXi::LinSpaced(A.size(), 0, A.size() - 1)
        .reshaped(A.cols(), A.rows()).transpose().reshaped(A.size(), 1);
  std::cout << p.indices().transpose() << "\n\n";
  /*
   * output: 0 2 4 1 3 5
   */
  Eigen::VectorXf pa = p * A.reshaped(A.size(), 1);
  std::cout << pa.transpose() << "\n\n";
  /*
   * output: 1 4 2 5 3 6
   */
  Eigen::MatrixXi P = p.toDenseMatrix();
  std::cout << P << "\n\n";
  /*
   * output:
   * 1 0 0 0 0 0
   * 0 0 0 1 0 0
   * 0 1 0 0 0 0
   * 0 0 0 0 1 0
   * 0 0 1 0 0 0
   * 0 0 0 0 0 1
   */
  Eigen::VectorXf pa = p * A.reshaped(A.size(), 1);
  std::cout << pa.transpose() << "\n\n";
  /*
   * output: 1 4 2 5 3 6
   */

请注意，

reshaped

是一个相当慢的操作，因为它的索引操作涉及整数除法和取模。如果您需要速度，您可能需要手动构建排列，这不太难，并使用

Eigen::Map

重塑矩阵。

  Eigen::MatrixXf A(3, 2);
  [...]
  PermutationMatrixXi p(A.size());
  Eigen::VectorXi& pi = p.indices();
  for(int col = 0, cols = A.cols(); col < cols; ++col)
    for(int row = 0, rows = A.rows(); row < rows; ++row)
      pi[col * rows + row] = row * cols + col;
  assert(A.innerStride() == 1 && A.outerStride() == A.rows());
  Eigen::VectorXf pa = p * Eigen::VectorXf::Map(A.data(), A.size());

注意断言：

Map

仅适用于密集矩阵或列块；不是更大矩阵的任意块。如果需要，评估成一个新的密集矩阵。