A spectral radius that grows first

A rate that is known in advance measured the three classical stationary iterations against their spectral radii, on the model problem where those radii are known in closed form, and found the rates exact. Jacobi’s ρ is 0.99211 at n = 24 and the residual falls by 0.99211 a step, and the site checked it against a formula rather than against a fit.

That worked because the model problem’s iteration matrix is symmetric. This essay is about what the spectral radius promises when it is not.

ρ(A) < 1   ⟹   Aᵏ → 0

is true for every matrix, and it is the whole of what the spectral radius guarantees. It is a statement about a limit. It says nothing at all about the first hundred steps.

What is being drawn #

The matrix is the one from the eigenvalues that are not there: ρ = 0.8 on the diagonal, m above it, everything else zero. Its spectrum is {0.8} with multiplicity six, at every m.

At m = 0 it is normal, and ‖Aᵏ‖ = 0.8ᵏ exactly — asserted as an equality, and asserted to be monotone at every one of the first sixty steps. There is no transient, and that is the control.

At m = 2 the norm rises by four orders of magnitude first.

The eigenvalues are the same. The characteristic polynomial is the same. Every quantity a spectrum can offer is the same, and the behaviour is not.

Where the growth comes from #

Aᵏ for this matrix has a closed form. Write A = 0.8·I + m·N with N the shift; N is nilpotent with N⁶ = 0, and the identity commutes with N, so the binomial theorem applies and terminates:

Aᵏ  =  Σ_{j=0}^{5}  C(k, j) · 0.8^{k−j} · mʲ · Nʲ

The largest term is the last: C(k,5)·0.8^{k−5}·m⁵, a polynomial of degree 5 in k multiplied by a geometric decay. A polynomial times a decaying exponential rises and then falls, and the turning point is at about k = (n − 1)/(1 − ρ), which for n = 6 and ρ = 0.8 is 25. The measured peak is at step 24.

So the transient is not a numerical artefact and it is not a subtlety about norms. It is the binomial coefficients, and there are n − 1 of them.

Two routes to the peak #

The peak is a fact about the powers. There is a second quantity that bounds it, computed from somewhere else entirely.

The Kreiss constant is

K  =  sup_{|z| > 1}  (|z| − 1) · ‖(zI − A)⁻¹‖₂

— a supremum over the region outside the unit circle of how large the resolvent gets, weighted by the distance from the circle. It is exactly the quantity the pseudospectrum picture in the previous essay is a plot of, integrated in a particular way, and it never multiplies two matrices together.

The Kreiss matrix theorem brackets the peak with it:

K   ≤   sup_k ‖Aᵏ‖   ≤   e · n · K

Measured here: K = 6.76·10³, peak = 1.98·10⁴, and e·n·K = 1.10·10⁵. Both inequalities hold, the lower one with a factor of 2.9 and the upper one with a factor of 5.6.

Two routes to one number, one through the complex plane and one through repeated multiplication, sharing no arithmetic. That is the site’s oldest habit, applied to a quantity nobody thinks of as having two routes.

What the sampling can and cannot break #

The supremum in K is over an unbounded region and the computation samples it: circles of radius 1 + δ for δ over six decades, and forty-eight angles on each.

So what is computed is a lower bound on a supremum, and the direction matters. K ≤ peak is the half the sampling cannot break — a sampled K is smaller than the true K, so if the true one is below the peak the sampled one certainly is. The other half, peak ≤ e·n·K, is the theorem, and a sampled K makes it harder to satisfy rather than easier.

One detail of that sampling is worth recording because it looked like a finding and was not. For a normal matrix the ratio δ/(δ + dist) climbs to 1 only as δ → ∞, so a grid stopping at δ = 2 returns 0.909 for a constant that is exactly 1. The first version of this computation did stop there and reported 0.909 as the normal matrix’s Kreiss constant. It is an artefact of the grid and not of the matrix, and it is the one place in this essay where the control could have been read as a result.

And the transient has no ceiling #

One matrix with a transient is an example. What makes it a statement is that the peak has no bound in terms of the spectrum.

assertTheTransientHasNoCeiling measures the peak at m = 1 to 5 and finds 637, 1.98·10⁴, 1.50·10⁵, 6.29·10⁵ and 1.95·10⁶ — a factor of three thousand across a sweep in which the spectral radius is asserted to be 0.8 at every point.

The peak position does not move at all: step 24 at every m, because (n − 1)/(1 − ρ) has no m in it. So the parameter changes how high the excursion goes and not how long it lasts, which is the separation the closed form above predicts.

Which norm, and why it does not matter #

A reader’s first objection to a transient is that a norm is a choice, and a different norm might not show one. It is a good objection and the answer is precise.

For any matrix with ρ(A) < 1 there is a norm in which ‖A‖ < 1, so the powers decay monotonically in it. The construction is standard: diagonalise where that is possible and use the eigenvector basis, or take ‖x‖_S = ‖S⁻¹x‖₂ for a suitable S.

The catch is what S costs. The norm in which this matrix contracts is one whose unit ball is enormously elongated, and the equivalence constants between it and the Euclidean norm are exactly the size of the transient — κ(S) ≈ 10⁴ here. So the statement “there is a norm in which it decays” is true and buys nothing: converting a bound in that norm back into a statement about the quantity anybody measures multiplies it by κ(S), and κ(S) is the peak.

This is the same trade as diagonalising a non-normal matrix at all. Bauer–Fike bounds eigenvalue movement by κ(V)·‖E‖, and κ(V) is infinite for a defective matrix, which is why the previous essay needed the pseudospectrum. Here κ(S) is finite and is the answer to the question, which is a slightly better situation and not a different one.

The transient is basis-dependent and its size is exactly the price of the basis change that removes it. That is not a way out; it is the same number written twice.

The continuous version, which is where it is famous #

Everything above is about Aᵏ. The analogue for e^{tA} is the same phenomenon with the unit circle replaced by the imaginary axis, and it is where the subject is usually met.

e^{tA} → 0 if every eigenvalue has negative real part. The transient is bounded below by the continuous Kreiss constant

K  =  sup_{Re z > 0}  Re(z) · ‖(zI − A)⁻¹‖₂

and above by e·n·K, the same bracket with the same constant. The matrix here is a discrete example because the powers are cheap to compute exactly and the picture has integers on its horizontal axis; nothing about the mechanism is discrete.

The best-known instance is fluid: the linearised Navier–Stokes operator for plane Couette flow has every eigenvalue in the stable half plane at every Reynolds number, and the flow becomes turbulent anyway. The transient growth of e^{tA} reaches a factor proportional to Re², which is large enough at any interesting Reynolds number for a disturbance to leave the linear regime — so the eigenvalue analysis is correct and predicts the wrong thing, which is the site’s wrong-blame verdict applied to a whole literature.

What it costs in practice #

Three consequences, and they are the reason this is not a curiosity about matrix powers.

A stopping rule fires early. An iteration x_{k+1} = Ax_k + b monitored by its residual, stopped when the residual has fallen by a factor of 10⁻⁶, will not stop during the transient — it will stop after it, having spent twenty-four steps going backwards. Worse, an iteration monitored by the change between iterates can be told to stop at the peak, where consecutive iterates are large and close.

A linearised stability analysis is wrong about the physics. A flow whose linearised operator has every eigenvalue in the left half plane is asymptotically stable, and the transient growth of e^{tA} is what lets a small disturbance grow by a factor of a thousand first — long enough for nonlinear terms to take over and for the flow to become turbulent. That is the standard modern explanation for subcritical transition in shear flows, and it is this figure with e^{tA} in place of Aᵏ.

And an error bound built on the radius is not a bound. ‖Aᵏ‖ ≤ Cρᵏ is true for every matrix with some C, and C for this matrix at m = 4 is above 10⁶. A statement of that form with the constant left unnamed is not a statement about anything.

The relationship to the previous essay #

The pseudospectrum and the transient are the same object measured two ways, and the Kreiss constant is the bridge.

If Λ_ε(A) sticks out past the unit circle by an amount d, then K ≥ d/ε — because there is a point z outside the circle with ‖(zI − A)⁻¹‖ ≥ 1/ε and |z| − 1 ≥ d. So a picture that reaches outside the circle is a picture of a lower bound on the transient, readable by eye.

That is why the two essays are the same anchor. One draws the region; the other draws what the region costs.

What to measure instead #

If the spectral radius does not answer the question, something has to, and there are three candidates with different costs.

The norm of the first few powers. Directly measured, exactly right, and costs a matrix multiplication per power. On a small matrix that is the answer; on a large one it is the simulation the analysis was supposed to replace.

The numerical abscissa or its discrete analogue. ‖A‖₂ itself bounds the first step, and the largest eigenvalue of (A + Aᵀ)/2 gives the initial slope of the continuous version. Cheap, and it describes only the beginning: it says whether the growth starts, not how far it goes.

The Kreiss constant. Expensive relative to the abscissa and far cheaper than the powers on a large matrix, because a resolvent norm at one point is one linear solve rather than a product of k matrices, and the number of points needed does not grow with n.

The three answer different questions and a reader deciding between them is deciding how much of the curve they need. The spectral radius describes the right-hand end, the abscissa the left-hand end, and the Kreiss constant the peak between them — which is the part with the transient in it.

Why this is not a small-matrix curiosity #

The matrices in this essay are 6×6, which invites the reading that a real problem would be better behaved. The opposite is true, and the closed form says why.

The peak is set by the binomial coefficients up to C(k, n−1), so it grows with the size of the Jordan structure, and it turns over at (n − 1)/(1 − ρ), so it lasts longer the larger that structure is and the closer ρ is to one. A discretisation refined by a factor of two has a larger n and a ρ closer to 1, and both push in the same direction.

What is small here is n, which makes the transient smaller and shorter than a realistic one. A 6×6 reaching a factor of twenty thousand is the conservative version of the phenomenon rather than an exaggeration of it.

The stopping rule this breaks #

Worth spelling out, because it is the failure a reader is most likely to meet.

An iterative solver monitored by the relative change between iterates stops when ‖x_{k+1} − x_k‖ ≤ tol·‖x_k‖. During a transient the iterates are large and the differences between them are a small fraction of that largeness, so the test can fire at the top of the excursion — with the error at its worst — and the solver reports success.

A solver monitored by the residual does not have that failure and has a different one: the residual grows during the transient, so a run that would eventually converge looks like a run that is diverging, and a code that aborts on a rising residual aborts on a problem it would have solved.

Both are defensible rules and both are wrong here for the same reason: they read a local property of the sequence as a statement about where it is going, and for a non-normal operator the local property does not carry that information for the first (n − 1)/(1 − ρ) steps.

What is worth carrying #

ρ(A) < 1 is a statement about a limit and this site’s asymptotic verdict exists for exactly this shape of claim. The powers decay in the end, and “in the end” is doing all the work.

The growth is the binomial coefficients, so its height is set by how far the matrix is from normal and its duration by (n − 1)/(1 − ρ). Neither is visible in the spectrum.

And the peak is bracketed by a quantity computed in the complex plane, which is worth having because it is the only thing in this essay that can be computed without multiplying the matrix by itself twenty-four times — and on a large matrix that is the difference between a diagnostic and a simulation.