重用Eigen::SimplicialLLT的符号分解
Re-use Eigen::SimplicialLLT's symbolic decomposition
我对 Eigen 库的 API 有点挣扎,即用于稀疏矩阵的 Cholesky 分解的 SimplicialLLT class。
我有三个矩阵需要分解并稍后用于求解许多方程系统(只改变右侧) - 因此我只想分解一次这些矩阵然后重新使用它们。此外,它们都具有相同的稀疏模式,所以我只想进行一次符号分解,然后将其用于所有三个矩阵的数值分解。根据文档,这正是 SimplicialLLT::analyzePattern
和 SimplicialLLT::factor
方法的用途。 但是,我似乎无法找到一种方法将所有三个因素都保留在内存中。
这是我的代码:
我的 class 中有这些成员变量 我想填充这些因素:
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 在这方面没有任何变化)。 这里的问题是,随后对 SimplicialLLT 的同一实例进行因式分解的调用将覆盖之前计算的因数,同时,我无法找到一种方法来保留它的副本。 我查看了源代码,但我不得不承认这并没有多大帮助,因为我看不到任何复制底层数据结构的简单方法。在我看来这是一个相当常见的用法,所以我觉得我遗漏了一些明显的东西,请帮忙。
我试过的:
我尝试简单地使用 choleskyB = choleskyA
希望默认的复制构造函数能够完成任务,但我发现基础 classes 被设计为不可复制的。
我可以从 choleskyA
获取 L 和 U 矩阵(它们有一个 getter),复制它们并仅存储它们,然后基本上复制粘贴 [的内容=50=]()(下面复制粘贴)直接用之前存的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 时才有意义吗?)。我希望这不是我需要做的...
最后一点 - 添加必要的 getters/setters 到 Eigen classes 并编译 "my own" Eigen 不是一个选项(好吧,不是一个好的)因为这段代码将(希望)进一步重新分发为开源,所以会很麻烦。
这是一个非常不寻常的模式。在实践中,符号分解与数值分解相比非常便宜,所以我不确定它是否值得费心。解决此模式的最干净的解决方案是让 SimplicialL?LT 可复制。
我对 Eigen 库的 API 有点挣扎,即用于稀疏矩阵的 Cholesky 分解的 SimplicialLLT class。
我有三个矩阵需要分解并稍后用于求解许多方程系统(只改变右侧) - 因此我只想分解一次这些矩阵然后重新使用它们。此外,它们都具有相同的稀疏模式,所以我只想进行一次符号分解,然后将其用于所有三个矩阵的数值分解。根据文档,这正是 SimplicialLLT::analyzePattern
和 SimplicialLLT::factor
方法的用途。 但是,我似乎无法找到一种方法将所有三个因素都保留在内存中。
这是我的代码:
我的 class 中有这些成员变量 我想填充这些因素:
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 在这方面没有任何变化)。 这里的问题是,随后对 SimplicialLLT 的同一实例进行因式分解的调用将覆盖之前计算的因数,同时,我无法找到一种方法来保留它的副本。 我查看了源代码,但我不得不承认这并没有多大帮助,因为我看不到任何复制底层数据结构的简单方法。在我看来这是一个相当常见的用法,所以我觉得我遗漏了一些明显的东西,请帮忙。
我试过的:
我尝试简单地使用 choleskyB = choleskyA
希望默认的复制构造函数能够完成任务,但我发现基础 classes 被设计为不可复制的。
我可以从 choleskyA
获取 L 和 U 矩阵(它们有一个 getter),复制它们并仅存储它们,然后基本上复制粘贴 [的内容=50=]()(下面复制粘贴)直接用之前存的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 时才有意义吗?)。我希望这不是我需要做的...
最后一点 - 添加必要的 getters/setters 到 Eigen classes 并编译 "my own" Eigen 不是一个选项(好吧,不是一个好的)因为这段代码将(希望)进一步重新分发为开源,所以会很麻烦。
这是一个非常不寻常的模式。在实践中,符号分解与数值分解相比非常便宜,所以我不确定它是否值得费心。解决此模式的最干净的解决方案是让 SimplicialL?LT 可复制。