Изменить / уменьшить собственную перестановочную матрицу

У меня проблемы с решением того, что я считаю довольно простой проблемой. Основная проблема в том, что я хочу изменить Eigen PermutationMatrix, но я не знаю как.

Я делаю QR-разложение какой-то матрицы X используя библиотеку C ++ Eigen. Я делаю это на матрицах с недостатком ранга и мне нужен конкретный вывод. Конкретно мне нужно

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

Проблема в том, что с помощью Eigen::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 R_inv, из которых я извлек 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;