如何在Eigen库中利用稀疏分解提取稀疏矩阵的行基？

行空间的基对应矩阵中一组线性无关的行，要获取这些行的索引，可以借助Eigen的稀疏分解方法，以下两种方案可行：

方案一：利用SparseQR分解（基于矩阵转置） #

原矩阵A的行空间等价于其转置矩阵A^T的列空间。通过对A^T做稀疏QR分解，可找到A^T列空间基对应的列索引——这些索引就是原矩阵A的行基索引。

Eigen的SparseQR分解形式为：A^T = Q * R * P，其中P是列置换矩阵，它会将A^T的列重新排序，使得R的前rank个主对角线元素非零（rank为矩阵的秩）。P的逆置换中前rank个位置对应的索引，就是A^T列空间基的列索引，即原矩阵A的行基索引。

代码示例：

#include <Eigen/SparseQR>
#include <Eigen/SparseCore>
#include <vector>

int main() {
    Eigen::SparseMatrix<double> A;
    // 假设A已完成初始化

    Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> qr;
    qr.compute(A.transpose()); // 对A的转置执行QR分解
    int rank = qr.rank();
    const auto& permutation = qr.colsPermutation(); // 获取列置换矩阵P

    // 提取行基索引：取P逆置换的前rank个元素
    std::vector<int> row_basis_indices;
    row_basis_indices.reserve(rank);
    for (int i = 0; i < rank; ++i) {
        row_basis_indices.push_back(permutation.inverse()(i));
    }

    // row_basis_indices即为构成行空间基的行索引
    return 0;
}

方案二：利用SparseLU分解（直接处理原矩阵） #

Eigen的SparseLU会对原矩阵做带行/列置换的LU分解：A = P * L * U * Q，其中P是行置换矩阵，U是上三角矩阵。U的非零行数等于矩阵的秩，这些行对应P*A中的线性无关行，因此原矩阵A中对应的行就是P置换后的前rank个行索引。

代码示例：

#include <Eigen/SparseLU>
#include <Eigen/SparseCore>
#include <vector>

int main() {
    Eigen::SparseMatrix<double> A;
    // 假设A已完成初始化

    Eigen::SparseLU<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> lu;
    lu.compute(A);
    int rank = lu.rank();
    const auto& row_perm = lu.rowPermutation(); // 获取行置换矩阵P

    // 提取行基索引：取row_perm前rank个元素，对应原矩阵的行索引
    std::vector<int> row_basis_indices;
    row_basis_indices.reserve(rank);
    for (int i = 0; i < rank; ++i) {
        row_basis_indices.push_back(row_perm.indices()(i));
    }

    // row_basis_indices即为构成行空间基的行索引
    return 0;
}

注意事项 #

• 两种方案得到的行基索引可能不同，但对应的行都能张成原矩阵的行空间（行空间的基不唯一）。
• 矩阵的秩计算依赖于分解过程中的阈值，若需自定义阈值，可通过qr.setPivotThreshold(threshold)或lu.setPivotThreshold(threshold)设置。

内容的提问来源于stack exchange，提问作者Pew