如何在c++中创建交换矩阵(使用Eigen或Armadillo?)

Question

我想知道如何使用Eigen或Armadillo库在c++中创建交换矩阵(参见https://en.wikipedia.org/wiki/Commutation_matrix#MATLAB)。wiki页面上有一些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

您可以从引入reshaped方法的Eigen-3.4开始逐步将其转换为特征代码。

Eigen::MatrixXf A(3, 2);
A <<
1.f, 4.f,
2.f, 5.f,
3.f, 6.f;
std::cout << A << "nn";
/*
* 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() << "nn";
/*
* output: 0 2 4 1 3 5
*/
Eigen::VectorXf pa = p * A.reshaped(A.size(), 1);
std::cout << pa.transpose() << "nn";
/*
* output: 1 4 2 5 3 6
*/
Eigen::MatrixXi P = p.toDenseMatrix();
std::cout << P << "nn";
/*
* 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() << "nn";
/*
* 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只适用于密集矩阵或列块;而不是大矩阵的任意块。如果需要，求值成一个新的密集矩阵。