How can I solve linear equation AX=b for sparse ma

Question

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.

Answer 1

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).

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.

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

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;

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.

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.