基于Eigen的稀疏SPD矩阵分解后求解Lx=b及Px=b的效率问询

Great question! Let's break this down clearly for Eigen's LLt and LDLt solvers for sparse SPD matrices:

1. 求解Lx=b的效率问题 #

Short answer: Yes, this approach is absolutely efficient—and it's exactly how you should do it.

Here's why:

• After calling solver.compute(A), the solver has already computed the sparse L factor (and D for LDLt) along with the permutation matrix P. The solver.matrixL() method returns a direct view of this precomputed lower triangular matrix, no extra computation needed.
• Wrapping it with TriangularView<Lower>() tells Eigen to use its optimized triangular matrix solver, which leverages the sparse structure of L to perform the solve in O(nnz(L)) time (where nnz is the number of non-zero elements). This is the optimal complexity for solving a triangular system with a precomputed factor.
• There's no redundant overhead here—Eigen doesn't reprocess the matrix or copy data unnecessarily; it uses the existing decomposition data directly.

Just a quick note: For LLt, L is a lower triangular matrix with non-unit diagonal entries, while for LDLt, L is a unit lower triangular matrix (diagonal entries are 1). The TriangularView<Lower>() handles both cases correctly, so your code works for both solvers.

2. 求解Px=b的操作方法 #

To solve the permutation system Px = b, you need to compute x = P⁻¹b (since multiplying both sides by P⁻¹ gives x). Eigen makes this straightforward with the solver's built-in permutation accessor:

// Get the permutation matrix from the solver
const auto& P = solver.permutation();
// Solve Px = b -> x = P⁻¹b
VectorXd x = P.inverse() * b;

Alternatively, since the inverse of a permutation matrix is its transpose (a property of orthogonal matrices), you can also write:

VectorXd x = P.transpose() * b;

Both approaches are equally efficient—permutation operations are O(n) time (just reordering elements of b), which is negligible compared to the decomposition or triangular solves.

A quick reminder: In the full decomposition A = P⁻¹LDLᵀP, the permutation P is used to reduce fill-in during the decomposition. When you call solver.solve(b), Eigen automatically handles applying P and P⁻¹ internally, but if you need to work with P directly, the above code is the way to go.

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