SimplicialLLT возвращает неверные холерные факторы

Question

Я хочу использовать Eigen для вычисления холеского разложения разреженной матрицы.Однако результат неверен, и я не могу найти причину для этого.Как я могу получить правильный ответ?

И существуют ли в Eigen специальные процедуры, использующие структуру разреженных матриц для повышения производительности (например, для полосчатых матриц, как в примере ниже или треугольных матриц)?

#include <iostream>
#include <Eigen/Sparse>
#include <Eigen/Dense>

int main() 
{

    // create sparse Matrix
    int n = 5; 
    std::vector<Eigen::Triplet<double> > ijv; 
    for(int i = 0; i < n; i++)
    {
        ijv.push_back(Eigen::Triplet<double>(i,i,1));
        if(i < n-1)
        {
            ijv.push_back(Eigen::Triplet<double>(i+1,i,-0.9));
        }
    }
    Eigen::SparseMatrix<double> X(n,n);
    X.setFromTriplets(ijv.begin(), ijv.end());
    Eigen::SparseMatrix<double> XX = X * X.transpose();

    // Cholesky decomposition 
    Eigen::SimplicialLLT <Eigen::SparseMatrix<double> > cholesky;
    cholesky.analyzePattern(XX);
    cholesky.factorize(XX);

    std::cout << Eigen::MatrixXd(XX) << std::endl;
    std::cout << Eigen::MatrixXd(cholesky.matrixL()) << std::endl;

}

Матрицы выглядят следующим образом:

Вход XX:

   1 -0.9    0    0    0
-0.9 1.81 -0.9    0    0
   0 -0.9 1.81 -0.9    0
   0    0 -0.9 1.81 -0.9
   0    0    0 -0.9 1.81

Выход (cholesky.matrixL()):

  1.34536         0         0         0         0
-0.668965   1.16726         0         0         0
        0 -0.771039    1.1025         0         0
        0         0         0         1         0
        0         0 -0.816329      -0.9  0.577587

Как это должно выглядеть (X):

   1    0    0    0   0
-0.9    1    0    0   0
   0 -0.9    1    0   0
   0    0 -0.9    1   0
   0    0    0 -0.9   1

Answer 1

SimplicialLLT вычисляет разложение Холецкого с точностью до перестановки. Используя Eigen::NaturalOrdering, используется перестановка идентификаторов, и вы получите ожидаемый результат:

Eigen::SimplicialLLT <Eigen::SparseMatrix<double>, Eigen::Lower, Eigen::NaturalOrdering<int>> cholesky;
cholesky.compute(XX);

Answer 2

Не забывайте, что SimplicialLLT не разлагает A = L * L^T, но P * A * P^T = L * T^T с P матрицей перестановок. Если вам нужно P для идентификации, используйте NaturalOrdering:

Eigen::SimplicialLLT<Eigen::SparseMatrix<double>, Eigen::Lower, Eigen::NaturalOrdering<int> > cholesky;