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

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?

• Seems to me like Matlab shouldn't even work with a matrix that large, but I'm not an experienced user. Are you sure the inverse is actually working? – EMP Sep 1 '19 at 9:01
• yes,I am sure, I implement this experiments in software "ifiss",the website is:IFISS Software personalpages.manchester.ac.uk/staff/david.silvester/ifiss/… the finite element method and fast iterative methods in incompressible fluid dynamic. A very efficient software for iterative methods. – sunshine Sep 1 '19 at 9:32
• I dont know if I'd use Matlab and fast and efficient in the same sentence but leaving that aside, I think we need more information. What are the various grid inputs, what do these mean? Have you tested the iterative solver? Are you sure it works? And I still don't believe that matlabs direct solve can actually handle a matrix that large that quickly. How do you know it properly inverted it? Did you test it? – EMP Sep 1 '19 at 12:13
• @EMP side note: there is no matrix inversion happening during A\f in Matlab. – Anton Menshov Sep 1 '19 at 15:13
• Anton thanks for the info. I still have trouble believing that a matrix that large is bring solved that fast in Matlab unless it has some verrrrry special properties. – EMP Sep 1 '19 at 19:07

There are several things to consider in this experiment:

### Why Matlab sparse direct might be "so fast":

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

• Great answer. I also think a key to OP's numerical experiments is that he seems to be actually forming the matrix $A$ and computing the matrix-vector product directly. For a Poisson problem, matvecs can be done very quickly, and that is also ignoring other direct techniques like Fast Poisson Solvers that may work depending on his domain. – whpowell96 Sep 1 '19 at 15:45
• Excellent answer,thanks for all.@whpowell96,@Anton, the reason is the special sparsity pattern in 2D poisson, and I have tried in 3D, and matlab A\b failed out of memory , so iterative method is necessary. – sunshine Sep 2 '19 at 8:48
• From my personal experience I found out that the Poisson problem is a very easy problem to solve for iterative and direct solvers. When moving on to more harder problems like structural problems, the comparison between direct and iterative changes. – vydesaster Sep 4 '19 at 17:32
• @vydesaster yeah, poisson matrix has very special structure so that matlab A\b can solve it fastly in 2 dimension, but failed in 3D for very large system size. – sunshine Sep 6 '19 at 0:40

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