简单LLT返回错误的Cholesky因子 [英] SimplicialLLT returns wrong cholesky factors
本文介绍了简单LLT返回错误的Cholesky因子的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!
问题描述
我想使用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
推荐答案
不要忘记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;
这篇关于简单LLT返回错误的Cholesky因子的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!
查看全文