稀疏矩阵乘积的特征求解器答案

【问题标题】：Eigen solver for sparse matrix product稀疏矩阵乘积的特征求解器
【发布时间】：2018-11-15 19:25:04
【问题描述】：

我想使用 Eigen 来计算 L_1^{-1}L_2，其中 L_1 和 L_2 都是下三角矩阵，并且在 Eigen 中存储为面向列的稀疏矩阵。我尝试了 Eigen 三角求解器。但是，这需要 L_2 是密集的。

【问题讨论】：

稀疏矩阵的逆矩阵通常是稠密的，因此L_1^{-1}*L_2 也将是稠密的。您确定需要存储矩阵，而不仅仅是动态计算产品/求解系统吗？
@chtz 感谢您的回复。我实际上保证 L_1^{-1} 也是稀疏的。

【解决方案1】：

solve 方法实际上并没有为稀疏 rhs 重载，但是您可以像这样使用solveInPlace 方法（我实际上并没有尝试过）：

Eigen::SparseMatrix<double> foo(Eigen::SparseMatrix<double> const& L1, Eigen::SparseMatrix<double> const& L2)
{
    Eigen::SparseMatrix<double> res = L2;
    L1.triangularView<Eigen::Lower>().solveInPlace(res);
    return res;
}

您仍然应该考虑是否真的需要完整的矩阵。

【讨论】：

感谢您的回复。实际上，我需要计算 spd 矩阵 K = LL^T 的逆，其中 L 非常稀疏。我的解决方案是调用三角求解两次：L_inv = L.triangularView().solve(Identity) 和 L.transpose().triangularView.solve(L_inv)。这是在 Eigen 中正确的做法吗？
你可以计算L_inv.transpose()*L_inv而不是第二次求解
但是如果 L 是稀疏的，L.transpose().triangularView.solve(L_inv) 不比 L_inv.transpose()* L_inv?
我想不出稀疏求解比稀疏乘积更快的好例子。要确定您的情况是否发生这种情况，您需要进行基准测试！
感谢您的帮助！