GMRES and Nonsymmetric Systems

Where conjugate gradient runs out of road

The previous guide handed you the conjugate gradient method, and it felt almost magical: solve A x = b using only matrix-vector products, store just a handful of vectors, and watch the error melt away. But read the fine print. Conjugate gradient leans on the matrix being symmetric positive definite (SPD) — symmetric so that A^T equals A, positive definite so that x^T A x > 0 for every nonzero x. That second property is what makes 'solve A x = b' equivalent to 'minimize a bowl-shaped energy', and the bowl is what CG slides down.

Now meet the systems that break this. Discretize a flow that drifts in one direction — say heat carried downstream by a moving fluid — and the matrix picks up a lopsided term: A is no longer symmetric, because the physics has a preferred direction and the math remembers it. The same happens for many circuit equations, for the linearization inside matrix-free nonlinear solvers, and across most of engineering. The energy bowl tilts into a saddle or a spiral, there is no minimum to slide toward, and conjugate gradient simply has no theorem to stand on.

The Krylov idea survives — but the cheap recurrence does not

What carries over is the engine, not the trick. The Krylov subspace is the same: starting from the residual r_0 = b - A x_0, you keep multiplying by A to build the growing space spanned by r_0, A r_0, A^2 r_0, and so on up to A^{m-1} r_0. After m steps this space K_m holds every vector you can possibly reach with m matrix-vector products, and the best approximate solution will live inside it. That part needs no symmetry at all; it is pure linear algebra.

What does not carry over is conjugate gradient's beautiful secret: a short, three-term recurrence. For an SPD matrix, each new search direction needs to be made orthogonal only to the previous one or two — the older ones take care of themselves, thanks to symmetry — so CG runs on a fixed, tiny amount of memory no matter how many steps it takes. That economy is a direct gift of A = A^T. Strip away the symmetry and the gift vanishes: a new direction is no longer automatically orthogonal to the old ones, and there is no shortcut to making it so.

So a nonsymmetric Krylov method faces a stark fork. Either keep a short recurrence and give up the guarantee that each step is genuinely optimal (this is the family of methods like BiCGStab), or insist on a genuinely optimal step and pay by storing and orthogonalizing against every previous direction. GMRES takes the second road. It refuses to forget — and that refusal is both its strength and its cost.

What GMRES actually minimizes

Here is the whole idea of GMRES in one sentence: among all candidate solutions that live in the Krylov subspace K_m, pick the one whose residual is smallest in the ordinary 2-norm. The name spells it out — Generalized Minimal RESidual. CG minimized a special energy norm that only made sense for an SPD matrix; GMRES drops that and minimizes the honest, universal quantity ||b - A x||_2, the length of the leftover. Because that quantity makes sense for any square matrix, GMRES works for any nonsingular A, symmetric or not.

But minimizing over a subspace needs a tidy description of that subspace, and the raw Krylov vectors r_0, A r_0, A^2 r_0, ... are a terrible basis: they all tilt toward the dominant eigendirection and quickly become nearly parallel — numerically almost identical, like trying to describe a room with three rulers all pointing the same way. The cure is the Arnoldi iteration: it builds an orthonormal basis q_1, q_2, ..., q_m for the same Krylov space, one new perpendicular vector per step, using modified Gram-Schmidt to subtract off the parts that point in directions you already have.

Arnoldi does something lovely on the side: as it orthogonalizes, it records the coefficients in a small (m+1)-by-m matrix that is upper Hessenberg — upper triangular plus one extra subdiagonal. This little matrix H is a compressed shadow of the giant A, projected onto the Krylov space. The minimization over the huge space then collapses to a tiny least-squares problem involving H: minimize ||beta e_1 - H y||_2 over the short vector y, where beta = ||r_0||. You solve a problem of size m, not of size n, and m is usually in the dozens while n can be in the millions.

GMRES, one cycle of m steps:
  r_0 = b - A x_0 ;  beta = ||r_0||_2 ;  q_1 = r_0 / beta
  for j = 1, 2, ..., m:
      w = A q_j                          # the ONLY matrix-vector product
      for i = 1..j:  h_ij = q_i . w ;  w = w - h_ij q_i   # orthogonalize
      h_{j+1,j} = ||w||_2 ;  q_{j+1} = w / h_{j+1,j}
      solve  min_y || beta e_1 - H y ||_2   # tiny (j+1)-by-j least squares
      stop when residual ||beta e_1 - H y|| is small enough
  x_m = x_0 + Q_m y                        # combine the orthonormal basis

The catch: cost that grows every step, and the restart

Now the honest accounting. Because GMRES keeps every basis vector, step m must orthogonalize the newest vector against all m previous ones. So the work at step m grows like m, the total work after k steps grows like O(k^2 n), and the memory grows linearly — you are holding k full-length vectors at once. For a matrix with millions of rows, holding even a few hundred such vectors is real memory, and the orthogonalization can come to dominate the matrix-vector products. CG had none of this; GMRES pays for generality with a bill that climbs every iteration.

The standard escape is restarted GMRES, written GMRES(k): run k steps, form the best solution x_k in that small Krylov space, then throw the whole basis away and start fresh with x_k as the new initial guess. Memory and per-cycle cost are now capped at k. But there is no free lunch — discarding the accumulated subspace can stall convergence or even freeze it entirely, because the method forgets directions it had paid to discover. Choosing k is a genuine tuning knob: too small and it stagnates, too large and you are back to drowning in memory.

What governs convergence — and what comes next

For CG the convergence story was clean: faster when the SPD matrix's eigenvalues are clustered, slower when they are spread out, all captured by the condition number. For GMRES the story is genuinely murkier, and this is one of the honest hard truths of the field. A nonsymmetric matrix need not even have a sensible set of orthogonal eigenvectors, so eigenvalues alone do not determine how fast GMRES converges — you can build matrices with identical eigenvalues yet wildly different convergence. The real driver involves how the eigenvectors themselves are conditioned, and in the worst case GMRES can sit nearly flat for many steps and then plunge.

The practical takeaway is the same one the next guide will hammer home: raw GMRES on a hard system is often hopeless, and the cure is not a cleverer Krylov method but a preconditioner — an easily-invertible approximation M of A so that you actually solve the transformed system M^{-1} A x = M^{-1} b, whose eigenvalues are squeezed into a tight cluster near 1. A well-preconditioned GMRES can converge in a dozen steps where the unpreconditioned version would never finish. The Krylov machinery is the chassis; the preconditioner is the engine that makes it move.

1. Is your matrix symmetric positive definite? If yes, use conjugate gradient — short recurrence, tiny memory, no contest. GMRES would only waste storage here.
2. If it is nonsymmetric, reach for GMRES — but expect to restart it as GMRES(k) to cap the growing memory and per-step orthogonalization cost.
3. Always supply a preconditioner. Unpreconditioned GMRES on a realistic problem is the usual reason people conclude 'iterative solvers do not work' — when in truth they skipped the step that makes them work.
4. Watch the residual norm, not the iteration count, against an honest stopping tolerance — and remember every digit you keep is bounded by conditioning times round-off, not by the solver's effort.