Every essay — page 9
The eigenvalue problem that is not linear
A damped structure does not produce Ax = λx. It produces (λ²M + λC + K)x = 0, where the matrix whose null vector is wanted is a function of the number being solved for — so there is nothing to factorise, an n × n problem has 2n eigenvalues, and every algorithm anybody runs is an algorithm for something else. What that substitution costs is measurable: the linearisation is solved to the rounding level at every stop of a sweep across which the answer loses eleven orders.
Reduction, and what a model is for
A state matrix is n × n and the thing anybody wants from it — the map from an input to an output — is a function of one complex variable. The gap between those two numbers is this field. A reduced model is an approximation of the function and not of the state, so its error is measurable as a norm of a function rather than as a residual, and one of the two ways of building it comes with a bound computable before the model exists. That bound turns out not to be an inequality: on ordinary problems it is attained. What makes any of it possible is that a Gramian's eigenvalues fall off a cliff, and the reason they do is a rational approximation problem with a closed-form rate.
A model that is a rational function
A state matrix has a hundred thousand rows and the thing anyone wants from it is a function of one complex variable. The number that says how much of that size was ever the complexity is a rank — and the rank a derivation writes down cannot be computed, while one built from samples alone can.
The bound that is known in advance
Almost every error in numerical linear algebra is measured after the fact. Balanced truncation has one that is computable before the reduced model exists, in a norm of a function rather than of a residual — and on ordinary problems it is not an upper bound that is loose. It is attained.
The product nobody had to form
The Hankel singular values are the square roots of the eigenvalues of PQ. Form that product and half of them stop existing, at a floor one number predicts — and the fix is the one the least-squares field has had since its first essay, arriving in a place with no least-squares problem in it.
Exact at the points that were named
Balanced truncation asks for nothing and bounds everything, at a cost no large model can pay. The other kind of reduction asks for r numbers, costs r solves, is exact at every one of them — and bounds nothing anywhere else. That trade is the whole of large-scale model reduction.
A basis that is the same subspace and not the same thing
The interpolation conditions are conditions on a subspace, so any basis of it will do. The one a derivation writes down reaches a condition number of 7.7·10⁹ by its eighth vector, and the rate at which it gets there is set by a number the user chose with no information.
Interpolating at the model’s own poles
One choice of interpolation points is not arbitrary — the mirrored poles of the model about to be built. It is a fixed point rather than a guess, and when it is reached it beats a method costing O(n³) — by 0.4 per cent, which is the honest size of the whole contest.
Why a Gramian can be truncated at all
Every method in this field rests on one fact nobody states the reason for — the eigenvalues of a Gramian fall off a cliff. The equation defining it has a rank-one right-hand side and no low-rank structure anywhere — and the answer's decay is a rational approximation problem with a closed-form rate.
Where to put the poles of a rational function
Three times in one field the same question arrives from different directions — ADI shifts, rational approximation of a square root, the decay of a Gramian — and it has one answer. Cluster them geometrically towards wherever the function is difficult, and the alternative that looks reasonable costs orders.
A model that cannot be run
A stable system, reduced by matching its transfer function at three points exactly, comes back with a pole in the right half plane at four of nineteen placements — and matches at all three points to 1.3·10⁻¹⁴ while doing it. The construction did what it promised.
A model with no matrices behind it
Twenty-four numbers — values of a transfer function at twenty-four points — produce a sixth-order model of a forty-state system that passes through every sample to 10⁻¹² and matches the function it never saw to 10⁻⁹. Noise of 10⁻¹⁰ on those numbers takes the rank decision's gap from 10⁸ to 4.
A condition number that is not the model's
κ(P)κ(Q) is quoted as the reason one route to the Hankel singular values fails, and it is. It is also read as a measure of how reducible a model is, and it is not — the norm of the Gramian does not move at all as the McMillan degree runs from four to fourteen, and the condition number wanders over a factor of thirty-five with no trend.
Bracketing an error nobody can measure
The error of a low-rank Gramian factor is the one quantity a caller cannot compute, because computing it needs the Gramian the factor exists to avoid forming. Two numbers that can be computed sit either side of it — a rational factor known before the run, and a residual known after — and they stay a factor of four apart across a fourfold change of size.
The definition asks for more of what defeats it
The rank of a p × p Hankel matrix of Markov parameters resolves a degree of p, and it needs 2p parameters to do it. Those parameters grow like the norm of the state matrix raised to their index, so the count that buys resolution is the same count that buys dynamic range. One model, four run lengths, and a spread that runs from 10¹¹ to 10¹⁰².
The state that is removed is not a mode
Balanced truncation removes one state and pays exactly twice one Hankel singular value. The natural reading is that the state removed is the model's least important mode and that the σ is that mode's own size — and on five systems that reading over-estimates by between 1.21 and 7.46, never once under.
Half the conditions and a certificate
A two-sided interpolatory reduction of a stable convection–diffusion system came back unstable at four placements of nineteen. The one-sided reduction built from the same kind of solves came back stable at all 152 placements measured across eight Péclet numbers, and not by luck: every pole of a Galerkin model lies inside the numerical range of the operator, whose edge here is exactly the diffusion term's largest eigenvalue, −0.986 at Péclet 10. The price is the slope conditions, and spent as six one-sided shifts instead of three two-sided ones it beats the stable two-sided model at seven placements of ten.
A certificate written in coordinates
Every one-sided reduced model of the convection–diffusion system was stable, because its poles cannot leave the numerical range and the range sat in the left half plane. Rescale the state by a diagonal matrix and nothing about the system changes — not its eigenvalues, not its transfer function, not one pole of any two-sided model — but the numerical range moves. At Péclet 40 a rescaling with condition number 260 pushes its edge to +6.6 and ten placements of nineteen return unstable one-sided models. Projecting in the inner product of a Lyapunov solution puts the certificate back in any coordinates.
The inner product the mesh already computed
A Galerkin reduced model is certified stable when the operator's numerical range sits in the left half plane, and the certificate belongs to the coordinates. The coordinates that broke it before were a constructed rescaling. Grade a convection–diffusion mesh towards its outflow boundary layer — the grading anyone resolving the layer would choose — and nodal values break it too: the range's edge is past zero at a ratio of ten and reaches +8.9 at a thousand, and at Péclet 100 seven reduced models come back unstable. Weight the projection by the cell sizes the discretisation already computed and every model at every grading is stable, with the range's edge back at −0.098.
The smaller cell downstream
Nodal values on a convection–diffusion mesh graded towards its outflow layer put the numerical range in the right half plane, and Galerkin reduced models then fail. A code that grades towards the layer usually also upwinds, and upwinding adds diffusion. It does not give the certificate back: at Péclet 100 and a grading of a thousand the range's edge falls from 8.94 to 6.78 and stays positive, and 18 reduced models fail where 23 did. Each row of the upwind stencil's symmetric part sums to half the difference between the inverse widths of the cell after the node and the cell before it, so wherever the cells shrink downstream the dissipation upwinding adds is outweighed by the coupling to the smaller cell next door.
A mode that rings is counted twice
Balanced truncation's bound, twice the sum of the discarded Hankel singular values, was attained at every order on the heat model — a formula for the error rather than a bound on it. On six lightly damped oscillators it is loose at every order but the last, by factors between 1.3 and 5, at every damping from 0.7 to 0.01. The reason is visible in the singular values: a mode that rings contributes two of them, nearly equal, and removing the mode costs about what removing one would. The last mode's bound over its error is 2.02, 2.00 as the damping falls to 0.03 and 0.01. The equality was a property of modes with one singular value each.
Where a ringing mode stops being a pair
A lightly damped mode puts two nearly equal Hankel singular values into balanced truncation, and the reduction error was read as twice the larger of the first pair cut — right to two per cent at damping 0.01, wrong by half at 0.1. Read instead as the pair's two values added, the error is right to within six per cent at every order up to damping 0.05, and the pair's split turns out to be twice the damping. What it cannot do is say when the reading fails. The split predicts nothing: the widest split at damping 0.1 has the best reading. The failures come from somewhere the pair does not look — the other dropped modes' response under the peak, and a cut between two singular values that agree to under one per cent, which left a reduced model with a pole a thousand times closer to the axis than any pole of the system it came from.
A tie is a gap smaller than its neighbours
A balanced truncation that cut between two Hankel singular values agreeing to under one per cent left a reduced model with a pole near the axis and twice the error the σ promised, and the warning drawn from it was a rule: never cut inside a tie. Over 1,452 reductions of twelve oscillators, all ten doubled errors are at ties — but three of them have gaps of about two per cent, which a one-per-cent rule lets through, and what they share is a gap under a fifth of the gaps beside it. Moving every such cut one order down never doubles and halves each doubled error, for a median cost of one per cent. It is not free everywhere: at damping 0.5 it costs up to two and a half times. And the near-integrator was never the harm: 32 reductions have one and 22 of them read their σ correctly. The cheapest test of all is one evaluation of the error at zero frequency, which picks out the ten doubled cuts and nothing else.
An error estimate made of samples
A Loewner model built from samples of a transfer function comes with no error bound, and the number every code reads to choose its order — the pencil's next singular value — is the obvious stand-in. On four systems sampled at 48 points it is the error of the order-k model for the first two orders, within a factor of 2.5, and then falls away from it by about a decade an order: by the fifth it is optimistic by 109, 181 and 457, and on one system by 6,700 at the sixth. An AAA fit of the same degree to the same samples has an error within a factor of six of the Loewner model's at every order, so the shortfall is not the projection's. The model's misfit at its own samples is optimistic the same way. What works is more samples: the error at the 47 midpoints between the samples is within a factor of 2.2 of the true error at every order on every system, and the four midpoints nearest the poles alone within a factor of three.
Iterating, instead of factorising
At scale nobody factorises, and the reason is not accuracy — it is that the factors of a sparse matrix are not sparse. What replaces elimination is a sequence of approximations, and the question changes shape: not what the residual of the factorisation is, but how fast the error falls and whether the thing you can measure tracks the thing you cannot. On the model problem every rate is known in closed form before anything runs.