重用Eigen :: SimplicialLLT的符号分解 [英] Re-use Eigen::SimplicialLLT's symbolic decomposition

问题描述

我在使用本征库的API（即用于稀疏矩阵的Cholesky因子分解的SimplicialLLT类）方面有些挣扎。
我需要分解三个矩阵，然后用它们来求解许多方程组（仅更改右侧）-因此，我只想分解一次这些矩阵，然后重新使用它们。而且，它们都具有相同的稀疏模式，因此我只想进行一次符号分解，然后将其用于所有三个矩阵的数值分解。根据文档，这正是 SimplicialLLT :: analyzePattern 和 SimplicialLLT :: factor 方法的用途。 但是，我似乎找不到一种将所有三个因素都保留在内存中的方法。
这是我的代码：

I am struggling a bit with the API of the Eigen Library, namely the SimplicialLLT class for Cholesky factorization of sparse matrices. I have three matrices that I need to factor and later use to solve many equation systems (changing only the right side) - therefore I would like to factor these matrices only once and then just re-use them. Moreover, they all have the same sparcity pattern, so I would like to do the symbolic decomposition only once and then use it for the numerical decomposition for all three matrices. According to the documentation, this is exactly what the SimplicialLLT::analyzePattern and SimplicialLLT::factor methods are for. However, I can't seem to find a way to keep all three factors in the memory. This is my code:

我的班级中有这些成员变量，我想补充以下因素：

I have these member variables in my class I would like to fill with the factors:

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

然后我创建三个稀疏矩阵A，B和C，并希望将它们分解为因数：

Then I create the three sparse matrices A, B and C and want to factor them:

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个向量：

And later I can use them for solutions over and over again, changing just the b vectors of right sides:

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

此方法有效（我认为），但对analyticPattern的第二次和第三次调用是多余的。我想做的事情是这样的：

This works (I think), but the second and third call to analyzePattern are redundant. What I would like to do is something like:

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

但这不是当前API的选项（我们使用Eigen 3.2.3，但如果正确地看到在3.3.2中这方面没有变化）。 这里的问题是，随后在SimplicialLLT的同一实例上进行分解的调用将覆盖先前计算出的因子，同时，我找不到一种方法来保留它的副本。我查看了源代码，但是我不得不承认这并没有太大帮助，因为我看不到任何简单的方法来复制基础数据结构。在我看来，这是一种相当常见的用法，所以我觉得我缺少明显的东西，请帮忙。

But that is not an option with the current API (we use Eigen 3.2.3, but if I see correctly there is no change in this regard in 3.3.2). The problem here is that the subsequent calls to factorize on the same instance of SimplicialLLT will overwrite the previously computed factor and at the same time, I can't find a way to make a copy of it to keep. I took a look at the sources but I have to admit that didn't help much as I can't see any simple way to copy the underlying data structures. It seems to me like a rather common usage, so I feel like I am missing something obvious, please help.

我尝试过的方法：

我尝试使用简单的 choleskyB = choleskyA 希望默认的复制构造函数能够完成任务，但是我发现

I tried using simply choleskyB = choleskyA hoping that the default copy constructor will get things done, but I have found out that the base classes are designed to be non-copyable.

我可以从获取L和U矩阵（它们有一个吸气剂）。 choleskyA ，复制它们并仅存储这些内容，然后基本上复制粘贴SimplicialCholeskyBase :: __ solve_impl（）的内容（下面复制粘贴），以编写使用先前存储的方法来解决自己的问题

I can get the L and U matrices (there's a getter for them) from choleskyA, make a copy of them and store only those and then basically copy-paste the content of SimplicialCholeskyBase::_solve_impl() (copy-pasted below) to write the method for solving myself using the previously stored L and U directly.

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;
}

...但这是一个丑陋的解决方案，我可能会搞砸因为我对过程没有很好的了解（因为我正在执行LLT，所以我不需要上面的代码中的m_diag，对吗？仅在使用LDLT时才有意义）。我希望这不是我需要做的...

...but that's quite an ugly solution plus I would probably screw it up since I don't have that good understanding of the process (I don't need the m_diag from the above code since I am doing LLT, right? that would be relevant only if I was using LDLT?). I hope this is not what I need to do...

最后一点-向Eigen类中添加必要的getter / setter并编译我自己的 Eigen是不是一种选择（当然，不是一种好的选择），因为此代码将（希望）作为开源进一步分发，因此会很麻烦。

A final note - adding the necessary getters/setters to the Eigen classes and compiling "my own" Eigen is not an option (well, not a good one) as this code will (hopefully) be further redistributed as open source, so it would be troublesome.

推荐答案

这是一个非常不寻常的模式。实际上，与数字分解相比，符号分解非常便宜，因此我不确定是否值得为此烦恼。解决此模式的最干净方法是让SimplicialL？LT可复制。

This is a quite unusual pattern. In practice the symbolic factorization is very cheap compared to the numerical factorization, so I'm not sure it's worth bothering much. The cleanest solution to address this pattern would be to let SimplicialL?LT to be copiable.

这篇关于重用Eigen :: SimplicialLLT的符号分解的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！