使用 Eigen 库执行稀疏LU 并显示 L&U?

Use Eigen library to perform sparseLU and display L & U?

本文关键字:显示 LU Eigen 执行 使用      更新时间:2023-10-16

我是特征的新手,我正在处理稀疏LU问题。我发现如果我创建一个向量b(n) Eigen可以计算Ax=b方程的x(n)

问题:

  1. 如何显示L &U是原始矩阵A的因式分解结果?

  2. 如何在特征中插入非零?现在我只是用一个小的稀疏矩阵来测试所以我一个接一个地插入非零矩阵,但是如果我有一个大的矩阵,我怎么把矩阵输入到我的程序中呢?

我意识到这个问题很久以前就被问过了。显然,参考Eigen文档:

矩阵L的表达式,内部存储为超节点。这个表达式唯一可用的操作是三角解

所以没有办法将它转换成一个实际的稀疏矩阵来显示它。Eigen::FullPivLU执行密集分解,在这里对我们没有用处。在一个大的稀疏矩阵上使用它,当我们试图将其转换为密集矩阵时,我们将很快耗尽内存,并且计算分解所需的时间将增加几个数量级。

另一种解决方案是使用来自Suite Sparse的CSparse库:

extern "C" { // we are in C++ now, since you are using Eigen
#include <csparse/cs.h>
}
const cs *p_matrix = ...; // perhaps possible to use Eigen::internal::viewAsCholmod()
css *p_symbolic_decomposition;
csn *p_factor;
p_symbolic_decomposition = cs_sqr(2, p_matrix, 0); // 1 = ordering A + AT, 2 = ATA
p_factor = cs_lu(p_matrix, m_p_symbolic_decomposition, 1.0); // tol = 1.0 for ATA ordering, or use A + AT with a small tol if the matrix has amostly symmetric nonzero pattern and large enough entries on its diagonal
// calculate ordering, symbolic decomposition and numerical decomposition
cs *L = p_factor->L, *U = p_factor->U;
// there they are (perhaps can use Eigen::internal::viewAsEigen())
cs_sfree(p_symbolic_decomposition); cs_nfree(p_factor);
// clean up (deletes the L and U matrices)

请注意,尽管它不像某些特征函数那样使用显式向量化,但它仍然相当快。CSparse也非常紧凑,它只有一个头文件和大约30个没有外部依赖的.c文件。易于合并到任何c++项目中。实际上不需要包含所有的Suite Sparse

如果您将使用Eigen::FullPivLU::matrixLU()到原始矩阵,您将收到LU分解矩阵。如果需要分别显示L和U,可以使用triangularView<模式>方法。在Eigen wiki中你可以找到一个很好的例子。在矩阵中插入非零取决于数字,而不是你想要的数字。Eigen具有方便的语法,因此您可以轻松地在循环中插入值:

for(int i=0;i<size;i++)
  {
  for(int j=size;j>someNumber;j--)
    {
    matrix(i,j)=yourClass.getNextNumber();
    }
  }