Why iterative method: AMG preconditioned PCG is slower than Matlab direct method 'A\b'?

Question

Recently, I have met a question that

a saying goes that for large linear system: iterative methods are required because of memory problem of direct methods.

But when I implement some experiments in Matlab, I find a strange phenomenon: I used the Poisson equation in 2D and use Q1 element, grid input is 10, and I got a system: $$Ax = f$$ where $A$ is $1050625\times 1050625$, which is large and sparse.

In principle, we should use iterative methods, such as PCG or minres built-in Matlab, with AMG preconditioner. But when I input A\f in the command window, Matlab direct method only costs $4.588234$ seconds, which is fast.

Then, I want to test PCG with AMG preconditioner. I find the time for setup of the AMG preconditioner to be very very long. First, I thought this is because the size of system is not large enough, then I use the grid input =11, and the memory breaks down, the system matrix $A$ cannot be generated, so in my PC, I can not get the result that iterative method is better than direct method for large sparse system. Why? Is it the reason that the size is not enough? but the larger system cannot be generated in our personal computer.

I also believe that for large sparse system, iterative methods are necessary, but the numerical results are contrary: Matlab's A\b is so fast.

How should we understand the saying "iterative methods are better than direct methods for large sparse system"? Could you please give me some suggestions?

Answer 1

There are several things to consider in this experiment:

Why Matlab sparse direct might be "so fast":

(for your particular test)

• In 2D (of course, problem-dependent), your matrix $A$ arising after FEM discretization, after some reorderings might appear to be "close to banded" structure. The smaller is the bandwidth of $A$, the more efficient a sparse direct method could handle it. And if the bandwidth is small enough, Matlab probably would choose a specialized banded solver.
• I am not aware of what particular implementation of algebraic multigrid (AMG) preconditioner you are using. While there might be internal performance problems inside it, the AMG itself might be an overkill for your problem.
• sparse linear solvers are now a thing. There was a lot of progress in the last years, and at least part of the effect you are seeing should be attributed to it.

Why one might want to use an iterative solver, no matter what:

• In general, even the best sparse direct solver during factorization would still generate the fill-in. Thus, the matrix after the factorization would take [significantly] more space compared to the original one, which is certainly an issue for large problems.
• Acceleration of iterative solvers is a more mature field. In different application areas, there is a multitude of fast algorithms (based on, say, Fast Fourier Transform – FFT, Fast Multipole Method – FMM, and others) which would accelerate the matrix-vector product, naturally fitting the iterative linear solver route.

Things to try:

• Try an iterative solver with a vanilla diagonal preconditioner on your problem. See, what is the number of iterations and decide whether a better preconditioner (like AMG) is even worth considering.
• Consider a test in 3-D if this is in the interest of your application area.
• Check the bandwidth and sparsity pattern of your matrix $A$ arising from 2-D Poisson. It might happen, that you are solving a very special case, not a general sparse matrix.

---

How to understand the "saying":

for large linear system: iterative methods are required because of memory problem of direct methods.

You may want to critically look at it. Pretty much as at any general blanket statement.

While iterative methods have their advantages, certain problems would still call for direct methods. Moreover, there are plenty of fast direct methods, where the factorization would be also accelerated (say, Hierarchical matrices applied to FEM). The field of unaccelerated sparse linear solvers became much more promising in the last 15 years, I would say. So, that saying might have been a dogma 20 years ago, now it is at least much weaker.

Answer 2

Thanks for all your attention. below is the reply from a professor:

The MATLAB sparse solver is a very efficient way of solving linear systems associated with the two-dimensional Laplacian operator. One reason for this is that the CHOLMOD solver is very effectively multithreaded so it can use all available processors in the solution process. For example, my Apple laptop is an I9 six-core architecture and I can see that all six are fully used when I solve the problem you discuss below. In contrast the AMG grid setup is interpreted code and, as you have observed, is extremely slow in a MATLAB environment. It is, however, memory efficient.

I have tried the numerical experiments in 3D, using 5 points difference to discretize the poisson equation:

$$\left\{\begin{array}{l}{-\Delta u=f}, \quad {(x, y,z) \in G=(-1,1)^3} , \\ {u=g,\quad (x, y,z) \in \partial G.}\end{array}\right.$$

when the system size become 1,000,000 X 1,000,000, the matlab command A\b is out of memory. the Matlab code is as follows:

%%  poisson in 2D and 3D 5 points difference matrix 
clc;clear;
n=10;
e=ones(n,1);
B = [-1 2 -1].*e;
d = [-1 0 1];
Tn = spdiags(B,d,n,n);
e=ones(n-1,1);
I = speye(n);
%  2D
Tn_I = kron(Tn,I);
I_Tn = kron(I,Tn);
A = Tn_I+I_Tn;
%  3D
A = kron(Tn_I,I)+kron(I,Tn_I)+kron(I,I_Tn);
b = sum(A,2);
tic;
A\b;toc