The product nobody had to form

The bound that is known in advance is a sum of Hankel singular values, and it is only worth as much as the smallest of them that can be computed. This essay is about how many of them can be, and the answer is: half, in digits, if the object is formed the way the definition writes it.

The definition writes σₖ = √λₖ(PQ), where P and Q are the two Gramians. So the obvious route is to solve the two Lyapunov equations, multiply, eigendecompose, and take square roots. It is four lines and it is what the derivation says.

What goes wrong, and where exactly #

The two curves agree for the first six values — to fourteen digits, which is not a coincidence but the statement that the two routes are the same algebra — and then the upper one stops descending. It flattens at 9.7·10⁻¹⁰ and stays there while the lower one continues to 10⁻¹⁶.

That floor is not a fudge factor and it is not a property of this model. It is σ₁√u, and the prediction is 2.3·10⁻⁹ against a measured 9.7·10⁻¹⁰.

Not a property of this model is a claim about other models, and it is checkable. Six of them, two families, three sizes and three McMillan degrees:

   model         σ₁         floor        σ₁√u       floor ⁄ σ₁√u
  heat n=20    5.65·10⁻²   1.94·10⁻¹⁰   5.95·10⁻¹⁰      0.326
  heat n=30    5.57·10⁻²   1.80·10⁻¹⁰   5.87·10⁻¹⁰      0.307
  heat n=40    5.53·10⁻²   1.19·10⁻¹⁰   5.82·10⁻¹⁰      0.204
  lowDeg r=4   1.44·10⁻¹   4.21·10⁻¹⁰   1.52·10⁻⁹       0.278
  lowDeg r=6   1.87·10⁻¹   6.43·10⁻¹⁰   1.98·10⁻⁹       0.326
  lowDeg r=10  1.82·10⁻¹   4.35·10⁻¹⁰   1.92·10⁻⁹       0.226

The form holds and the constant is not one. The ratio sits between 0.20 and 0.33 across a family whose σ decay smoothly and one with an exactly low degree, over a factor of three in σ₁ — so σ₁√u is the right quantity rather than merely the right shape, and the floor is about 0.28 σ₁√u rather than σ₁√u itself.

The direction of that correction is the comfortable one: the product route resolves three to five times further than the argument promises. What it does not change is the headline, because the quantity in the exponent is a square root — the floor is eight or nine digits below σ₁ on every model in the table, and a factor of three inside a square root moves that by half a digit.

The comfortable direction is not universal, though, and the table above cannot see why. Every row of it moves one thing: three sizes at a fixed degree, then three degrees at a fixed size. A real model has no reason to hold either fixed, and when the degree rises with the state dimension the ratio leaves the 0.20–0.33 band in both directions.

So 0.28 is a fit to that table and not a constant of the law. Across the same family with the degree following the size, the ratio runs 4.1, 0.61, 0.35, 0.35 at ten, fourteen, twenty-six and thirty-four states — an order of magnitude, straddling one. The half of the reading that survives is the one that was always the point: the floor is at σ₁√u to within an order of magnitude, which is what a square-root law can promise and all it can promise. The half that does not survive is the reassurance — at fourteen states the route resolves less far than the argument promises, not three to five times further, and a caller who had budgeted for the comfortable direction would be short by a factor of four.

That this figure’s own generator asserts only 0.05 < settled/floor < 20 rather than the band the table quotes is the same fact seen from the machinery: whoever wrote the assertion had already found that the ratio moves further than a fitted constant suggests.

That last observation is the reason a √u law is worth stating at all, and it deserves a sentence because it cuts both ways.

A prediction that is out by a factor of three is normally a poor prediction. Here the quantity predicted is a floor in σ, and σ is the square root of what the arithmetic actually damages, so a factor of three in the floor is a factor of nine in the eigenvalue error underneath it, and a factor of nine in an absolute error at the scale of σ₁² is unremarkable — it is a small polynomial in n from the Schur factorisation, exactly where one would expect it.

Read the other way: the square root that makes the floor so high is also what makes it so insensitive. A route that resolves σ down to σ₁√u would need the eigenvalue error to fall by a factor of 10⁸ to gain four digits, and no amount of care in the product’s formation delivers that. That is the whole argument for not forming it, and the measured constant of 0.28 is the size of what careful formation does buy: half a digit.

So the table above is a table of the improvement available inside the wrong route, and it is half a digit against the eight the right route gives. That asymmetry — the repair being worth eight digits and the tuning worth half of one — is the shape this collection keeps finding wherever a method squares something, and it is why the road that squares the problem ends with a change of algorithm rather than with advice about how to form the normal equations carefully.

The argument is two lines. An eigenvalue of PQ is σ², and a computed eigenvalue of a matrix of that kind carries an absolute error of order u‖PQ‖ = uσ₁². Taking the square root turns an absolute error of uσ₁² on σ² into a floor of σ₁√u on σ. So a route that forms the product can resolve Hankel singular values down to σ₁√u and no further — eight digits below the largest, where the other route reaches sixteen.

The condition numbers say the same thing more bluntly. On this model κ(P) and κ(Q) are each about 10⁹, so their product is 1.8·10³⁶ and their geometric mean is 1.3·10¹⁸. One of those is past 1/u = 4.5·10¹⁵, which is the point at which nothing small survives at all, and the other is not.

The shape, which this site has met before #

Two constructions, algebraically identical, one of which forms a product and squares the conditioning. There is an essay about that and it is the road that squares the problem: the normal equations AᵀAx = Aᵀb solve exactly the same least-squares problem as a QR factorisation, and below a computable value of ε they break outright.

What is worth noticing is that there is no least-squares problem here. There is no overdetermined system, no residual to minimise, no A to be tall. The argument transfers anyway, because it was never about least squares: it is about the difference between working with a matrix and working with its Gram product, and it applies wherever somebody forms one.

The repair transfers too, and it is the same repair. Do not form the product. Factor P = SSᵀ and Q = RRᵀ, take the SVD of RᵀS, and the σₖ come back as its singular values. The matrix RᵀS has condition number √(κ(P)κ(Q)) rather than κ(P)κ(Q), so its smallest singular value is resolved relative to its largest at u rather than at √u, which is exactly the difference between eight digits and sixteen.

That is the same accounting two condition numbers of one matrix does for a scaled system: one object, two ways of measuring how hard it is, and the answer to which one governs is decided by what the algorithm actually forms.

The factors, and why they are not Cholesky factors #

P is symmetric positive semi-definite, so the textbook square root is a Cholesky factor. It is the wrong tool here and the reason is the same fact that makes the whole field possible.

A Gramian whose Hankel singular values fall off a cliff is numerically singular by construction. That is not a defect: it is the property being exploited. A Cholesky factorisation of a numerically singular matrix either fails outright at a non-positive pivot or returns a factor whose trailing columns are noise, and neither outcome is a good foundation for a bound.

So the square root here comes from a symmetric eigendecomposition — P = VΛVᵀ, S = V√Λ, with the negative eigenvalues rounding produces clamped at zero and reported. That is the honest version of the same clamp: the clamp still happens, and the size of what was clamped is a number the routine returns rather than a decision it hides. The largest negative eigenvalue on this model is at the level of u‖P‖, which is what it should be.

What the difference is worth #

The floor is not academic, because every σ below it is a term in the error bound. A bound of $2\sum_{k>r}\sigma_k$ computed from the product route has, past r = 6, a tail made of numbers that are all 9.7·10⁻¹⁰ because the routine cannot see below that. The bound it reports is therefore too large — by a factor that grows with how many terms are floored — and a user reading it concludes that no further truncation is worth making when in fact the error would keep falling by orders.

That is a specific and quiet failure. Nothing errors. The σ come back in decreasing order, the bound is a decreasing function of r, the reduced models are correct at every order. The only thing wrong is that the prediction stops improving six orders before the reality does, and the two are distinguishable only by computing them another way.

It is worth putting beside an estimate that can be fooled, which is the same shape one field over: a quantity that is right until it is not, with nothing in the output to mark the transition.

Where the two routes agree, and why that matters #

They agree for the first six values to fourteen digits. That is the half of the measurement that makes it a measurement rather than a complaint.

If the product route were simply worse, the two curves would part immediately and the comparison would be about accuracy in general. They do not part immediately: above σ₁√u the product route is correct to the last digit, because there the absolute error uσ₁² is small relative to σ² and the square root does not amplify it. The routes part at exactly the value the argument says they should, and that agreement above the floor is what turns “the product is worse” into “the product resolves half the digits, and here is which half”.

This is the pattern two Gram–Schmidts established for this site and that six routes to one spectrum generalised: the interesting statement about two algebraically identical procedures is not that one is better, it is where they stop agreeing and whether that place was predicted.

The two bases, and why the same number governs both #

There is a second place the square roots appear, and it is worth following because it explains a constraint on the whole method that would otherwise look arbitrary.

The reduced model is a projection with two bases, W and T, and the square-root construction builds them as $W = RU\Sigma^{-1/2}$ and $T = SV\Sigma^{-1/2}$, where U and V are the singular vectors of RᵀS and Σ holds the σ being kept. The $\Sigma^{-1/2}$ is what makes WᵀT the identity, and it is where the conditioning of the construction lives: the smallest σ kept sets how large that scaling is.

So the smallest σ kept decides two things at once. It decides the error, through the bound, and it decides how well conditioned the transformation that produces the model is. Truncating at a σ that has not yet fallen is both inaccurate and unstable, and the two cannot be traded against each other — there is no order at which the model is poor but reliably computed, or accurate but shakily formed.

Measured, the biorthogonality ‖WᵀT − I‖ on this model runs 2·10⁻¹⁶ at r = 1 and 9·10⁻¹⁶ at r = 7, climbing to 9·10⁻¹⁴ at r = 9 where the σ being kept have reached the Gramians’ own accuracy. That is the same curve as everything else in this essay, seen from a third angle: the point at which the numbers stop being numbers is one point, and every quantity in the method reaches it together.

The order of the two decisions #

There is a practical ordering here worth stating, because it is not the one a reader arrives with.

The first decision is which route to take to the σ, and it is free: the square-root route costs one extra SVD of an n × n matrix on top of two Lyapunov solves that are already O(n³). It buys eight digits of the bound.

The second decision is where to truncate, and it is the one everybody thinks about. It is made by reading the σ curve — and the curve is only readable down to whatever the first decision left.

So the free decision governs the visible one, which is the arrangement this site keeps finding. The scaling that buys ten orders has the same shape in the polynomial field: two lines of scaling, chosen before anything is computed, worth more than every subsequent choice put together.

What it does not fix #

The square-root route is not exact and does not claim to be. Its own floor is the accuracy of the Gramians, which is where the σ curve above finally flattens at about 10⁻¹⁶ relative — the level at which the Lyapunov solves were performed. Below that the numbers are rounding and the sweep in the bound that is known in advance stops there deliberately rather than drawing an inequality between two noise floors.

And the conditioning of the Gramians themselves is not addressed by any of this. κ(P) ≈ 10⁹ is a property of the model — it is the Hankel decay, viewed from a different angle — and no choice of route makes it smaller. What the route decides is whether the algorithm pays it once or twice.

What a caller can check #

None of the above requires trusting this essay. There are two cheap tests that separate the routes on any model, and both are worth running once.

The first is the one drawn here: compute the σ both ways and look for the value at which they part. If the product route’s tail flattens, it has flattened at σ₁√u and the values below that are not values. A tail that is flat over more than two indices is the signature, and it is visible without knowing what the right answer is.

The second is to check the bound against the model it bounds, at one order. Build the reduced model, form the difference system, and measure ‖H − Hᵣ‖∞. If the measurement is far below the bound at an order where the bound’s tail is made of floored values, the bound is the thing that is wrong — which is the opposite of the usual reading, where a measurement below a bound is taken as the bound being conservative. Here it is diagnostic.

Both tests cost one extra computation each on a quantity already in hand, which is the ratio this field keeps producing: the check is cheap and the failure it catches is silent.

Why the definition writes it the other way #

It is worth asking why every derivation defines the σ through PQ if forming PQ is the thing not to do, and the answer is not carelessness.

The definition is a statement about invariance. Under a change of state coordinates x → Tx the Gramians transform as P → TPTᵀ and Q → T⁻ᵀQT⁻¹, so neither of them individually means anything — P alone can be made the identity by a choice of coordinates. The product PQ transforms as T(PQ)T⁻¹, a similarity, so its eigenvalues do not move. The σ are invariants of the system and PQ is the object that displays that, which is exactly what a definition is for.

The square-root construction computes the same invariants without ever writing down the object that makes their invariance obvious. RᵀS is not coordinate-free — it depends on which square roots were taken — but its singular values are, because a different choice of square roots differs by an orthogonal factor on each side.

That is a shape this site has met before and it is worth naming: the object that makes a property evident and the object worth computing are different objects. The number that decides nothing is the same observation about a determinant, which is the natural definition of singularity and a hopeless test for it, and the inverse that is never formed is the same about A⁻¹. Three instances now, in three fields, and the pattern is reliable enough to use as a prompt: when a quantity is defined through a product or a determinant, ask what it would take to compute it without forming one.

The prompt has a corollary that is easy to miss. The definition is not wrong and does not need replacing — it is doing a different job, and a derivation that avoided PQ to be numerically careful would be a worse derivation. What has to be separate is the sentence that says and therefore the algorithm is, which is the sentence this essay is about and which no amount of care in the derivation supplies.

The one-line summary #

Two Lyapunov solves, then either an eigendecomposition of PQ or an SVD of RᵀS. The first is what the definition says and resolves eight digits of Hankel singular value; the second is one extra factorisation and resolves sixteen. The floor of the first is σ₁√u, predicted and measured, and every value below it is a term in the bound the whole method exists to compute.