如何解决稀疏矩阵的线性方程AX = b [英] How can I solve linear equation AX=b for sparse matrix

问题描述

我有一个稀疏矩阵A（120,000 * 120,000）和一个向量b（120,000），我想使用本征库求解线性系统AX = b。我尝试按照文档进行操作，但始终出现错误。我还尝试将矩阵更改为稠密并求解系统

I have sparse matrix A(120,000*120,000) and a vector b(120,000) and I would like to solve the linear system AX=b using Eigen library. I tried to follow the documentation but I always have an error. I also tried change the matrix to dense and solve the system

    Eigen::MatrixXd H(N,N);
    HH =Eigen:: MatrixXd(A);
    Eigen::ColPivHouseholderQR<Eigen::MatrixXd> dec(H);
    Eigen::VectorXd RR = dec.solve(b);  

但是我遇到了内存错误。
请帮助我找到解决此问题的方法。

but I got a memory error. Please help me find the way to solve this problem.

推荐答案

对于稀疏系统，我们通常使用迭代方法。原因之一是LU，QR ...因式分解等直接求解器会填充您的初始矩阵（从某种意义上说，最初是零分量被非零分量替换）。在大多数情况下，生成的密集矩阵不再适合内存（-> 您的内存错误）。简而言之，迭代求解器不会更改您的稀疏模式，因为它们仅涉及矩阵矢量乘积（无填充）。

For sparse systems, we generally use iterative methods. One reason is that direct solvers like LU, QR... factorizations fill your initial matrix (fill in the sense that initially zero components are replaced with nonzero ones). In most of the cases the resulting dense matrix can not fit into the memory anymore (-> your memory error). In short, iterative solvers do not change your sparsity pattern as they only involve matrix-vector products (no filling).

也就是说，您必须首先知道您的系统是否是在这种情况下，您可以使用共轭梯度法。否则，您必须使用非对称系统的方法，例如 BiCGSTAB 或 GMRES 。

That said, you must first know if your system is symmetric positive definite (aka SPD) in such case you can use the conjugate gradient method. Otherwise you must use methods for unsymmetric systems like BiCGSTAB or GMRES.

您必须知道大多数时候我们使用<一个href = https://en.wikipedia.org/wiki/Preconditioner rel = nofollow noreferrer>前置条件，尤其是在您的系统状况不佳的情况下。

You must know that most of the time we use a preconditioner, especially if your system is ill-conditioned.

看我发现的Eigen文档（不带前置条件afaik的示例）：

Looking at Eigen doc I found(example without preconditioner afaik):

  int n = 10000;
  VectorXd x(n), b(n);
  SparseMatrix<double> A(n,n);
  /* ... fill A and b ... */ 
  BiCGSTAB<SparseMatrix<double> > solver;
  solver.compute(A);
  x = solver.solve(b);
  std::cout << "#iterations:     " << solver.iterations() << std::endl;
  std::cout << "estimated error: " << solver.error()      << std::endl;

这可能是一个不错的开始（如果矩阵是SPD，请使用共轭梯度法）。请注意，这里没有前置条件，因此收敛肯定会很慢。

which is maybe a good start (use Conjugate Gradient Method if your matrix is SPD). Note that here, there is no preconditioner, hence convergence will be certainly rather slow.

回顾：大矩阵->迭代方法+前置条件。

Recap: big matrix -> iterative method + preconditioner.

这实际上是第一个基本/幼稚的解释，您可以在萨阿德的书：稀疏线性系统的迭代方法。同样，这个主题非常重要，您可以找到很多与此主题相关的其他书籍。

This is really a first "basic/naive" explanation, you can found futher information about the theory in Saad's book: Iterative Methods for Sparse Linear Systems. Again this topic is huge you can found plenty of other books on this subject.

我不是Eigen用户，但是有此处的预处理器。
->不完整的LU（不对称系统）或不完整的Cholesky（SPD系统）通常是通用的预处理器，有待测试。

I am not an Eigen user, however there are precondtioners here. -> Incomplete LU (unsymmetric systems) or Incomplete Cholesky (SPD systems) are generally good all purpose preconditioners, there are the firsts to test.

这篇关于如何解决稀疏矩阵的线性方程AX = b的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！