Eigen LDLT Cholesky 分解就地

Question

我试图让 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？

Answer 1

一种选择是只分配一次 LDLT 求解器，然后调用计算方法：

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

如果这也不是一个选项，您可以对 ldlt 对象存储的矩阵进行 const 转换：

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

玩A，然后：

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

这很难看，但只要 MatType 是主要列，这应该可以工作。