特征LDLT Cholesky分解到位

Eigen LDLT Cholesky decomposition in-place

本文关键字：分解 Cholesky LDLT 特征      更新时间：2023-10-16

我正在尝试使用Eigen3来求解具有局部Cholesky分解的线性系统A * X = B。我负担不起在堆栈上推送任何A大小的临时文件，但在此过程中我可以自由销毁A。

不幸的是，

A.llt().solveInPlace(B);

因为A.llt()在堆栈上隐式地推送A大小的临时矩阵。对于LLT的情况，我可以访问必要的功能，比如：

// solve A * X = B in-place for positive-definite A
template <typename AType, typename BType>
void AllInPlaceSolve(AType& A, BType& B)
{
    typedef Eigen::internal::LLT_Traits<AType, Eigen::Upper> TraitsType;
    TraitsType::inplace_decomposition(A);
    TraitsType::getL(A).solveInPlace(B);
    TraitsType::getU(A).solveInPlace(B);
}

这很好，但我担心：

  • 我的矩阵A可能只是半正定的，在这种情况下需要LDLT分解
  • LLT分解为系统的解计算不必要的sqrt()

我找不到一种方法来钩住Eigen的LDLT功能，类似于上面的代码，因为代码的结构非常不同。

所以我的问题是：有没有一种方法可以使用Eigen3来求解线性系统，使用LDLT分解，使用的划痕空间不比对角矩阵D多？

一个选项是只分配一次LDLT解算器，并调用计算方法：

LDLT<MatType> ldlt(size);
// ...
ldlt.compute(A);
x = ldlt.solve(b);

如果这也不是一个选项，您可以constcast ldlt对象存储的矩阵：

LDLT<MatType> ldlt(MatType::Identity(size,size));
MatType& A = const_cast<MatType&>(ldlt.matrixLDLT());

播放A，然后：

ldlt.compute(A);
x = ldlt.solve(b);

这很难看，但只要MatType是列主，这就应该有效。

相关文章：