Решение для Lx = b и Px = b, когда A = LLt
Я разлагаю разреженную матрицу SPD A, используя Eigen. Это будет либо разложение LLt, либо LDLt (Холецкий), поэтому мы можем предположить, что матрица будет разложена как 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 эффективно и стабильно.
Решение
Посмотри как _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
потому что порядок матриц меняется на обратный при взятии обратного произведения.
Другие решения
Других решений пока нет …