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

Question

我正在使用 Eigen 分解稀疏 SPD 矩阵 A。它将是 LLt 或 LDLt 分解（Cholesky），因此我们可以假设矩阵将分解为A = P-1 LDLt P其中 P 是置换矩阵，L 是下三角形，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);

请注意，您需要乘以Pand not Pinv，因为 的倒数 A = P^-1 L D L^t P是

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

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