Решение для 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 эффективно и стабильно.


c++ eigen sparse-matrix linear-algebra matrix-decomposition
person yannick    schedule 21.03.2018    source источник

Ответы (1)


arrow_upward
0
arrow_downward

Посмотрите, как _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

потому что порядок матриц меняется на противоположный при использовании инверсии продукта.

person chtz    schedule 23.03.2018