SimplicialLLT 返回错误的cholesky 因子

Posted 2023-02-22

【中文标题】SimplicialLLT 返回错误的cholesky 因子

【英文标题】：SimplicialLLT returns wrong cholesky factors

【发布时间】：2018-11-09 12:02:03

【问题描述】：

我想使用 Eigen 来计算稀疏矩阵的 Cholesky 分解。但是，结果不正确，我无法找到原因。如何获得正确答案？

是否在 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

【参考方案1】：

SimplicialLLT 计算 Cholesky 分解直至排列。使用Eigen::NaturalOrdering 使用身份排列，您会得到您期望的结果：

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

【参考方案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;

【讨论】：

感谢您的回答。添加NaturalOrdering 会导致编译错误error: type/value mismatch at argument 2 in template parameter list for ‘template<class _MatrixType, int _UpLo, class _Ordering> class Eigen::SimplicialLLT’。我想我还必须提供int _UpLo 才能使其正常工作。我试图在文档中找出如何做到这一点，它说the triangular part that will be used for the computations. It can be Lower or Upper. Default is Lower.，但不是如何实际添加它。你也可以帮忙吗？

请从下面的答案中查看示例代码。

对，有错别字。