c++ - 当 A=LLt 时求解 Lx=b 和 Px=b

Question

我正在使用 Eigen 分解稀疏 SPD 矩阵 A。它将是 LLt 或 LDLt 分解 (Cholesky)，因此我们可以假设矩阵将分解为D 对角线(可能是恒等式)。如果我这样做

SolverClassName<SparseMatrix<double> > solver;
solver.compute(A);

要解决 Lx=b 那么执行以下操作是否有效？

solver.matrixL().TriangularView<Lower>().solve(b)

类似地，要解决 Px=b 那么执行以下操作是否有效？

solver.permutationPinv()*b

我想这样做是为了高效稳定地计算bt A-1 b。

Answer 1

看看 _solve_impl 是如何为 SimplicialCholesky 实现的。本质上，你可以简单地写:

Eigen::VectorXd x = solver.permutationP()*b; // P not Pinv!
solver.matrixL().solveInPlace(x); // matrixL is already a triangularView

// depending on LLt or LDLt use either:
double res_llt = x.squaredNorm();
double res_ldlt = x.dot(solver.vectorD().asDiagonal().inverse()*x);

请注意，您需要乘以 P 而不是 Pinv，因为A = P^-1 L D L^t P 是

P^-1 L^-t D^-1 L^-1 P

因为在对乘积求逆时矩阵的顺序会颠倒。