A geometric series whose ratio is a matrix

The geometric series is the yardstick against which every other series is measured. Its whole theory fits in a line: $1 + r + r^2 + \cdots$ converges exactly when $|r| < 1$, and then its sum is $1/(1 - r)$.

The algebra behind that line never uses the fact that $r$ is a number. Multiply the partial sum $S_n = 1 + r + \cdots + r^n$ by $1 - r$ and everything cancels except $1 - r^{n+1}$. That works equally well when $r$ is a square matrix $A$ and $1$ is the identity matrix $I$:

$$(I - A)(I + A + A^2 + \cdots + A^n) = I - A^{n+1}.$$

So if the powers $A^n$ shrink to the zero matrix, the partial sums converge, and their limit is $(I - A)^{-1}$. The matrix geometric series — called the Neumann series, after Carl Neumann, who used it for integral equations in 1877 — inverts $I - A$ by adding up powers of $A$.

The question is when the powers shrink. For a number it is whether $|r| < 1$. For a matrix there are two candidate answers, and the obvious one is wrong.

The scalar case, for comparison #

With a number, every step shrinks the remaining gap by the same factor. The total after $n$ terms falls short of the limit by exactly $r^{n+1}/(1 - r)$, a quantity that decreases steadily from the first term. There is nothing to be surprised by, and no detour. The square that half, a quarter and an eighth fill shows it at a glance: each piece is a fixed fraction of the one before, and the uncut corner shrinks in proportion.

A matrix does not act by one factor. It stretches some directions and shrinks others, and — the source of everything that follows — it can turn one direction into another. The size of a matrix, its norm, is the most it stretches any vector. Its eigenvalues are the factors by which it stretches the special directions it leaves alone. For a symmetric matrix the two agree: the norm is the largest eigenvalue in absolute value, because the direction of greatest stretch is itself an eigenvector. For a general matrix they need not, and the gap between them is the whole subject of this essay.

Powers that grow before they shrink #

The three matrices are:

$$\begin{pmatrix} 0.9 & 0 \\ 0 & 0.8 \end{pmatrix}, \qquad \begin{pmatrix} 0.9 & 4 \\ 0 & 0.8 \end{pmatrix}, \qquad \begin{pmatrix} 0.95 & 1 \\ 0 & 0.95 \end{pmatrix}.$$

The first two have the same eigenvalues, 0.9 and 0.8, because the eigenvalues of a triangular matrix are its diagonal entries. The first is diagonal, so it simply multiplies the two coordinates by 0.9 and 0.8, and its powers shrink from the start. The second has a 4 above the diagonal: it takes the second coordinate, shrinks it by 0.8, and adds four times it to the first. That shear is what the eigenvalues cannot see. The powers of the second matrix have top-right entry $4(0.9^n - 0.8^n)/(0.9 - 0.8) = 40(0.9^n - 0.8^n)$ — a difference of two decaying exponentials, multiplied by forty, which rises to about 10 before the decay wins. The matrix itself has norm about 4.2, and its sixth power stretches some vector more than ten times.

The third matrix has a repeated eigenvalue, 0.95, and a 1 above the diagonal. It cannot be diagonalised at all — it is a Jordan block — and its powers have top-right entry $n \cdot 0.95^{n-1}$, because multiplying by it adds a copy of the second coordinate into the first at every step, and those copies accumulate. That grows like $n$ for a while and then decays like $0.95^n$, peaking near $n = 19$ at about 7.6. Here the transient comes not from a large off-diagonal entry but from two eigenvalues that have merged.

In both cases the eventual decay rate is the largest eigenvalue — 0.9 per step for the second matrix, 0.95 for the third — exactly as for the diagonal one. The eigenvalues decide where the powers go. They say nothing about how far the powers travel first.

The series still arrives #

The partial sums converge for all three, because the powers eventually shrink geometrically, and a series whose terms eventually shrink geometrically converges by the comparison test. The transient costs a little: the shear’s partial sums lag behind those of the diagonal matrix for the first dozen terms and then close at the same rate. And the limit itself carries the mark of the shear. The inverse of $I - A$ for the diagonal matrix is $\mathrm{diag}(10, 5)$; for the shear it has a top-right entry of $4/(0.1 \times 0.2) = 200$. A small perturbation in the second coordinate is amplified two hundred times in the first.

That is the practical content of the picture. $I - A$ is invertible, and the series finds its inverse, whenever the eigenvalues of $A$ are inside the unit circle — but how large the inverse is, and how long the series takes to settle, depend on more than the eigenvalues.

For a number the size of the sum is also a function of the ratio alone: $1/(1 - r)$ is large only when $r$ is close to 1. For a matrix the inverse can be large even when every eigenvalue is far from 1 — the shear’s eigenvalues are 0.9 and 0.8, as far from 1 as the diagonal matrix’s, and its inverse is forty times larger in one entry. The amplification that makes the powers grow is the same amplification that makes the sum large.

The scalar shadow of a repeated eigenvalue #

The Jordan block’s powers, $n \cdot 0.95^{n-1}$ in the corner, are a familiar kind of term: a polynomial times a geometric factor. The scalar series $\sum n r^n$ has exactly that form, and it was the first test case for the comparison with a geometric series.

The bars rise before they fall — $1/2$, $2/4$, $3/8$ are $0.5$, $0.5$, $0.375$ — and they are eventually dominated by a geometric curve with any ratio between $1/2$ and $1$. The Jordan block does the same thing in matrix form: its corner entry is $n$ times a geometric term, rises while the $n$ is winning, and is eventually held under a slightly slower geometric curve. A $k \times k$ Jordan block contributes terms up to $n^{k-1} r^n$, so larger blocks rise higher and longer before they fall. Every matrix is similar to a combination of Jordan blocks — up to a change of coordinates it is a stack of them along the diagonal — which is why this is the general mechanism: the powers of any matrix are sums of polynomials times powers of its eigenvalues, and the polynomials are what the eigenvalues cannot see.

Gelfand’s formula #

The decisive quantity is the spectral radius $\rho(A)$: the largest absolute value of an eigenvalue. The series converges when $\rho(A) < 1$ and diverges when $\rho(A) > 1$, whatever the norm. The bridge between the spectral radius and the size of the powers is a formula of Israel Gelfand’s:

$$\rho(A) = \lim_{n \to \infty} \|A^n\|^{1/n}.$$

Whatever norm is used, the $n$-th root of the size of the $n$-th power converges to the spectral radius. For a single step the size may be large — the shear has norm 4.2 — but the $n$-th root averages the stretching over $n$ steps, and in the long run only the eigenvalues survive. The curves fall slowly: after 200 steps the shear’s is about 0.92 and the Jordan block’s about 0.98, each still a little above its limit, because a polynomial factor $n^{k}$ contributes $n^{k/n}$ to the $n$-th root, which tends to 1 only slowly.

Gelfand’s formula is what makes the convergence theorem work. If $\rho(A) < 1$, pick a number $q$ strictly between $\rho(A)$ and 1; eventually $\|A^n\|^{1/n} < q$, so $\|A^n\| < q^n$, and the terms of the series are held under a convergent geometric series. The transient is a finite number of terms at the start, and a finite number of terms cannot stop a series converging.

Inside, on and outside the circle #

The three cases are the scalar trichotomy, relocated. A rotation matrix scaled by $\lambda$ has eigenvalues $\lambda e^{\pm 0.3i}$, of absolute value $\lambda$. At $\lambda = 0.95$ the powers spiral inwards and the sums converge. At $\lambda = 1$ the powers go round the circle for ever without shrinking; the partial sums are bounded — the rotation’s powers average out — but they oscillate and have no limit, the matrix analogue of $1 - 1 + 1 - \cdots$. At $\lambda = 1.05$ the powers spiral outwards and the sums grow like $1.05^n$. The boundary case is exactly where the eigenvalues touch the unit circle, and nothing about the norm locates it.

The middle case deserves a second look, because it is where scalar intuition says the least. For a number on the unit circle other than 1 itself, the partial sums of $\sum r^n$ are bounded and oscillate, and methods that average them assign a value anyway — the Cesàro averages converge to $1/(1 - r)$. The rotation behaves identically: its partial sums trace a bounded loop, their averages converge to $(I - A)^{-1}$, and the inverse exists because 1 is not an eigenvalue. Only when an eigenvalue is exactly 1 does $I - A$ fail to be invertible, and then no summation method can rescue a formula whose answer does not exist.

Where the series is used #

The Neumann series appears wherever a quantity is built up by repeated application of a linear rule, and each application contributes less than the last.

Absorbing Markov chains. In a chain that stops, the transitions among the transient states form a matrix $Q$ whose powers give the chance of still wandering after $n$ steps. The expected number of visits to each state before absorption is $I + Q + Q^2 + \cdots = (I - Q)^{-1}$, the chain’s fundamental matrix. The series converges because absorption is certain, which is to say $\rho(Q) < 1$. Its partial sums have a meaning of their own: the first $n$ terms count the expected visits within the first $n$ steps, so the series is literally the running tally of time spent, and its limit the total.

Input–output economics. Wassily Leontief modelled an economy as a matrix $A$ whose entry $a_{ij}$ is the amount of good $i$ used to make one unit of good $j$. To deliver a final demand $d$, the economy must produce $d$, plus the inputs to make $d$, which is $Ad$, plus the inputs to make those, $A^2 d$, and so on: total output $(I + A + A^2 + \cdots)d = (I - A)^{-1} d$. The economy is viable exactly when the spectral radius of $A$ is below one — when each round of inputs is, in the long run, smaller than the last.

Iterative solvers. To solve a linear system $x = Ax + b$, start from any guess and repeat $x \leftarrow Ax + b$. The iterates are the partial sums of $\sum A^k b$, and they converge to the solution exactly when $\rho(A) < 1$. The transient in the figures is a real phenomenon in practice: an iteration can appear to diverge for dozens of steps and then converge, and the eigenvalues alone cannot say how long that stretch will last.

What the eigenvalues do not tell #

The picture of powers rising before they fall is the phenomenon numerical analysts call transient growth, and it has a precise explanation: the matrix’s eigenvectors are nearly parallel. For the shear, the eigenvectors for 0.9 and 0.8 are $(1, 0)$ and $(-40, 1)$ scaled — almost the same direction. A vector that is a small difference of two nearly parallel eigenvectors has large components along each; as the powers shrink those components at different rates, the difference is no longer small, and the vector grows before both components die away. Where the eigenvectors are perpendicular, as for a symmetric matrix, nothing like this can happen, and the norm and the spectral radius agree.

Gershgorin’s discs give one way to bound the eigenvalues from the entries, and so to guarantee convergence cheaply: if every disc lies inside the unit circle, the series converges. The shear’s discs are centred at 0.9 and 0.8 with radii 4 and 0, so the first disc reaches far outside the circle and the test says nothing, although the series converges. The test is sufficient, not necessary, for the same reason the norm is: it cannot distinguish a shear from a stretch.

What the figures leave out #

Only two-by-two matrices are drawn. The phenomena are clearest there, but they grow with size. A Jordan block of size $k$ has powers with entries of order $n^{k-1}\rho^n$, and for matrices of size a hundred, as arise in discretised differential equations, transient growth by factors of millions is possible with every eigenvalue comfortably inside the unit circle.

The matrices were chosen to show a transient. Most matrices picked at random, with eigenvalues inside the unit circle, have mild transients or none; the shear with a 4 above the diagonal and the Jordan block are extreme on purpose. The claim is not that transients are typical but that the eigenvalues cannot rule them out.

The norm is the operator 2-norm. Every size in the figures is the largest stretch of a unit vector, computed exactly for $2 \times 2$ matrices. Other norms give other numbers, and Gelfand’s formula says they all agree in the limit; the heights of the transient peaks depend on the norm, and the conclusion that there is a transient does not.

Floating point is used throughout. The powers and partial sums are computed in double precision, which for these small, well-scaled matrices is accurate far beyond the plots’ resolution. For badly conditioned matrices the computed powers themselves can be dominated by rounding error, which is one more thing the eigenvalues do not warn about.

Still open: how large the transient can be #

The general question — given the eigenvalues, how large can the powers get before they decay? — has no answer in terms of the eigenvalues alone, and the tools that do answer it are the subject of active work. Pseudospectra, the sets of complex numbers that become eigenvalues under small perturbations of the matrix, give bounds on transient growth in both directions, and Lloyd Trefethen and Mark Embree’s theory of them has turned “the eigenvalues are misleading” into quantitative statements. But sharp bounds for specific classes of matrices, and a good understanding of transient growth in the linearised equations of fluid flow, where it is thought to be the route by which smooth flow becomes turbulent below the threshold that eigenvalues predict, are still being worked out.

A cleaner mathematical question sits underneath. For a matrix of size $N$ with spectral radius $\rho < 1$ and norm at most some bound, how large can $\sup_n \|A^n\|$ be? The Kreiss matrix theorem answers it up to a factor that grows with $N$, and whether that factor can be removed — whether transient growth is controlled by a quantity independent of the matrix’s size — was settled only for special classes; in general the dependence on $N$ is known to be necessary, and its exact form is still studied.

The ratio that was never a number #

The formula $1/(1 - r)$ survives the passage to matrices intact, as $(I - A)^{-1}$, and so does the condition for it, provided the right quantity is compared with one. It is not the matrix’s size, which can be large while the series converges, but the largest of its eigenvalues — the number the powers eventually decay by. Between the first power and the eventual decay lies a region the eigenvalues cannot describe, where the powers may grow by a factor of ten or of a million, and where the geometric series, the most predictable object in analysis, takes a detour before it arrives.

The detour is the reason the matrix series is worth knowing separately from the scalar one. In every application above — a chain wandering before absorption, an economy’s rounds of inputs, an iteration closing on a solution — the answer is the same inverse, and the practical question is what happens on the way. The eigenvalues guarantee the destination. The shape of the matrix, and in particular how nearly its eigenvectors line up, decides the journey. For a number there is no journey to speak of, and that is the one thing the scalar yardstick could never have taught.