A problem with infinitely many eigenvalues

A quadratic eigenvalue problem is still a polynomial. It has 2n eigenvalues, a linearisation, and a finite answer that a finite computation can return. Everything in this field so far has been about what that answer costs.

Take one more step — let the matrix depend on λ through something that is not a polynomial —

T(λ) = A − λI + γe^−λI,

which is what a system with a delay produces, and three things stop being true at once.

There is no linearisation, because a 2n × 2n matrix has 2n eigenvalues and this problem has infinitely many. There is no characteristic polynomial, so the exact-degree route that counted the eigenvalues at infinity has nothing to interpolate. And “compute the spectrum” is not a request that can be granted at all.

The ground truth is a Lambert function #

The site’s habit is to measure against something exactly known — the standing this Hilbert inverse has and the reason it is used — and it would be easy to assume that a problem with infinitely many eigenvalues cannot supply one. This one can.

A is the model matrix, whose eigenvalues μ are the closed form used everywhere here, and the nonlinearity is a scalar times the identity. So the problem decouples in A’s eigenvector basis and each μ contributes one scalar equation

μ − λ + γe^−λ = 0.

Rearranged, (λ − μ)e^(λ−μ) = γe^−μ, which is the defining equation of the Lambert W function — and W has infinitely many branches, so

λ = μ + W_k(γe^−μ),      k ∈ ℤ.

The whole spectrum, in closed form, countably infinite, lying along two curves that run off to the left. Each branch value is computed by Newton in complex arithmetic and checked twice: by substitution into the Lambert identity, and by asking whether T(λ) is actually singular — measured at 9.9·10⁻¹⁵ over twenty-eight values, which is the check that matters, because it is about the original problem rather than about the identity used to solve it.

The restriction that makes this available is stated rather than hidden. A general delay problem has a matrix coefficient on the exponential, does not decouple, and has no such formula. What is bought by the restriction is the same thing the Hilbert matrix buys elsewhere on this site: a family where the answer is known and the computation can therefore be graded rather than compared.

The question has to change #

With infinitely many eigenvalues, “find the spectrum” is not merely expensive — it is ill posed. Every practical method for a nonlinear eigenvalue problem answers one of two narrower questions:

• Which eigenvalue is near here? — Newton, or an inverse iteration, from a starting guess.
• Which eigenvalues are inside this region? — a contour method.

Both are worth measuring, and the second is the one with an integer in it.

Newton, and the branch it does not find #

Newton’s method for det T(λ) = 0 needs no determinant. The logarithmic derivative of a determinant is a trace,

(log det T)′ = tr(T(λ)⁻¹T′(λ)),

so the Newton step is λ ← λ − 1/tr(T⁻¹T′), which costs n complex solves. Started 0.15 + 0.15i away from a known branch value it converges quadratically — measured order 1.91 over ten runs — and lands on an eigenvalue of the problem to twelve digits every time.

Five of the ten land on a different eigenvalue from the one they were started beside.

That is not a failure of the method and it is the central fact about the shape of the problem. Neighbouring branches for neighbouring μ are 0.22 apart in this family, and a start 0.21 away is not near enough to decide anything. There is no “the” eigenvalue here, and which of infinitely many is returned is a property of where the method was pointed rather than of the problem.

Compare that with the polynomial case, where the answer is a set and a solver returns all of it — and with the six routes that differ in accuracy while agreeing about which numbers they are computing. The whole apparatus of a small residual is not a small error assumes there is a right answer to be compared against; here the right answer depends on the question, and the question is a starting point.

The count is an integral and the integral is an integer #

The one thing that is finite is the number of eigenvalues inside a contour, and it is available without finding any of them. By the argument principle applied to det T,

N = (1/2πi) ∮ tr(T(z)⁻¹T′(z)) dz.

On a circle, with the trapezoidal rule, an analytic periodic integrand converges exponentially in the number of quadrature points — the same phenomenon that makes the trapezoidal rule the right choice for a periodic function, and the reason contour methods are practical at all.

Measured on a circle of radius two, where two eigenvalues lie inside:

4 points      5.7·10⁻²
8             1.5·10⁻²
16            5.3·10⁻⁴
32            3.6·10⁻⁷
128           1.4·10⁻¹³

Each doubling buys more than the last, which is what exponential convergence looks like on an arithmetic axis.

Then the answer is rounded to an integer, and the rounding is where the trap is. At eight points on that circle the rounded count is already right while the integral is 0.06 away — so a coarse quadrature can be correct by luck. And on a circle of radius four, which passes closer to the spectrum, eight points return a count of minus two where the answer is four.

A negative number of eigenvalues. From a quantity that is an integer whenever it has converged, by a routine that has not failed, thrown or warned. The only evidence available either way is the distance from the nearest integer — which is why that distance is the quantity drawn, and why a contour method that does not report it is reporting half of its answer.

The convergence rate is set by how far the contour is from the nearest eigenvalue, not by how many are inside. A circle passing close to one has a nearly singular integrand and needs several times as many points for the same accuracy, which is visible on the slider: at radius 1.5, which runs close to an eigenvalue at 1.587, a hundred and twenty-eight points reach only 7.6·10⁻⁴.

Two routes to the trace #

Every number above rests on tr(T(z)⁻¹T′(z)) being computed correctly in complex arithmetic, and this site does not take that on trust.

There is no complex kernel here and this field did not write one. A complex n × n solve is a real 2n × 2n solve — [[Re, −Im], [Im, Re]] against [xᵣ; xᵢ] — so the whole of the complex linear algebra is an embedding and a call to the same pivoted elimination elimination is a sequence of choices describes. The embedding doubles the size and squares nothing.

The second route is a closed form: since T is A plus a scalar times the identity, T⁻¹T′ has eigenvalues −(1 + γe^−z)/(μⱼ − z + γe^−z), and the trace is a sum of n scalars. The two agree to 5.6·10⁻¹⁵ at sixteen points around a contour. One is n solves of a 2n × 2n real system; the other is n divisions. Nothing is shared, which is what makes the agreement worth printing.

The ceiling nothing announces #

The practical contour method — Beyn’s — does more than count. It probes with a random n × ℓ block V and reads the eigenvalues out of the rank of ∮T(z)⁻¹V dz. That rank is at most ℓ.

Measured with four eigenvalues inside: the rank is 1, 2, 3, 4, 4 at ℓ = 1, 2, 3, 4, 5. A caller who asked for two gets two — with small residuals, each of them a genuine eigenvalue, and nothing in the returned object to say there were four.

And there is a third ceiling, which a reader would not predict. The block is n × ℓ, so its rank cannot pass n either. On a circle of radius five, which contains twelve eigenvalues of this four-by-four problem, the rank stops at four however many probes are used: the first moment simply cannot represent more than n eigenvalues, and recovering the rest needs higher moments rather than a wider block.

That is the same shape as the sketch that is spent: a random probe carries a fixed budget of information, and asking it for more than it holds returns something plausible rather than an error. Three ceilings, one number, and the returned object distinguishes none of them. The counting integral is the check — it costs one more pass round the same contour with the same integrand — and it is the difference between a method that answers the question and a method that answers as much of it as it happened to have room for.

Why there is no linearisation, precisely #

The impossibility is worth one paragraph of proof rather than an assertion, because it is what divides this essay from the six before it.

A linearisation of a matrix polynomial of degree d is a pencil of size dn whose eigenvalues are exactly the polynomial’s, elementary divisors included. Any pencil of any finite size has finitely many eigenvalues. T here has infinitely many — the Lambert branches are indexed by ℤ, and each of the n values of μ contributes one per branch. So no pencil of any size has T’s spectrum, and the device the whole field has been built on is simply unavailable.

What replaces it, in the methods that do exist, is approximation rather than equivalence: interpolate T by a polynomial or a rational function of some degree on a region of interest, linearise that, and accept that the eigenvalues found are the eigenvalues of the approximant. Which turns the problem into two error terms rather than one — the approximation of T, and the arithmetic of the linearised approximant — and the field’s whole ledger has to be redone with an extra column in it. That is a real subject and it is not this essay’s; what this essay establishes is the boundary, and the two questions that survive it.

The difference between infinitely many and merely a lot #

It is tempting to file this as “very many eigenvalues” and reach for the machinery that handles large spectra: compute a few, restart, filter. That machinery — which restarting is a filter sets out — rests on there being a finite set with a distinguished part of it, a largest, a smallest, a nearest to a shift.

Here there is no largest. Follow the branches leftwards and the eigenvalues run to infinity in both imaginary directions, with |λ| unbounded and no last one. “The dominant eigenvalue” exists only relative to a region; “the ten smallest” exists only relative to a point. Every request has to name a region before it means anything, which is why both surviving questions above have a where in them — near this point, or inside this circle.

That is the structural change, and it is bigger than the loss of the linearisation. A finite spectrum can be described; an infinite one can only be interrogated.

What this field ends on #

The polynomial field opened with an object that is not a matrix and ends with one that is not a polynomial either, and the arc is the same one this site keeps walking: each time the object gets further from a matrix, a quantity that was a property becomes a decision.

The number of eigenvalues was 2n, a fact. Then it was a degree, when M went singular. Now it is a contour integral rounded to an integer, with a convergence that has to be checked and a probe count that has to be chosen. An eigenvalue count that cannot be slightly wrong has a sequel and the sequel is that the count can be wrong by any amount, including a negative one, and that the only defence is to print the distance from the integer rather than the integer.

What a caller has to decide #

Collecting what this field’s last essay leaves in the reader’s hands, because it is more than the polynomial case leaves and the difference is the subject.

A region. There is no default. “The eigenvalues” is not a request, and every method needs either a point to start near or a contour to integrate around. For the delay problem the physically interesting region is the right half plane — an eigenvalue there is an instability — which is a region with an infinite boundary and therefore needs a truncation whose omission has to be argued for separately.

A quadrature resolution, and evidence for it. The distance from the nearest integer is that evidence, it is free, and it is the difference between a count and a guess.

And a probe count, which cannot be chosen without the count. Three ceilings bind it and only one of them is under the caller’s control.

Three decisions, none of which exists for a matrix eigenvalue problem, all of which have to be made before a number appears. That is what it costs to leave polynomials behind, and it is worth knowing before a model with a delay in it is written down rather than after.

At other settings #