# Direct vs Iterative solvers choice

Is there any other reason except “the big size of matrix” that makes me prefer the use of iterative solvers than direct ones, for (linear algebraic systems)? Thanks

• Solvers for what? Jun 28, 2016 at 15:42
• linear algebraic systems @Emily
– MBM
Jun 28, 2016 at 15:49
• Related scicomp.SE question. Jun 28, 2016 at 15:51

The "big matrix" problem is certainly one reason.

Another is that iterative solvers can be run until the result is "good enough", so you are checking the results every k-th iteration.

Another reason is that the explicit matrix may not be available, so you can get Ax for any specified x, but not A itself.

Big picture:

• Gaussian elimination with partial pivoting in theory has some severe stability problems. It is technically backward stable in a fixed dimension $n$, but the constant multiplying factor is something like $2^n$. Thus in dimension greater than about $50$, machine epsilon size errors in theory could scale up to order 1 errors (roughly speaking "no correct digits" in the solution). There is an explicit example of this given in Trefethen and Bau (it is something like the identity with the last column instead filled with 1s). Surprisingly, this phenomenon is quite rare in applications; I have yet to see a good explanation for why.
• You may not be aware that there are two different aspects to the "the matrix is too big" problem. One is "I can't assemble the matrix in full form". Think about it: suppose you have 8 GB of RAM, we know a double precision number is 8 bytes, so you can hold about a billion double precision numbers in RAM at a time. Thus the largest matrices you can store in full form with fast access are about 50000 by 50000. If you need bigger, then $A$ is only available to you as "a black box that computes $Ax$ given $x$", so iterative methods are all you can do.
• The other aspect of the "the matrix is too big" problem is that you might not have $O(n^3)$ time to burn. Taking those 50000 by 50000 matrices for illustration, running Gaussian elimination on them takes about $10^{14}$ operations. The number of floating point operations we can do in a second varies by machine and how nicely the problem parallelizes, but if you take it to be 10 billion, Gaussian elimination on a 50000 by 50000 matrix would take something like 2 hours. Faster than that would be good.
• Even when you can assemble $A$, sometimes it naturally comes to you as a black box that computes $Ax$ given $x$, and assembling it as a matrix may cost you additional time. This happened to me: I had a linear function acting on Fourier coefficients which was defined in real space, so just to compute the matrix for it, I had to do $2n$ $n$-dimensional FFTs. In my case $n$ was just $60$ so this wasn't so bad, but it would've been much worse with large $n$.
• Least squares problems are often very badly conditioned when approached with Gaussian elimination on the normal equations. They can be so badly conditioned that the only good options among "direct" methods are based on QR decomposition (algebraically reduce $(A^T A)^{-1} A^T b$ to $R^{-1} Q^T b$) and SVD (algebraically reduce to $V \Sigma^{-1} U^T b$). And even that SVD option is not really direct, because there is no direct method for the SVD itself.

The small picture boils down to "what properties of $A$ can we exploit to fix problems like the ones above?" For example, if $A$ is quite sparse, $Ax$ may be quite fast to compute, so even $n$ matrix multiplications could turn out to be much faster than $n^3$ time (even though the worst case for $n$ matrix multiplications is $n^3$ time).

• With respect to point 1: see this paper, among a number of other references. Jun 29, 2016 at 2:12