Restarting is a filter

The Lanczos essay ends with three costs, all of them growing with the number of steps. The basis is m vectors and they have to be kept. Orthogonality is lost geometrically over about ten consecutive steps and never recovered. And the method returns twenty-five extra copies of thirteen eigenvalues in eighty steps on a forty-by-forty matrix, each copy accurate to 1.9·10⁻⁸ and none of them a second eigenvalue.

Restarting is the standard answer to all three, and the standard description of it — run m steps, keep the k best pairs, start again — makes it sound like a housekeeping measure. It is not. The restart is where the method does its aiming, and the aiming is a filter in exactly the sense the regularisation field’s filter factors are one — the same polynomial a stopping test is a race is about, aimed at a different object.

How well it separates is a function of how many shifts it is given, and the function is not gradual.

Across one to six shifts the worst wanted amplification reads 1.26, 1.13, 1, 0.803, 0.32 and 0.279, the best discarded one reads 1.18, 0.973, 0.802, 0.548, 0.00213 and 1.18·10⁻⁵, and their ratio — which is the whole of what a restart buys — reads 1.07, 1.16, 1.25, 1.47, 150 and 23,600.

Four gentle steps and then two cliffs. From one shift to four the separation improves by 37 per cent in total; from four to five it improves by a factor of a hundred, and from five to six by another hundred and fifty-seven. The polynomial is a product of (λ − shift) factors, so each extra shift multiplies every direction by its distance from a new root — and once the roots surround the unwanted part of the spectrum, one more of them costs the discarded directions two orders of magnitude apiece.

And the wanted directions stop being amplified at three shifts. The first two columns read 1.26 and 1.13 — above one, so the filter is growing the directions it means to keep, which is harmless and is not what it was sold as — and from three shifts on they read 1, 0.803, 0.32 and 0.279. A restart with enough shifts damps everything and merely damps the unwanted part far harder, which is a different mechanism from the one the word “filter” suggests.

One number degrades as the rest improve, and it is worth naming. The agreement between the measured amplification and the polynomial evaluated directly reads 2.2·10⁻¹³, 1.2·10⁻¹², 3.91·10⁻¹¹, 3.33·10⁻¹¹, 2.5·10⁻¹⁰ and 2.76·10⁻⁹ across the six — four orders of magnitude worse at six shifts than at one. That is the polynomial’s own dynamic range: a product whose value spans 10⁵ across the spectrum cannot be evaluated to the accuracy of one whose value spans 1.1, and the two-route check degrades accordingly rather than the method doing so.

The mechanism, written out #

Run m = k + p steps of Lanczos from a starting vector v. That produces m Ritz values; keep the k wanted ones — the largest, here — and call the other p unwanted. Then form

v⁺  =  Πⱼ (A − θⱼ I) v  /  ‖·‖

with the θⱼ the unwanted values, and start again from v⁺.

Write v in A’s eigenbasis as Σ cᵢ xᵢ. Each factor (A − θI) multiplies the ith component by (λᵢ − θ), so the whole product multiplies it by Πⱼ(λᵢ − θⱼ) — a polynomial evaluated at λᵢ, and the same polynomial for every component. Directions whose eigenvalues sit near a discarded Ritz value are multiplied by something small; directions near the wanted values are multiplied by something large.

That is what the figure shows and it is checkable because the matrix was built from its own eigenvectors: withSpectrum takes a list of eigenvalues and a fixed random orthogonal V and returns VΛVᵀ. So the components can be measured directly rather than inferred — the exact ground truth habit applied to a filter — and the assertion is that every one of them is scaled by the polynomial times a single shared normalisation — checked at all twenty eigenvalues, not at the extremes.

What the implicit method does instead, and why the difference is not visible here #

The version above forms the product explicitly: p products with A per cycle, and a starting vector that has been multiplied by a polynomial in A.

The implicitly restarted method reaches the same vector without forming any of it. It applies p steps of shifted QR to the m×m tridiagonal, which — by the implicit Q theorem — performs the same filtering on the Krylov basis at a cost that does not involve A at all. That is cheaper and it is the same vector.

Drawing the explicit version is a deliberate choice, and this site has made it before. The Francis double shift draws the implicit step and then checks it against the explicitly formed (A − μI)(A − μ̄I) = QR, entry by entry, agreeing in absolute value to 2.2·10⁻¹⁵. The explicit form is the argument; the implicit form is an implementation of it, and separating the two is what lets a reader see what the machinery is for before meeting the machinery.

What the filter buys #

Eight vectors. An unrestarted run reaching the same accuracy on this matrix keeps every vector it has generated, and the whole reason the method exists is that at n = 10⁶ the difference between eight vectors and eighty is the difference between a method that runs and one that does not.

The comparison is drawn in products with A rather than in cycles, because that is the currency both methods spend: a cycle of m = 8 steps costs eight products plus p = 4 for the filter, so twelve cycles is 144 products. An unrestarted run of 144 steps on a 40×40 matrix would exhaust the space five times over — which is the other half of why restarting is not merely a memory measure.

Two properties of the convergence are worth separating.

The residual bound is free. βₘ|sₘ,ₖ| needs the last entry of the tridiagonal’s own eigenvector and no product with A at all. It is what a real code stops on, and it is what is plotted; the true error is drawn beside it only because this matrix was constructed and the answer is known.

And the progress is not monotone at every cycle, which is why the assertion is stated as a ratio across the whole run rather than step by step. A restart is not required to improve the worst wanted Ritz pair at every cycle, and claiming it does would be claiming something false. What it does instead, past the point of convergence, is worse than non-monotone and is the subject of the next section.

Past convergence the bound diverges, and the answer does not #

Two endpoints are not a trajectory. Read cycle by cycle, and carried out to twenty:

The method converges at cycle three. The true error is at the rounding level from there on and stays between 1.8·10⁻¹⁵ and 1.4·10⁻¹⁴ for the remaining seventeen cycles — a factor of eight, with no trend in it. Whatever happens after cycle three does not happen to the answer.

The bound bottoms out at cycle five and then rises, monotonically, for the rest of the run. By cycle twelve it is sixteen times worse than its own minimum; by cycle twenty it is sixteen thousand times worse. That is not non-monotonicity — it is a steady divergence, and a caller watching it without a reference solution would read it as a method coming apart.

The mechanism is the one this collection keeps finding at the bottom of a converged iteration. Once v spans the wanted invariant subspace, the p unwanted Ritz values are extracted from a Krylov space that is numerically deficient, so the filter’s roots are rounding noise; and βₘ, the last off-diagonal of a tridiagonal whose recurrence should have terminated, is the size of a residual that ought to be zero. The bound βₘ|sₘ,ₖ| is then a ratio of two noise quantities, and nothing holds it down.

So the twelve-cycle figure runs seven cycles past its own best answer, and the 1.1·10⁻¹³ quoted from it is sixteen times worse than the 6.6·10⁻¹⁵ the run reached at cycle five. Neither number is wrong and the pair of them is misleading, in the direction this collection is usually alert to and this time reversed: the reported quantity is pessimistic and the answer is better than it says.

It also puts a number on the parameter the essay says is chosen as a memory budget. Twelve cycles was picked to fill a figure; five is what the problem needed, and the four products a cycle past that are spent making the reported quantity worse. On a forty-by-forty matrix that is 96 wasted products out of 236 — forty per cent of the run — and on the n = 10⁶ problems the method exists for it is forty per cent of everything. The cycle count is not a memory budget at all; it is a stopping rule, and it has to be written as one.

The practical instruction is the ordinary one and the measurement is why it is not optional. Stop when the bound stops falling, not when it reaches a tolerance — because past convergence it will not reach one again, and a loop written as while bound > tol on this problem at tol = 10⁻¹⁵ runs forever while holding an answer correct to fifteen digits from its third cycle. assertTheBoundDivergesPastConvergenceAndTheAnswerDoesNot holds the early convergence, the minimum inside six cycles, the monotone rise from the tenth onwards, the four orders of divergence, and the flatness of the true error across the whole run.

The ghosts, as a side effect #

lanczos.js measures what the plain three-term recurrence does over eighty steps: orthogonality lost by twelve orders of magnitude, and twenty-five duplicate eigenvalues that are accurate to 1.9·10⁻⁸ and are not there.

A restarted method with full reorthogonalisation inside each cycle cannot do that, and the reason is arithmetic rather than clever: the basis is eight vectors, reorthogonalising against eight vectors costs eight projections a step, and orthogonality lost geometrically over ten steps has no ten steps to be lost over. The ghosts are a consequence of a long recurrence, and there is no long recurrence.

So restarting solves three problems with one mechanism, and only one of the three is what it is usually introduced for. That is worth noticing because it changes what the parameters are chosen for: k and p are usually described as a memory budget, and they are also the degree of the filter and the length of the recurrence.

Choosing the shifts, which is choosing the polynomial #

The shifts here are the unwanted Ritz values, which is the standard choice and is called exact shifts. Read as a filter, it says: put the polynomial’s roots where the method’s own current estimate says the uninteresting part of the spectrum is.

That is a good choice and it is not the only one, and the filter reading is what makes the alternatives legible rather than arbitrary. Shifts placed on a line segment covering the unwanted region — Leja points, Chebyshev points — are choosing a polynomial small on an interval rather than at m points, and they behave differently when the unwanted spectrum is a continuum rather than a few isolated values. None of that is measured here; what is measured is that the mechanism is a polynomial, at which point the design space is the space of polynomials and not a list of tricks.

The measurement in the figure says how well the exact-shift choice does its job on this matrix: the worst wanted direction is amplified by 1.005 and the best discarded one by 0.802, a ratio of 1.25 in one cycle. That is not a dramatic separation — and it does not need to be, because the cycle repeats and the ratio compounds. Twelve cycles of 1.25 is a factor of fourteen, and the bound in the figure above falls by thirteen orders of magnitude over those twelve cycles, because the Ritz values sharpen as the vector improves and each cycle’s polynomial is better aimed than the last one’s.

The connection this essay exists to make #

Three filters have now been measured on this site and they act on three different objects.

The first two are the previous field’s subject and they suppress components of an answer. This one suppresses components of an input, and what it is aiming at is not accuracy but attention — which directions of the space the method will spend its next m steps on.

They are the same mathematics in the sense that matters: a polynomial in A applied to a vector, with the design question being where to put its roots. The site’s identical-algebra-different-arithmetic thread is usually about two computations of one quantity; this is one computation appearing in two fields with different names, which is the shape a synthesis phase is for.

What a restart cannot do #

One limitation belongs here rather than in the what is left, because it is the subject of the next essay and because it is a consequence of the filter reading rather than an unrelated fact.

Everything above multiplies the components of a single starting vector. Whatever the polynomial does, the space the next cycle explores is the Krylov space of one vector — and a repeated eigenvalue’s eigenspace contributes exactly one direction to that, no matter what v is. So a restart can aim the method at a part of the spectrum and cannot make it see a multiplicity.

The refusal this essay’s family runs makes the same point from the other side: fed a starting vector whose component along the largest eigenvalue has been removed entirely, the claim that a Krylov method finds that eigenvalue anyway is required to fail. It does. A filter can suppress a direction, and it cannot create one.

What the machinery had to be checked against #

ritzPairs writes the Lanczos recurrence out a second time, because a restart has to start from a given vector and this site’s own lanczos starts from a seeded random one. A second body for a routine the site already has is exactly what kit_dup_check exists to be suspicious of, and the suspicion is right: a copy that drifts is how a measurement ends up describing a bug.

So it is held against lanczos.js on the one input where the two must agree — same matrix, same starting vector, same number of steps, full reorthogonalisation in both — and the Ritz values agree to 5.3·10⁻¹⁵. That check is cheap and it is load-bearing: without it the convergence in the figure above could be measuring the copy rather than the method.

The filter measurement needs a second piece of machinery for the same reason. The components of a vector along each eigenvector are only available because withSpectrum builds the matrix from a prescribed spectrum and a fixed random orthogonal V — the same move exact.js makes with the Hilbert inverse and fft.js with the Kac–Murdock–Szegő one. On a measured matrix the components would have to be estimated by an eigensolver, and comparing an eigensolver’s output against a polynomial would be comparing two approximations.

The residual bound, which is what a real code stops on #

Every convergence claim in this essay is drawn against βₘ|sₘ,ₖ| rather than against the true error, and the reason is that the bound is free: it needs the last entry of the tridiagonal’s own eigenvector and no product with A at all. A code that has no exact answer to compare against — every code — stops on that number.

The true error is drawn beside it only because this matrix was constructed. What the pair shows is how much slack the bound has: at the last cycle the bound is 1.1·10⁻¹³ and the error is 7.1·10⁻¹⁵, so the quantity a real run reads is about fifteen times pessimistic there. That is the sort of factor a stopping tolerance has to be set with in mind, and it is the sort of number that is almost never reported beside a convergence curve.

Choosing k and p, which is choosing three things at once #

The parameters are usually presented as a memory budget: keep k, work in k + p, and pick p as large as the machine allows. The filter reading says that the same two numbers are doing two more jobs.

p is the degree of the filter. Each shift is a root, so a larger p is a sharper separation between the wanted and unwanted parts of the spectrum per cycle — and costs p products with A per cycle to apply.

k + p is the length of the recurrence, which is what decides how much orthogonality is lost inside a cycle before it is thrown away. A restarted method with a very long cycle has the ghost problem back.

And k is what the answer is. Asking for four eigenvalues and keeping four is a different method from asking for four and keeping eight, because the four extra Ritz pairs carry information the next cycle’s filter is aimed with.

Nothing here measures the trade between them, and the natural experiment — the same total number of products with A, spent as many short cycles or few long ones — is a figure this field could have and does not.

What is left #

Thick restarting, which keeps the wanted Ritz vectors explicitly rather than re-deriving them from a filtered starting vector, and is what most current codes do. It reaches the same subspace by a different arithmetic route, and the interesting comparison would be their behaviour under rounding rather than their exact-arithmetic equivalence.

Restarting for interior eigenvalues, where the exact-shift choice is much less obviously right — the unwanted spectrum surrounds the wanted values on both sides, and a polynomial small everywhere except a window in the middle is a harder object than one small at one end.

And the count of restarts as a parameter. Everything here fixes k, p and the number of cycles by hand. The step count became a regularisation parameter in another field this phase; whether the cycle count behaves the same way — with an interior optimum, and a rule for finding it — is a measurement nobody here has made.

The same subspace applied to a function #

A restart chooses a polynomial to apply to a starting vector. Every Krylov method is doing that; for a matrix function the polynomial is chosen to approximate f on the spectrum, and the machinery is otherwise identical.