重用特征：：辛普利夏LLT的符号分解

Question

我在使用特征库的API时有点困难，即用于稀疏矩阵Cholesky因子分解的SimplicialLLT类。我有三个矩阵，我需要因子，然后用来求解许多方程组(只改变右侧)——因此，我只想因子这些矩阵一次，然后重复使用它们。此外，它们都有相同的稀疏模式，所以我只想做一次符号分解，然后用它来对所有三个矩阵进行数值分解。根据文档，这正是SimplicialLLT::analyzePattern和SimplicialLLT::factor方法的用途然而，我似乎找不到把这三个因素都记在记忆里的方法这是我的代码：

我在班上有这些成员变量，我想用因子来填充：

Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> choleskyA;
Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> choleskyB;
Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> choleskyC;

然后我创建了三个稀疏矩阵A、B和C，并想对它们进行因子化：

choleskyA.analyzePattern(A);
choleskyA.factorize(A);
choleskyB.analyzePattern(B); // this has already been done!
choleskyB.factorize(B);
choleskyC.analyzePattern(C); // this has already been done!
choleskyC.factorize(C);

后来我可以一次又一次地用它们来求解，只改变右边的b向量：

xA = choleskyA.solve(bA);
xB = choleskyB.solve(bB);
xC = choleskyC.solve(bC);

这是有效的(我认为)，但第二次和第三次调用analyzePattern是多余的。我想做的是：

choleskyA.analyzePattern(A);
choleskyA.factorize(A);
choleskyB = choleskyA.factorize(B);
choleskyC = choleskyA.factorize(C);

但对于当前的API来说，这不是一个选项(我们使用Eigen 3.2.3，但如果我看得正确的话，在3.3.2中没有这方面的变化)。这里的问题是，对SimpliciaLLT的同一实例进行因子分解的后续调用将覆盖之前计算的因子，同时，我找不到方法来复制它以保留我看了一下源代码，但我不得不承认，这并没有多大帮助，因为我看不到任何简单的方法来复制底层数据结构。在我看来，这是一个相当常见的用法，所以我觉得我错过了一些明显的东西，请帮忙。

我尝试过的：

我试着简单地使用choleskyB = choleskyA，希望默认的复制构造函数能完成任务，但我发现基类被设计成不可复制的。

我可以从choleskyA中获得L和U矩阵(它们有一个getter)，复制它们并只存储它们，然后基本上复制粘贴SimplicialCholeskyBase:：_solve_impl()的内容(复制粘贴在下面)，直接使用之前存储的L和U编写解算方法。

template<typename Rhs,typename Dest>
void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
{
eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
eigen_assert(m_matrix.rows()==b.rows());
if(m_info!=Success)
return;
if(m_P.size()>0)
dest = m_P * b;
else
dest = b;
if(m_matrix.nonZeros()>0) // otherwise L==I
derived().matrixL().solveInPlace(dest);
if(m_diag.size()>0)
dest = m_diag.asDiagonal().inverse() * dest;
if (m_matrix.nonZeros()>0) // otherwise U==I
derived().matrixU().solveInPlace(dest);
if(m_P.size()>0)
dest = m_Pinv * dest;
}

但这是一个相当丑陋的解决方案，而且我可能会把它搞砸，因为我对这个过程没有很好的理解(我不需要上面代码中的m_diag，因为我在做LLT，对吧？只有当我使用LDLT时，这才有意义？)。我希望这不是我需要做的…

最后要注意的是，向Eigen类添加必要的getter/setter并编译"我自己的"Eigen不是一个选项(好吧，不是一个好的选项)，因为这段代码将(希望)作为开源进一步重新分发，所以会很麻烦。

Answer 1

这是一个非常不寻常的模式。在实践中，与数字因子分解相比，符号因子分解非常便宜，所以我不确定它是否值得花太多精力。解决这种模式的最干净的解决方案是让SimplicialL？LT是可复制的。