How large is large for direct solvers?

Question

Let us say I want to solve a large sparse linear system. It is said that iterative solvers should be better than direct solvers in this case. But how large is large? What is the exact threshold beyond which I must use iterative solvers? A thousand by thousand matrix? A million by million matrix? I understand this might depend on the particular situation but I want to have some idea of how large number we are talking about. Currently I have absolutely no idea what kind of size they are talking about when textbooks say "iterative methods are for large systems."

It is similar to having no idea of what temperatures should be regarded as hot and what temperatures should be regarded as cold. I want somebody tell me "You can quite safely assume that beyond 40ºC is hot in normal circumstances, even though people whose primary occupation is to burn things might consider even 100ºC as cold."

Let me make the question concrete. Please post the result of specific tests performed on direct and/or iterative methods. Please include information on the PDE (or other) problem, preconditioners, sizes, dimensions, etc, anything that might be important.

[Here is a frustrated rant: It feels like all textbooks and even some comments and answers here are trying to hide the exact numbers, similarly to some videos on the internet saying they have something important to say but never really say it, and in the end ask you to pay if you want the answer. They say "it depends on this depends on that" but never say any concrete numbers to at least give some rough idea of the scale. I am sure there is a good reason, but I suspect it is because the experts don't remember how they were when they first learned the subject. It is almost as if they only want to talk amongst themselves, and only care what the other experts think of what they say even when they are talking to newcomers.]

Answer 1

The answer here mostly depends on how good the preconditioner you have for the iterative solver. If you don't have a good preconditioner, direct methods tend to be the best until you run out of ram. However, for pdes and similar problems where you are solving a problem with discretization, geometric multi-grid methods where you precondition your big linear solve with an upscaled solution to a lower resolution solve can be extremely effective for medium sized data.

no matter which method you use, make sure you are linking against good blas libraries and have all your threading set up right, because if you mess that up, you can suffer by much more than the difference between a direct and iterative solver.

Answer 2

Here is a "real live" test that I did: A few weeks ago I did a simulation of a large grid of 1Ohm resistors. Both with an iterative relaxation solver. ( The code is gitlab: grid-of-1ohm-resistors The simulations is with python/numpy/scipy:

With the direct/sparse (scipy) approach I could simulate at grid of 2000x2000 elements (actually a quarter of the grid with that size since the problem is symmetric so it would be a real 4000x4000 grid).

The simulation (not optimized only used 1 core and also has some room for improvement) took about 390seconds on my AMD Ryzen 5 5500U 1.4Ghz CPU). And used up most of the 16GB RAM. (a 3000x3000 grid would fail). The result was with 14 digits precision. (The matrix of the 2000x2000 grid would be a 4000000x400000 matrix but extremely sparse so only 5 entries in each row so only about 1 in a million elements non-zero).

Compare this to the iterative: There you less RAM so you can tackle larger simulations and you can control the precision by deciding how much iterations you run. For lower precision you will get faster results but if you want the same precision you would need much longer computations. For a runtime of about 300 seconds you only get a precision of 4 digits. (5000 iterations). (Also not optimized, only using one core) If you want more digits you need more iterations.

So iterative relaxation method: good if you only need rough estimates and do not care for precision or if you have do much larger simulations.

	direct solution with sparse matrix	iterative solution with relaxation method
size on a 16GB RAM machine (room for optimization)	ca 4M	ca 100M
speed	fast	slow except if you trade precision for speed
precision	always high	dependent on number of iterations
easy to parallelize	no (unless your library for sparse matrices supports it)	yes

conclusion: direct calculation with sparse matrices is preferable in most cases.

Answer 3

This is mostly a response to your last edit:

"It feels like all textbooks and even some comments and answers here are trying to hide the exact numbers, similarly to some videos on the internet saying they have something important to say but never really say it, and in the end ask you to pay if you want the answer. They say "it depends on this depends on that" but never say any concrete numbers to at least give some rough idea of the scale. I am sure there is a good reason, but I suspect it is because the experts don't remember how they were when they first learned the subject. It is almost as if they only want to talk amongst themselves, and only care what the other experts think of what they say even when they are talking to newcomers."

A lot of knowledge we generate as experts are free to access, you can read both of my theses which are on preconditioners and iterative methods and some papers I wrote have preprints posted on arXiv. Here Prof. Wolfgang Bangerth is a big expert on numerical PDEs, iterative methods, and he is one of lead investigators of deal.ii (open source scientific computing software).

However, the question you ask is hard to answer succintly. For example, knl, Julian Roth and Prof. Bangerth mentioned it depends on the sparsity pattern, preconditioner and application domain. Which are all correct. On the other hand, some problems do not have good preconditioners (for example, Laplacian solved with multigrid method will always outperform direct solvers due to pure computational complexity advantages). Also many linear systems (especially coming from PDE domains) get more ill-conditioned as the system size grows, so if there is no preconditioner available to make the system well-conditioned then the iterative methods will not work. Then there are scalability and memory considerations, out-of-core computations, required level of precision and/or accuracy (precision and accuracy are different terms, they affect the decision differently)

Short answer, experts do not know the answer; we are searching for it. My recent survey was about problems from a large set of domains, PDEs, optimization, graph theory, ... What I found was mixed, but to a decent degree of certainty (as physicist say 1-$\sigma$), I can say that the most important thing is the target hardware and density of factors after factorization. I have access to decent hardware and factors fit into the memory so I found that over the testset of problems (obtained from SuiteSparse Matrix Collection) direct methods outperform the iterative methods. Exceptions being tridiagonal problems of certain domains, applying algebraic multigrid for some (nice) PDE problems, or some circuit design problems which lend themselves to iterative methods without preconditioning.