# How does gmres method iteration behave for this non-diagonalizable matrix?

Recently, I have been studied my lessons about gmres iteration, probably the most popular iteration method for general large sparse linear system of equations Ax=b. And the convergence is obtained under the condition of which the matrix A is diagonalizable, but in real applications, this is seldom the case where such a good diagonalizable matrix A. So, is there any other convergence results for GMRES iteration? or in other words, when we choose a preconditioner M for solving $$M^{-1}Ax=M^{-1}b,$$ which aspects should we consider to guarantee the better convergence for this preconditioned GMRES? Because $$M^{-1}A$$ can not be diagonalizable in general case, is it enough to cluster the spectral distributions of the $$M^{-1}A$$? What aspects of the preconditioned matrix $$M^{-1}A$$ should we discuss?

A strange example is as follows: $$A=\left(\begin{array}{cccccc}{1} & {c_{1}} & {} & {} & {} & {} \\ {0} & {1} & {c_{2}} & {} & {} & {} \\ {} & {} & {\ddots} & {\ddots} & {} & {} \\ {} & {} & {} & {1} & {c_{n-2}} & {} \\ {} & {} & {} & {} & {1} & {c_{n-1}} \\ {} & {} & {} & {} & {} & {1}\end{array}\right)\in\mathbb{R}^{n\times n},$$ with $$c_i\neq 0$$, we set $$b=A*[1, 1,...,1]^T$$ to solve this system $$Ax=b$$ using matlab command gmres. And we compare this result with matlab A\b, but gmres fails (my CPU memory is 8GB), A\b instead succeed. I just cannot understand why this happened, because we often say that for large sparse matrix, direct method can be beaten by iterative method. Can you give me some suggestions about gmres convergence and some explanations about the above example Ax=b, because all the eigenvalues of matrix A is unit, theoretically, gmres will converges very fast. Thanks very much.

Below is my simple test matlab code

%   just a simple test
clc;clear;
rng(0)% fix the random number
n=1e+5;
c = rand(n,1);% c_i
A = spdiags([ones(n,1)  c],[0 1],n,n);
b = A*ones(n,1);
%   mdirect method
A\b;
%   iterative method
x=gmres(A,b,[],[],n);


• Well for the direct method the LU factorisation of this matrix is fairly easy to work out ... Maybe MATLAB can detect this and take advantage to reduce the problem to at worst O(N**2) - and better if it can take advantage of the sparsity in the triangular matrix. Basically his is an extremely easy matrix for classical direct methods, an O(N) solution method is extremely easy to write. Nov 4 '19 at 8:52
• Thanks Professor. When we use a preconditioner, in general, if we can cluster the eigenvalues then, a better convergence can be obtained, right? It seems that in general case, only eigenvalues can be discussed about a preconditioner for $M^{-1}A$. Thanks professor. Nov 4 '19 at 9:19
• Firstly I am a mere Doctor. Secondly iterative methods are not my area at all, so I'll leave the GMRES side to somebody who really knows what they are saying. Nov 4 '19 at 9:52

Unfortunately, convergence of GMRES does not have a clear dependence on the distribution of eigenvalues.

It was proved by Greenbaum, Ptak and Strakos in 1996 that you can construct examples with an arbitrary spectrum and an arbitrary convergence history: that is, give me any $$n$$ nonzero complex numbers, and any decreasing sequence $$\|r_k\|$$, and I can construct a matrix $$A$$ with that spectrum and a vector $$b$$ so that gmres(A, b) achieves those residual norms at each step. Essentially everything you heard about the relationship between spectrum and GMRES convergence can fail spectacularly.

Some more references here http://www.cs.cas.cz/duintjertebbens/Talks/Liblice12.pdf .

• I think it's worth pointing out that if a matrix is normal then GMRES does have a clear dependence on the distribution of its eigenvalues. Furthermore, mild departure from normality is often still acceptable for understanding GMRES behavior based on eigenvalues - but ultimately it will fail. Nov 4 '19 at 21:50
• @Reid.Atcheson True, good point. Nov 4 '19 at 22:06
• That’s an astonishing theorem by Greenbaum et al. and the right answer to this question. But it reminds me of Trefethen’s inverse yogiisms (people.maths.ox.ac.uk/trefethen/berra_published.pdf) - statements that are true but may give the wrong impression. Mar 9 at 7:28

Estimates for GMRES convergence based on eigenvalue distribution often implicitly assume that the matrix is normal. Sometimes the convergence rate is still provable in an asymptotic sense in the non-normal case, but if the matrix is severely non-normal then the "pre-asymptotic" behavior will make such convergence rates never reachable in practice.

Your matrix is not normal, and its non-normality gets worse as the size grows.

• Normality alone is not enough, "fast convergence occurs when the eigenvalues of A are clustered away from the origin and A is not too far from normal" -- W GMRES. A simple testcase: poisson 1d - 2 I. Mar 17 at 13:58
• "Eigenvalue distribution" in comments above was meant to include potential clustering. Mar 17 at 18:41

Well, I can only propose a potential cause for why GMRES method fails for the problem you showed. I don't have enough reputation so I can't comment.

Since GMRES is using Anorldi Process to generate a set of orthogonal vectors in the Krylov subspace, if the matrix itself is very large, that means the Q and H matrix produced could be very large. So if you are solving a very large problem, you need to use the limited-GMRES (or limited-memory GMRES, can't recall the name exactly). But the key point is that the solver will reset the orthogonal matrix and restart the entire process to handle the memory limitation. This is probably the reason why your GMRES failed. Check if there is a limited option in the MatLab command.

A\b works straight away since this is a very sparse and bi-diagonal system, which is simply perfect for backward elimination.

Hope this helps. Cheers!

• Thanks Mr Sen, you mean the restarted GMRES, which is used to handle the memory problem. many thanks. Nov 8 '19 at 0:45
• @Zhen-WeiSun, sorry for the wrong terminology. That is the one I was referring too. Hope that would help.
– Sen
Nov 8 '19 at 2:02
• That's ok, I get it. hope to communicate with you often ;-). Nov 8 '19 at 2:22
• No problem. Cheers🍺
– Sen
Nov 8 '19 at 5:35

First some comments on why such matrices are hard for solvers. Notice that if $$U$$ is upper diagonal like

[[. 2 . . .]
[. . 2 . .]
[. . . 2 .]
[. . . . 2]
[. . . . .]]


then $$(I - U)^{-1} = I + U + U^2 + ... U^{n-1} \ (U^n = 0)$$. For this example, $$(I - U)^{-1} =$$

[[ 1  2  4  8 16]
[ .  1  2  4  8]
[ .  .  1  2  4]
[ .  .  .  1  2]
[ .  .  .  .  1]]


has corner $$2^{blksize - 1}$$. Huge entries will cause trouble for any solver -- either floating-point overflow or underflow, or extreme ill-conditioning. With $$U$$ random-Cauchy instead of random-normal (to amplify the problem), the singular values of $$I - U \ \ 300 \times 300$$ are

[6.7e-18 1.2e-06 2.1e-05 0.0025 0.0026 0.005 0.0063 0.0068 0.0083 0.011 ... 18 22 22 25 30 32 35 37 45 610]


and the pseudo-inverse is not triangular. (This is with python numpy and scipy, I don't have Matlab.)

Do Krylov methods work for at all for matrices which cannot be diagonalized, such as these $$I - U$$ with eigenvalues all 1 ? What-is-the-principle-behind-the-convergence-of-krylov-subspace-methods-for-solving linear systems of equations ?, says clearly "IF A is diagonalizable ...". But if not, how can one go about finding a polynomial for which $$p(A) \, b \approx A^{-1} \, b$$ ? No idea -- over to experts.

See also the many papers by Trefethen et al. and the book Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators 2005, 606 pages.

• +1 for pseudospectra, this is the way out for non-normal operators. Mar 11 at 18:16