3

我无法解决我认为应该是一个相当简单的问题。基本问题是我想修改 Eigen PermutationMatrix 但我不知道如何。

我正在X使用 C++ Eigen 库对一些矩阵进行 QR 分解。我在秩不足的矩阵上做这个，我需要一个特定的输出。具体来说，我需要

R^{-1} * t(R^{-1})

问题是 usingEigen::ColPivHouseholderQR返回一个置换版本的R. 这很容易解决X满级时的问题，但我想要最快的解决方案来解决等级不足的情况。让我演示一下：

using namespace Eigen;

// Do QR 
ColPivHouseholderQR<MatrixXd> PQR(X);
// Get permutation matrix
ColPivHouseholderQR<MatrixXd>::PermutationType Pmat(PQR.colsPermutation());
int r(PQR.rank());
int p(X.cols());

// Notice I'm only getting r elements, so R_inv is r by r
MatrixXd R_inv = PQR.matrixQR().topLeftCorner(r, r).triangularView<Upper>().solve(MatrixXd::Identity(r, r));

// This only works if r = p and X is full-rank
R_inv = Pmat * R_inv * Pmat;
XtX_inv = R_inv * R_inv.transpose();

所以基本问题是我想修改 Pmat 以便它只置换我从中提取的 R_inv 的 r 列PQR.matrixQR()。我的基本问题是我不知道如何使用 Eigen PermutationMatrix 修改工作，因为它似乎没有任何普通矩阵的方法或属性。

一种可能的解决方案如下：当我相乘时Pmat * MatrixXd::Identity(p, p)，我得到一个有用的矩阵。

例如，我得到类似的东西：

[0, 1, 0, 0,
 1, 0, 0, 0,
 0, 0, 0, 1,
 0, 0, 1, 0]

如果 p = 4 和 r = 3，那么我会喜欢这个子视图，在其中我删除前 r 列右侧的所有列，然后删除所有为 0 的行：

[0, 1, 0,
 1, 0, 0,
 0, 0, 1]

所以我可以做到以下几点：

P = Pmat * MatrixXd::Identity(p, p)
P.leftCols(p);
MatrixXd P = Pmat * Eigen::MatrixXd::Identity(p, p);

// https://stackoverflow.com/questions/41305178/removing-zero-columns-or-rows-using-eigen
// find non-zero columns:
Matrix<bool, 1, Dynamic> non_zeros = P.leftCols(r).cast<bool>().rowwise().any();

// allocate result matrix:
MatrixXd res(non_zeros.count(), r);

// fill result matrix:
Index j=0;
for(Index i=0; i<P.rows(); ++i)
{
  if(non_zeros(i))
    res.row(j++) = P.row(i).leftCols(r);
}

R_inv = res * R_inv * res;

XtX_inv = R_inv * R_inv.transpose();

但这似乎很昂贵，并且没有利用 Pmat 已经知道应该删除哪些 Pmat 行的事实。我猜有一种更简单的方法可以使用 Pmat。

有什么方法可以轻松修改 Eigen PermutationMatrix 以仅考虑未放置在第一个r位置之外的列？

任何帮助或提示将不胜感激。

我想出了另一个解决方案，它可能需要更少的计算。

// Get all column indices
ArrayXi Pmat_indices = Pmat.indices();
// Get the order for the columns you are keeping
ArrayXi Pmat_keep = Pmat_indices.head(r);
// Get the indices for columns you are discarding
ArrayXi Pmat_toss = Pmat_indices.tail(p - r);

// this code takes the indices you are keeping, and, while preserving order, keeps them in the range [0, r-1]
// For each one, see how many dropped indices are smaller, and subtract that difference
// Ex: p = 4, r = 2
// Pmat_indices = {3, 1, 0, 2}
// Pmat_keep = {3, 1}
// Pmat_toss = {0, 2}
// Now we go through each keeper, count how many in toss are smaller, and then modify accordingly
// 3 - 2 and 1 - 1
// Pmat_keep = {1, 0}
for(Index i=0; i<r; ++i)
{
  Pmat_keep(i) = Pmat_keep(i) - (Pmat_toss < Pmat_keep(i)).count();
}

// Now this will order just the first few columns in the right order
PermutationMatrix<Dynamic, Dynamic> P = PermutationWrapper<ArrayXi>(Pmat_keep);

R_inv = P * R_inv * P;
4

0 回答 0