Решить для обратного квадратного корня

У меня есть положительно определенная матрица A, из которой я уже рассчитал разложение Холецкого: A = LDL ^ T.
Для некоторого вектора x я бы хотел вычислить S ^ {- 1} x, где S — квадратный корень из A. На данный момент я делаю

Eigen::SelfadjointEigenSolver<Eigen::MatrixXd> es(A);
Eigen::MatrixXd Si(es.operatorInverseSqrt());
return Si*get_x();

Это стабильный способ сделать это вычисление? Я думал, что вычисление инверсий было плохой вещью вообще. Есть ли способ использовать уже выполненную декомпозицию LDLT? Я чувствую, что это возможно, потому что это то, что на самом деле происходит за кулисами LDLT::solve()!

2

Решение

Вот полный код для решения задачи для симметричной матрицы A и общей правой части b (вектор или матрица). Я не смог найти в Интернете ничего, с чем можно было бы поиграть (или просто скопировать-вставить), поэтому я написал одно.

Метод stable_cholesky_solver выполняет работу по поиску квадратного корня с использованием стабильной декомпозиции lldt() который использует поворот. main проверяет, что он делает все, что должен, и представляет способ достижения той же цели, используя менее стабильный (но более быстрый) llt() разложение.
Смотрите первые несколько строк документация чтобы понять мои обозначения L, P, D.

#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;

Matrix<double, Dynamic, Dynamic> stable_cholesky_solver(
LDLT<MatrixXd> ldltDecomp,
Matrix<double, Dynamic, Dynamic> A,
bool transpose = false )
{

// Preparations:

// For some reason if I sub it below I get error
Matrix<double, Dynamic, Dynamic> L = ldltDecomp.matrixL();

// Number of rows is all that matters, regardless if rhs is a
// matrix or a vector
int k = A.rows();

// Manually inverting D. This procedure has the advantage that
// D^{-1/2} can also be applied to matrices.
VectorXd diag;
diag.resize(k);
for( int i = 0 ; i < k ; ++i )
diag(i) = 1. / sqrt( ldltDecomp.vectorD()(i) ) ; // Manual inversion
DiagonalMatrix<double, Dynamic > sqrtInvD = diag.asDiagonal();

// The permutation "matrix" P
Transpositions<Dynamic> P = ldltDecomp.transpositionsP();

// Holds all the computations
Matrix<double, Dynamic, Dynamic> x;

// Now the actual computation

if( !transpose ) {

// x = PA
x = P * A;

// x = L^{-1}PA
x = L.triangularView<Lower>().solve<OnTheLeft>(x);

// x = D^{-1/2}L^{-1}PA
x = sqrtInvD * x;

} else {

// x = D^{-1/2}A
x = sqrtInvD * A;

// x = L^{-t}D^{-1/2}A
x = L.triangularView<Lower>().transpose().solve<OnTheLeft>(x);

// x = P^tL^{-t}D^{-1/2}A
x = P.transpose() * x;
}
return x;

}int main()
{

int k = 3; // Dimensionality

// Define, declare and enter the problem's data
MatrixXd A;
A.resize(k, k);
MatrixXd b;
b.resize(k, 2 );

A <<
13, 5, 7 ,
5 , 9, 3 ,
7 , 3, 11;
b <<
3, 3, 4,
1,-2, 9;

cout << "Here is the " << A.rows() << " by " << A.cols() << " matrix A:\n" << A << endl;
cout << "Here is the " << b.rows() << " by " << b.cols() << " matrix b:\n" << b << endl;
cout << "Let's solve Ax = b using different methods.\n" <<endl;

// Two placeholders that will be used throughout
MatrixXd L;
MatrixXd x;

// ldlt()
cout << "\n\nUsing the stable Cholesky decompostion ldlt()" << endl;

// The object that actually holds the entire decomposition
LDLT<MatrixXd> ldltDecomp = A.ldlt();

// Direct solution, using Eigen's routines
x = ldltDecomp.solve(b);
cout << "Direct x =\n" << x << endl;
cout << "Direct b =\n" << A*x << endl;

// Manual solution - implementing the process Eigen is taking, step
// by step (in the function defined above).
x = stable_cholesky_solver( ldltDecomp, b );
x = stable_cholesky_solver( ldltDecomp, x, true );

cout << "Manual x =\n" << x << endl;
cout << "Manual b =\n" << A*x << endl;// llt()
cout << "\n\nUsing the less stable, but faster Cholesky decomposition " << "without pivoting llt()" << endl;

// Here A.llt() is the object that actually holds the decomposition
// (like ldltDecomp before) but we only need the matrix L.
L = A.llt().matrixL();

x = L.triangularView<Lower>().solve<OnTheLeft>(b);
x = L.triangularView<Lower>().transpose().solve<OnTheLeft>(x);
cout << "Manual x =\n" << x << endl;
cout << "Manual b =\n" << A*x << endl;

}
0

Другие решения

Других решений пока нет …

По вопросам рекламы [email protected]