A Krylov space for a problem that is not linear

Every method in this field builds a subspace out of A, Ab, A²b, … and asks a small projected problem for its answer. The rate the condition number predicts prices the convergence, an orthogonalisation nobody calls one is about the sweep that keeps the basis usable, and an eigenvalue that arrives twice is what happens when it stops.

A quadratic eigenvalue problem has no A. It has three coefficient matrices, and the only matrix in sight is the 2n × 2n linearisation. So the obvious method is Arnoldi on that — and it stores vectors of length 2n whose bottom half is a scaled copy of the top half of an earlier one.

The recurrence that does not store the bottom half #

Write the linearised matrix as [[A, B], [I, 0]] with A = −M⁻¹C and B = −M⁻¹K. A Krylov vector is [qⱼ; pⱼ], and the recurrence’s next vector has top half Aqⱼ + Bpⱼ and bottom half qⱼ. The bottom half is therefore never new information: it is what the top half was, one step ago, with the same orthogonalisation coefficients applied to it.

The second-order Arnoldi recurrence stores the top halves and updates the bottoms with the same coefficients without orthogonalising them. Same subspace, n numbers a step instead of 2n.

Then the projection can be taken onto the quadratic rather than onto the linearisation:

Mₘ = QᵀMQ,   Cₘ = QᵀCQ,   Kₘ = QᵀKQ,

an m × m quadratic eigenvalue problem, solved by the machinery of the field’s first essay, with 2m eigenvalues.

So m stored vectors of length n yield 2m approximations, where Arnoldi’s m stored vectors of length 2n yield m. A factor of four in Ritz values per number held, exactly and by construction — measured at exactly four at every subspace size, because both halves of it are counting rather than arithmetic.

And it is closer, at every size #

The exactness of the storage claim makes it easy to assume the cheaper method must be worse somewhere. On this family it is not.

Measured on a chain of forty masses, distance from the dominant eigenvalue:

Closer at every size, and five times closer at the smallest. The mechanism is not subtle: more approximations out of one subspace is more chances at the eigenvalue wanted, and the projected quadratic extracts 2m of them where the projected linearisation extracts m.

The honest half is the absolute rate. Twenty-six vectors buy three digits, which is slow, and it is slow for a reason that has nothing to do with either method: the chain’s spectrum is an arc with no isolated extreme, and Krylov convergence is decided by separation. This is a comparison of prices at equal work, not a demonstration of a fast method.

The two extractions are different answers #

It would be tidy if the two methods produced the same Ritz values from the same subspace. They do not, and measured they differ by 9 per cent at six vectors and 32 per cent at twenty.

That is not a defect. Arnoldi’s Ritz values are the eigenvalues of the projected linearisation; the second-order route’s are the eigenvalues of the projected quadratic. Those are two different small problems asked of the same vectors, and there is no reason for their answers to coincide. What is required — and measured — is that both approach the same spectrum, and that the best Ritz value from each is a genuine eigenvalue of the quadratic to two digits or better at every size.

It is worth noticing which of the two is the more natural. Projecting the linearisation preserves the shape of the algorithm; projecting the quadratic preserves the shape of the problem. The small problem the second route solves is a quadratic eigenvalue problem, so everything the field has established about linearising, scaling and backward error applies to it — recursively, at size m, where it is cheap. A method whose subproblem is of the same kind as its problem is easier to reason about, and here it is also better.

What the subspace actually contains #

A step back, because the two extractions being different raises a fair question: what is the subspace, and is it the same one in both cases?

The linearised Krylov space is span{v, Av, A²v, …} in ℝ^{2n}, with v = [q₁; 0]. Write out its elements and every one of them is a pair whose top half lies in

span{q₁, Aq₁ + Bp₁, …}

and whose bottom half lies in the span of the previous top halves. So the projection of the linearised space onto its top block is exactly the space the second-order recurrence builds, and the second recurrence is not an approximation of the first — it is a representation of it, with the redundant half not stored.

That is what makes the storage claim exact rather than heuristic, and it is also what makes the orthogonality finding a real one. The two methods are working in the same space. What differs is the coordinates they hold it in, and one set of coordinates degrades.

Why this is where the field meets the rest of the site #

There is a practical reason this essay sits in the iterative field rather than in the polynomial one, and it is about the sizes involved.

Every measurement in the polynomial field so far is on a dense problem of size eight or twelve, factorised whole. The application the field comes from is not that: a finite-element model of a damped structure has n in the hundreds of thousands, M and K sparse, C often a low-rank modification of them, and the request is for twenty eigenvalues nearest zero. A 2n × 2n dense factorisation is out of the question, and everything the field has established has to survive the transition to a subspace method.

Most of it does, unchanged, and that is worth saying explicitly:

• The scaling matters as much. The projected quadratic inherits the coefficients’ imbalance, so the two lines of the scaling that buys ten orders go in front of the recurrence rather than in front of a factorisation.
• The backward error is the same quantity. η for the projected problem is η for a quadratic, and the residual for the large problem is the one that decides convergence.
• The count is still 2n. A subspace of dimension m gives 2m Ritz values out of 2n eigenvalues, and no subspace method changes the arithmetic of that.

What does not survive is the closed form: a real structure has no Lambert function and no decoupling sine basis, so the reference disappears and the residual becomes the only evidence. Which is exactly why the measurements in this field are made on a family that has one.

The shift, which is the other half of every real run #

Everything above starts the recurrence at the coefficients themselves, which converges to the largest eigenvalues — and the eigenvalues anybody wants from a damped structure are the smallest, since those are the modes that persist.

The standard repair is shift-and-invert: pick a σ near the region of interest and run the recurrence on Q(σ)⁻¹ rather than on the coefficients, which maps the eigenvalues near σ to the large end. That costs one factorisation of the n × n matrix Q(σ) — not of a 2n × 2n linearisation, which is the second place the second-order route saves — and every step after it is a solve with a factor already computed.

It also changes what the essay’s finding means in practice. With a shift, convergence to the wanted eigenvalues takes a handful of steps rather than twenty-six, which puts the useful subspace sizes below the size at which the basis collapses. The measurement above is therefore a statement about where the cliff is, made at a size where it is reachable; a shifted run in production stays on the safe side of it, and knows where the side is because someone measured.

The basis stops being a basis #

The price arrives, and it arrives suddenly.

Measured, ‖QᵀQ − I‖ for the second-order basis:

6 vectors     5.9·10⁻¹¹
10            6.1·10⁻⁷
14            1.6·10⁻³
20            1.35

Orthogonality gone entirely by twenty vectors. The linearised Arnoldi basis on the same problem stays at 2.4·10⁻¹³ throughout.

Both are orthogonalised by the same modified Gram–Schmidt sweep. So the difference is not the orthogonaliser — it is the recurrence. The second-order sweep orthogonalises the vector r against the qs, and applies the same coefficients to the auxiliary vector s without orthogonalising it. The auxiliary sequence therefore accumulates error freely, and it feeds the next residual. What is being orthogonalised at step j is a vector built out of an already-degraded auxiliary, and no amount of care in the sweep repairs an input.

This is the same finding an orthogonalisation nobody calls one makes about classical against modified Gram–Schmidt — same algebra, different arithmetic — one level up: same orthogonaliser, different recurrence, and a loss that the orthogonaliser cannot see.

And the answer keeps improving through the collapse #

Here is the part that makes it a measurement rather than a caution.

Across the range where orthogonality goes from 10⁻¹¹ to 1.35, the distance from the dominant eigenvalue falls from 1.4·10⁻² to 1.6·10⁻³. The answer gets better as the basis fails.

A run that watched only the answer would report steady progress at every step. There is no divergence, no stagnation, no residual plateau — the diagnostic that the residual the method reports recommends does not fire, because the Ritz residual is computed from the projected problem and the projected problem is consistent with itself.

What fails silently is everything downstream that assumes Q is orthonormal: the Ritz vectors are not orthogonal, the Ritz values may be repeated where the spectrum is not, and the projection is no longer a projection. That is precisely the situation an eigenvalue that arrives twice describes for Lanczos — a ghost, arriving because the basis lost the memory of a direction it had already covered — and it is the reason the fix in that essay is the fix here: re-orthogonalise, at a cost, or restart before the collapse.

What a code should do about it #

Two options, and the measurement prices both.

Restart before twenty. The collapse arrives at about the same number of vectors whatever the chain’s length — measured at n = 16, 24, 32, 40 and 56 — because it is driven by the recurrence’s own growth rather than by the size of the problem. So the safe subspace size is a property of the method, and it is small. That is the same discipline restarting is a filter sets out, with an unusually clean number attached to it.

Or orthogonalise the auxiliary sequence too, which is what the two-level variant of this recurrence does: keep an orthonormal basis for the pair rather than for the top halves, and represent both halves in it. That restores the orthogonality and costs back part of the storage saving — which is the trade in one sentence, and the honest way to end a comparison whose headline was that the cheap method wins.

The other thing a projected quadratic can do #

One consequence of projecting the problem rather than the matrix is worth drawing out, because it is not available to the linearised route at all.

The projected quadratic Mₘ, Cₘ, Kₘ inherits structure. If M, C and K are symmetric then so are their projections, and the small problem is a symmetric quadratic — which means the certificate of every eigenvalue real can be run on it for the cost of an m × m Cholesky, and the answer is a statement about the small problem that is often the statement wanted about the large one. If the coefficients are palindromic, the projection is palindromic, and the reciprocal pairing of a spectrum that comes in reciprocal pairs is available at size m.

The projected linearisation keeps none of that: the Hessenberg matrix an Arnoldi run produces is unsymmetric whatever the coefficients were, and the structure has been dissolved into it. So the second-order route’s real advantage may not be the factor of two in storage at all — it may be that the subproblem is still a member of the class the problem came from.

That is a claim this essay measures the storage half of and leaves the structural half of standing. It is the kind of loose end worth writing down rather than rounding off: the measurement here is about price, and the structural argument is about what a method can be told.

What is worth carrying into a large problem #

Three sentences, and the third is the one that is easy to lose.

The cheaper recurrence is genuinely cheaper and genuinely better at equal work, by a factor of two in storage, four in Ritz values per number held, and a measured factor of two to five in distance to the dominant eigenvalue at every subspace size tested.

Its basis fails before its answer does, at about twenty vectors whatever the size of the problem, and the failure is invisible in every quantity the method returns.

And the failure is a property of the recurrence rather than of the orthogonaliser. That is what makes it worth a measurement: a reader who knew that modified Gram–Schmidt is the good one, which this site established years of essays ago, would conclude that a basis built with it is safe. Here it is not, because what is being orthogonalised has already been damaged before the sweep sees it.

At other settings #