Numerical Linear Algebra
where the algebra stops being the arithmetic
The matrix that is a graph

No subtraction, and no waiting

The power iteration keeps a positive vector's small entries correct because it never subtracts, and it pays for that in steps. On a Markov chain made of two clusters joined by a coupling ε it needs about 7/ε of them — 6,900 at ε = 10⁻³, 69,000 at 10⁻⁴ — and then stalls at a fixed point of its own rounded arithmetic about u/(4ε) from the answer; at ε = 10⁻⁶ it is still 45 per cent wrong after 200,000. Gaussian elimination needs no waiting and subtracts: on the same chains its error grows like 1/ε, and on a chain whose entries fall to 10⁻¹⁸⁴ it returns the small ones as noise, eight of them negative. The Grassmann–Taksar–Heyman elimination does neither. Every entry of every chain comes back to a relative 10⁻¹⁵, in a third of n³ operations, because the stationary vector is well conditioned under the perturbations it makes.
21 min read 5 figures Measured, not assumedExact ground truthTwo routes to a numberClaims that can fail

Worth reading first: An eigenvector that must not change sign · A ranking that is an eigenvector.

An eigenvector that must not change sign found the Perron vector of a clique with a long tail coming back from a symmetric eigensolver with negative entries, where Perron’s theorem says every entry is positive — and found the power iteration returning every one of them, down to 10−2510^{-25}, correct to several digits. The reason was a sentence long: the power iteration multiplies a nonnegative vector by a nonnegative matrix, forms only sums of positive numbers, and so can neither produce a negative entry nor lose a small one to cancellation. “A property guaranteed by a theorem is preserved numerically only by an algorithm whose arithmetic cannot violate it.”

A ranking whose order is not determined then found the price: the power iteration’s change between steps understates its remaining error by 1/(1−λ2)1/(1 - \lambda_2), and a chain whose second eigenvalue is near one converges slowly. Both essays were about one route that never subtracts and one that does. There is a third, and it is the one every Markov-chain code uses.

Three routes to one vector

A Markov chain’s stationary vector π\pi satisfies πP=π\pi P = \pi with ∑πi=1\sum \pi_i = 1; it is the Perron vector of the transposed transition matrix, positive whenever the chain is irreducible.

The power iteration multiplies by PP and renormalises, every step, from the uniform vector. Gaussian elimination with partial pivoting solves (I−P)Tπ=0(I - P)^{\mathsf T}\pi = 0 with one equation replaced by ∑πi=1\sum \pi_i = 1, as any linear-algebra library would. And the Grassmann–Taksar–Heyman elimination is Gaussian elimination on I−PI - P rearranged so that no subtraction happens anywhere: each pivot, which ordinary elimination would compute as one minus a diagonal entry, is computed instead as the sum of the off-diagonal entries in its row — the same number, written as a sum of positive probabilities — and every update adds a product of probabilities to a probability.

Three states are enough to see the rearrangement. Take the chain whose rows are (0.5,0.3,0.2)(0.5, 0.3, 0.2), (0.1,0.6,0.3)(0.1, 0.6, 0.3) and (0.4,0.4,0.2)(0.4, 0.4, 0.2). GTH eliminates the last state first: the probability of leaving it for the other two is s3=0.4+0.4=0.8s_3 = 0.4 + 0.4 = 0.8, and every path through it is folded into the remaining chain — the entry from the first state to the second becomes 0.3+0.2×0.4/0.8=0.40.3 + 0.2 \times 0.4 / 0.8 = 0.4, from the second to the first 0.1+0.3×0.4/0.8=0.250.1 + 0.3 \times 0.4/0.8 = 0.25. Then the second state: its probability of leaving for the first is 0.250.25, read directly off the folded chain. Ordinary elimination would have reached the same pivot as 1−0.751 - 0.75, one minus the folded diagonal; GTH never forms the diagonal at all. Back substitution is the same kind of arithmetic: π1=1\pi_1 = 1, π2=0.4/0.25=1.6\pi_2 = 0.4/0.25 = 1.6, π3=(0.2+1.6×0.3)/0.8=0.85\pi_3 = (0.2 + 1.6 \times 0.3)/0.8 = 0.85, normalised to 20/6920/69, 32/6932/69 and 17/6917/69, which is the stationary vector exactly. Every number in the computation is a sum, product or quotient of positive probabilities.

Each route is checked against the exact answer. The chains are stored in double precision, and their stationary vector is solved again in rational arithmetic from the stored numbers, with each diagonal of I−PI - P taken as the exact sum of its row’s stored off-diagonal entries — the chain all three routes are actually given, since a stored row sums to one only to rounding and on a nearly decomposable chain that difference matters.

Two families of 24 states, each built to stress one property. Two dense clusters of twelve states joined by a coupling — every row of the first sends a fraction 3ε3\varepsilon of its probability to the second, every row of the second sends ε\varepsilon back — so that the chain’s second eigenvalue is within about 4ε4\varepsilon of one and mixing takes about 1/ε1/\varepsilon steps. And birth–death chains whose stationary entries fall by a factor ρ\rho per state, so that the last is ρ23\rho^{23} of the first: 10−2310^{-23} at ρ=0.1\rho = 0.1, 10−18410^{-184} at ρ=10−8\rho = 10^{-8}.

Every entry, every chain

The figure at the top of the page is the whole of the cluster family. GTH holds every entry to between 3⋅10−163\cdot10^{-16} and 7⋅10−167\cdot10^{-16} relative error at every coupling from 10−110^{-1} to 10−1210^{-12}. Gaussian elimination’s worst entry error grows as the coupling shrinks — 4⋅10−154\cdot10^{-15} at 10−110^{-1}, 2⋅10−122\cdot10^{-12} at 10−410^{-4}, 3⋅10−73\cdot10^{-7} at 10−1010^{-10}, 2⋅10−52\cdot10^{-5} at 10−1210^{-12}, close to u/εu/\varepsilon. The power iteration, run until it stops moving or for 200,000 steps, is accurate at the strong couplings and has not converged at all from 10−610^{-6} down.

What Gaussian elimination gets wrong on the clusters is specific, and it is not the small entries — on this family no entry is small. At ε=10−6\varepsilon = 10^{-6} every entry of the first cluster is too large by the same relative 8.08⋅10−118.08\cdot10^{-11} and every entry of the second too small by the same 2.69⋅10−112.69\cdot10^{-11}: a quarter of the probability and three quarters, and 2.0⋅10−112.0\cdot10^{-11} of it moved from one to the other. At 10−1010^{-10} the transfer is the other way and four thousand times larger. The shape of the vector inside each cluster is exact to rounding; what is wrong is only how much of the total each cluster holds. That is the quantity that depends on the coupling — the shares are set by the ratio of the rates between the clusters, 3ε3\varepsilon against ε\varepsilon — and it is the one an absolute error of a rounding in a transition probability of size ε\varepsilon moves by u/εu/\varepsilon.

Every entry's relative error in the stationary vector of a 24-state birth–death chain whose entries fall by 10⁻⁴ per state, by the GTH elimination and by Gaussian eliminationThe smallest stationary entry is 10·10⁻⁹³. GTH's worst entry error is 7.5·10⁻¹⁶; Gaussian elimination's is 1.7·10⁷⁵. The power iteration reaches 10⁻¹² on every entry in 350 steps.ρ = 10⁻⁴GTH, worst7.5·10⁻¹⁶elimination, worst1.7·10⁷⁵smallest entry10·10⁻⁹³481216202410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹staterelative error (capped at 1000)red: Gaussian elimination; green: GTH; dotted: every digit gonea subtraction-free elimination keeps the small entries
Fig. 1 Every entry’s relative error on a birth–death chain whose stationary entries fall by 10⁻⁴ per state. Green: GTH. Red: Gaussian elimination, capped at a thousand. Drag the decay from 10⁻¹ to 10⁻⁸.

The decaying chains separate the two eliminations entry by entry. Gaussian elimination computes each entry to an absolute accuracy of about a rounding of the largest, so its relative error climbs one decade per decade of the entry’s smallness: right to 10−1610^{-16} on the first state, wrong in every digit by the fifth or sixth at ρ=10−4\rho = 10^{-4}, and from there on returning numbers around 10−1710^{-17} — noise — where the true entries run down to 10−9210^{-92}. At ρ=0.1\rho = 0.1 eight of those noise entries are negative. Measured in absolute terms, Gaussian elimination is doing exactly what it promises: its largest error on any entry is 2.5⋅10−152.5\cdot10^{-15} of the largest entry at ρ=0.1\rho = 0.1 and 1.1⋅10−161.1\cdot10^{-16} at 10−410^{-4} — a rounding of the biggest number in the problem, spread over every entry. The last three entries it returns at ρ=10−4\rho = 10^{-4} are all 1.7⋅10−171.7\cdot10^{-17}, the same number three times: not three approximations to three different tiny probabilities, but the floor the absolute error leaves. GTH’s error is flat across the whole vector, under 10−1510^{-15} on an entry of 10−18410^{-184} as on an entry of a half, which is the property an eigenvector that must not change sign found in the power iteration, delivered by an elimination.

On this family the power iteration is also accurate, and quick: the chain mixes fast, and it reaches 10−1210^{-12} on every entry in 250 to 650 steps. The decaying chain is the case the earlier essay measured, and it confirms it.

The power iteration pays twice

The clusters are the case it did not measure.

The power iteration on the two-cluster chain: the largest relative error in any stationary entry against the number of steps, at four couplingsε 0.01: best 5.2·10⁻¹⁶ within 200000 steps; ε 0.001: best 3.3·10⁻¹⁵ within 200000 steps; ε 10⁻⁴: best 2.3·10⁻¹³ within 200000 steps; ε 10⁻⁶: best 0.67 within 200000 steps.110¹10²10³10⁴10⁵10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1power-iteration stepsworst relative entry errorε = 0.01ε = 0.001ε = 10⁻⁴ε = 10⁻⁶the plateau before the fall is the slow mode, about 1/ε steps longit pays in steps, and then in digits
Fig. 2 The power iteration’s worst relative entry error against its step count on the two-cluster chain, at four couplings. Each curve sits near one for about 1/ε steps, falls, and levels off.

From the uniform vector the clusters hold the wrong shares of probability — a half each, where the answer gives the second three quarters — and moving probability between clusters is exactly the slow mode. So the error sits near one for about 1/ε1/\varepsilon steps before it falls. The iteration reaches a relative 10−1210^{-12} on every entry after 680 steps at ε=10−2\varepsilon = 10^{-2}, 6,900 at 10−310^{-3} and 69,370 at 10−410^{-4}: between 6.8 and 6.9 times 1/ε1/\varepsilon. At 10−610^{-6}, after 200,000 steps, it is still 45 per cent wrong; at 10−810^{-8}, 99 per cent.

That is the first price, and it is the one the ranking essay’s stopping test was about. The second is where it stops.

Where the power iteration stalls on the two-cluster chain — the error at the fixed point of its own rounded map — and the steps it takes to reach 10⁻¹², against the coupling ε, with the lines u/(4ε) and 7/εε 0.1: 10⁻¹² after 60 steps, stalls at 6.3·10⁻¹⁶ by step 200000; ε 0.01: 10⁻¹² after 680 steps, stalls at 1.9·10⁻¹⁵ by step 6000; ε 0.001: 10⁻¹² after 6900 steps, stalls at 3.3·10⁻¹⁵ by step 13000; ε 10⁻⁴: 10⁻¹² after 69370 steps, stalls at 2.3·10⁻¹³ by step 79000. u is the unit roundoff; 4ε is the chain's spectral gap.10⁻⁴10⁻³10⁻²10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110³10⁶coupling εfloor (errors), or stepswhere it stallssteps to 10⁻¹²u/(4ε)7/εboth lines have slope −1; the stall varies by a factor of ten about its linea slow chain costs the power iteration both ways
Fig. 3 Where the power iteration stalls, and the steps it takes to reach 10⁻¹², against the coupling, with the lines u/(4ε) and 7/ε.

The power iteration in floating point is a deterministic map on floating-point vectors, and it eventually reaches a fixed point of that map and stops moving. The fixed point is not the stationary vector. At ε=10−4\varepsilon = 10^{-4} the iteration stalls by step 79,000 at a relative error of 2.3⋅10−132.3\cdot10^{-13} and never changes again; at 10−310^{-3}, at 3.3⋅10−153.3\cdot10^{-15}. Measured on nine chains — five couplings from 3⋅10−33\cdot10^{-3} to 3⋅10−53\cdot10^{-5} and three random cluster structures — the stall sits between 0.12 and 1.15 times u/(4ε)u/(4\varepsilon), where 4ε4\varepsilon is the chain’s spectral gap. Each step’s rounding nudges the vector by about uu, and in the slow direction the map only pulls it back by a fraction 4ε4\varepsilon per step, so the fixed point settles where the two balance.

So the power iteration’s subtraction-free arithmetic gives it entrywise accuracy only up to the rounding divided by the gap. On the earlier essay’s clique that was invisible, because that graph’s gap is large. On a nearly decomposable chain it is the whole story.

Its stopping test reads the change, not the error

A code running the power iteration does not have the exact answer to compare against. It stops when the vector stops changing, and a ranking whose order is not determined measured how badly that reads the error: the change between iterates understates the remaining error by 1/(1−λ2)1/(1 - \lambda_2). On these chains the second eigenvalue is 1−4ε1 - 4\varepsilon, so the factor is 1/(4ε)1/(4\varepsilon), and it is exact. Stopped at the first step where no entry changes by more than a relative 10−1210^{-12}, the iteration’s true error is 2.3⋅10−112.3\cdot10^{-11} at ε=10−2\varepsilon = 10^{-2}, 2.5⋅10−102.5\cdot10^{-10} at 10−310^{-3} and 2.5⋅10−92.5\cdot10^{-9} at 10−410^{-4} — 23, 249 and 2,500 times the tolerance, against 1/(4ε)1/(4\varepsilon) of 25, 250 and 2,500. A tolerance of 10−1210^{-12} buys nine correct digits at the weakest coupling measured and would buy three at 10−1010^{-10}.

A chain with no stationary vector found the same division necessary for a periodic chain, where the change never falls at all. Here it falls and lies. Neither failure is in the arithmetic; both are in reading a step’s change as a distance to the answer, which on a slowly mixing chain it is not, by exactly the mixing time.

Why GTH does not care

The last two results seem to point in opposite directions. Gaussian elimination’s error grows like 1/ε1/\varepsilon, which suggests a problem that is ill conditioned when the clusters are weakly coupled; GTH’s does not grow at all, which suggests one that is not. Both readings are right, about two different kinds of perturbation.

How far the stationary vector of the two-cluster chain moves, as a largest relative change in any entry per unit perturbation, when every transition probability is perturbed relatively and when it is perturbed absolutely, against the couplingε 0.1: relative 0.78, absolute 4.7; ε 0.01: relative 0.90, absolute 6.1; ε 0.001: relative 0.92, absolute 13; ε 10⁻⁴: relative 0.92, absolute 84; ε 10⁻⁶: relative 0.92, absolute 7900; ε 10⁻⁸: relative 0.92, absolute 7.9·10⁵; ε 10⁻¹⁰: relative 0.92, absolute 1.9·10⁹; ε 10⁻¹²: relative 0.92, absolute 8·10⁹. The perturbation is 10⁻¹⁰ times a standard normal, solved exactly.10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹10¹10³10⁵10⁷10⁹coupling εchange per unit perturbationrelative perturbationabsolute perturbationGTH's rounding is relative; elimination's is absolutethe problem is well conditioned in the right sense
Fig. 4 How far the exact stationary vector moves, per unit perturbation, when every transition probability is perturbed relatively (green) and absolutely (red), against the coupling.

Multiply every off-diagonal probability by 1+10−10z1 + 10^{-10} z, with zz standard normal, and solve exactly again: the stationary vector’s largest relative change is 0.780.78 to 0.920.92 times 10−1010^{-10} at every coupling from 10−110^{-1} to 10−1210^{-12}. Move every probability instead by 10−10z10^{-10} z times the largest one — an absolute perturbation, the same size everywhere — and the change grows like 1/ε1/\varepsilon, about 0.008/ε0.008/\varepsilon times the perturbation, until at 10−1210^{-12} the perturbation is larger than the coupling itself. On the decaying chains it is the same: relative sensitivity 3.8 to 3.9 at every decay, and absolute sensitivity 76 at ρ=0.1\rho = 0.1, 733 at 10−210^{-2}, 7.3⋅1047.3\cdot10^{4} at 10−410^{-4} and 7.2⋅1087.2\cdot10^{8} at 10−810^{-8} — between 7.2 and 7.6 times 1/ρ1/\rho. Each small entry is held up by a single small probability, the one that leads into it from the state before; an absolute perturbation of fixed size is a relative perturbation of 1/ρ1/\rho times that size on that probability, and the entry moves accordingly. A relative perturbation of the same probability moves it by the probability’s own relative change, which is why the deep tail costs GTH nothing.

The stationary vector is well conditioned under small relative changes to the transition probabilities, however slowly the chain mixes, and ill conditioned under small absolute changes. Which matters depends on what a method’s rounding amounts to. GTH forms every quantity as sums and products of probabilities; each of its roundings is a small relative change to a probability, and the measurement above says those move the answer by about their own size. Gaussian elimination subtracts, and its backward error is a small absolute perturbation of the whole of I−PI - P — the kind the vector is sensitive to as 1/ε1/\varepsilon. The same distinction runs through accurate is not a property of a method, where a bidiagonal singular-value algorithm computed tiny singular values to full relative accuracy because its rounding perturbed the matrix’s entries relatively and those singular values were relatively well conditioned. A Markov chain’s stationary vector is that situation for a nonsymmetric problem: the theorem that guarantees positivity also, in a quantitative form, guarantees insensitivity to relative change, and the elimination that perturbs only relatively inherits the guarantee.

The difference between the two eliminations is not their arithmetic’s precision. It is which of the two condition numbers their errors meet, which is the distinction two condition numbers of one matrix drew for a linear system’s scaling.

No waiting either

Arithmetic spent on the 24-state two-cluster chain against the coupling: the GTH elimination, Gaussian elimination, and the power iteration until it reaches 10⁻¹² or 200,000 stepsGTH about 4608 multiply-adds and Gaussian elimination about 9216, at every coupling. The power iteration spends 3.5·10⁴ at ε 0.1, 3.9·10⁵ at ε 0.01, 4·10⁶ at ε 0.001, 4·10⁷ at ε 10⁻⁴, 1.2·10⁸ (cap) at ε 10⁻⁶, 1.2·10⁸ (cap) at ε 10⁻⁸, 1.2·10⁸ (cap) at ε 10⁻¹⁰, 1.2·10⁸ (cap) at ε 10⁻¹².10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10³10⁴10⁵10⁶10⁷10⁸10⁹coupling εmultiply-addsGTHGaussian eliminationpower iterationpower-iteration points at the top hit the 200,000-step capthe accurate route is also the cheap one
Fig. 5 Multiply-adds each route spends on the 24-state cluster chain against the coupling: GTH, Gaussian elimination, and the power iteration until it reaches 10⁻¹² or 200,000 steps.

GTH is a direct method. It costs about n3/3n^3/3 multiply-adds whatever the chain — 4,608 at n=24n = 24 — half of Gaussian elimination’s, since it eliminates only one triangle. The power iteration costs n2n^2 per step, so at ε=10−4\varepsilon = 10^{-4} its 69,370 steps are about 40 million multiply-adds, nearly nine thousand times GTH’s, for an answer it stalls at 350 times less accurate. Even on the decaying chains, where the power iteration is at its best, its 250 to 650 steps are 140,000 to 370,000 multiply-adds — thirty to eighty times GTH’s — for a vector it stalls at between 1⋅10−151\cdot10^{-15} and 7⋅10−157\cdot10^{-15}, where GTH’s is under 8⋅10−168\cdot10^{-16}. The rate is the second eigenvalue is the reason it cannot be otherwise: an iterative method’s speed is set by the chain’s mixing, and an elimination’s is not.

That is the place GTH occupies. The power iteration never subtracts and has to wait for the chain to mix. Gaussian elimination never waits and subtracts. GTH is the one that does neither, and on every chain measured here it is both the most accurate route and the cheapest.

What two families of 24 states do not show

The chains are small and dense, so all three routes are run exactly as written, with no sparsity and no ordering decisions. On a large sparse chain GTH fills in like any elimination and the power iteration’s cost per step is proportional to the number of transitions, so the cost comparison turns on fill — the order decides the memory is the place that is priced. The power iteration is the plain one; an iteration that aggregates the clusters and iterates on the aggregated chain is the standard repair for slow mixing and is not measured. And every probability here is a normal double; the well on the far side of the band asked what happens to GTH’s guarantee when the small entries fall below the smallest normal number, and that question is untouched — every entry here, even 10−18410^{-184}, is normal.

Still open: the subnormal band, aggregation, and a third cluster

Below the smallest normal number. Push the decaying chain’s ratio to 10−1510^{-15} per state, so that its tail enters the subnormal range near state twenty. The prediction with a sign is that GTH’s relative accuracy fails exactly at the first entry below 2−10222^{-1022}, as the recurrence in the earlier essay did, and that its error there grows like the subnormal spacing over the entry — the band decides, not the elimination.

Aggregation. Iterative aggregation–disaggregation solves a small chain between clusters and the clusters separately, and is the standard cure for the power iteration’s 1/ε1/\varepsilon steps. The prediction is that on the two-cluster chain it reaches 10−1210^{-12} in a number of steps independent of ε\varepsilon — under thirty from 10−210^{-2} to 10−1210^{-12} — and that its stall point is independent of ε\varepsilon as well, because the aggregated chain is solved directly.

Three clusters in a line. Couple three clusters in a chain, the middle one to both ends by ε\varepsilon and the ends not at all. The prediction is that the power iteration then needs about 1/ε21/\varepsilon^2 steps rather than 1/ε1/\varepsilon — the slowest mode moves probability from one end to the other through the middle — while GTH’s error stays at 10−1510^{-15} and Gaussian elimination’s grows like u/ε2u/\varepsilon^2.

Where this goes next

1611162126313610⁻²²10⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹vertex, clique first then tail|entry| relative to the largestone rounding of the largest entrya bound that is provedPerron root11bracket, low11bracket, high11bracket width1.2·10⁻⁴smallest entry-1.5·10⁻²⁰entries below zero4every entry is positiveand the picture disagrees An eigenvector that must not change sign Perron's theorem says the leading eigenvector of a connected nonnegative matrix is strictly positive. On a clique with a long tail, four of its thirty-six entries come back negative — and beside them is the only two-sided bound in these essays that is proved rather than estimated. places settled, of 14the arithmetic14 of 14one missing link3 of 14α from 0.80 to 0.9011 of 14entries move by 4.05·10⁻¹¹entries move by 0.00443entries move by 0.00123three kinds of doubtcomparisons14smallest gap6.3·10⁻⁵largest gap0.0024arithmetic band4.1·10⁻¹¹one-link band0.0044α band0.0012the arithmetic settles nearly all of themthe data settles three A ranking whose order is not determined Of the fourteen adjacent comparisons in a top-fifteen, fourteen survive perturbing the arithmetic at the rounding level, eleven survive moving the teleportation parameter across its usual range, and three survive removing one link. The computation is the strongest part of the answer. 0173451688510211910⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1steps of the walkdistance from stationarya rate read as a timeλ₂ of the walk0.9steps measured122steps predicted131ratio0.93the dashed line is the eigenvaluethe curve is the walk The rate is the second eigenvalue A walk forgets where it started at a rate the graph's second eigenvalue names exactly. Across three orders of magnitude in the step count the prediction is five per cent high — and the published rate for PageRank is right for a reason nobody states, which is that a link graph is in pieces.

Named objects

A flat tag is an object no other essay names yet.

Componentwise condition numberGaussian eliminationMixing timePerron frobeniusPower iterationRelative accuracyStationary distribution