The bound that is known in advance

Almost every error on this site is measured after the fact. Run the method, compare against a closed form, report the difference. The three or four bounds this collection has met are each disqualified from being useful in a different way: the growth factor’s 2ⁿ⁻¹ is never attained, the conjugate-gradient rate describes a rate rather than a number, and Johnson–Lindenstrauss holds with a probability.

Balanced truncation has one that is none of those:

‖H − Hᵣ‖∞  ≤  2 (σᵣ₊₁ + … + σₙ).

It is computable before the reduced model exists. It bounds a norm of a function rather than the size of a residual. And — this is the part worth the essay — on ordinary problems it is not a loose upper bound. It is attained.

What the σ are #

They are the Hankel singular values, and getting to them costs two Lyapunov equations:

AP + PAᵀ + BBᵀ = 0,     AᵀQ + QA + CᵀC = 0,

whose solutions are the controllability and observability Gramians, followed by σₖ = √λₖ(PQ). Nothing new has to be written to solve them: an equation whose unknown is a matrix already built a Bartels–Stewart solver for AX + XB = C, and a Lyapunov equation is that with B = Aᵀ. This field asks the existing routine a question it was already able to answer, which is the whole point of having written it there.

What the two Gramians mean is worth one sentence each, because the bound is unreadable without them. P measures how much of the state space the input can reach; Q measures how much of it the output can see. A direction that is hard to reach and hard to see contributes almost nothing to the map from u to y, and is exactly what a reduced model should drop. σₖ = √λₖ(PQ) is the number that combines the two, and it is invariant under any change of state coordinates — which is what makes it a property of H rather than of how the model was written down, in the same way the McMillan degree is.

And an answer to check them against #

The model here is the second-difference operator, so A = VΛVᵀ with V orthogonal and Λ a closed form. In modal coordinates the Lyapunov equation becomes n² scalar equations:

P̃ᵢⱼ = −b̃ᵢ b̃ⱼ / (λᵢ + λⱼ),      b̃ = VᵀB.

That is a formula, not a solve — n² divisions and one similarity. The Bartels–Stewart route factorises two matrices and never mentions an eigenvalue, and the two agree to 3.5·10⁻¹⁴. The equation the solve claims to have solved is separately checked by substitution, at 1.1·10⁻¹⁴.

Both checks matter and they are different checks. The second says the solver did what it said. The first says the thing it solved is the object the formula describes. A field whose central claim is about a bound needs the numbers the bound is made of to be known rather than computed, and this is the strongest form of that available here.

The reduced model #

The balancing transformation is the change of coordinates in which P = Q = diag(σ) — every direction equally reachable and observable — and in those coordinates truncation is just dropping the last n − r rows and columns. In practice nobody forms it that way, for a reason that is the next essay and that has a familiar shape.

What comes back is a Petrov–Galerkin projection: Aᵣ = WᵀAT, Bᵣ = WᵀB, Cᵣ = CT with WᵀT = I. The two bases are built from square-root factors of the Gramians, and the number that decides both the accuracy and the conditioning of the construction is the same one — the smallest σ kept. That is worth noticing early, because it means truncating at a σ that has not yet fallen is both inaccurate and unstable, rather than one or the other, and the two failure modes cannot be traded against each other.

The measurement #

Take a twenty-state model whose transfer function has ten poles, reduce it at every order from one to seven, and measure ‖H − Hᵣ‖∞ properly — not by sampling, but by forming the difference as a system in its own right. H − Hᵣ is itself a transfer function, of the block-diagonal state matrix with the outputs subtracted, so every norm applies to it with no new code. That is what makes an error in this field measurable rather than estimated, and it is the reason the numbers below are worth reading.

r      measured error     2 Σ σ        ratio
1      1.421·10⁻²         1.421·10⁻²   1.0000
2      6.880·10⁻⁴         6.880·10⁻⁴   1.0000
3      3.953·10⁻⁵         3.953·10⁻⁵   1.0000
4      2.026·10⁻⁶         2.026·10⁻⁶   1.0000
5      9.944·10⁻⁸         9.944·10⁻⁸   1.0000
6      3.561·10⁻⁹         3.561·10⁻⁹   1.0000
7      1.181·10⁻¹⁰        1.181·10⁻¹⁰  1.0000

Seven orders of magnitude of error, and one number for the ratio. The bound is not a safety margin.

Where the equality comes from #

It is not a coincidence and it is not a property of this model. It comes from a one-state result and a sum, and separating the two is what makes the whole thing legible.

And it is not a property of one order either, which is worth establishing before the mechanism is given.

Equality, not inequality, at every order tested. The sweep reports bound ÷ error between 1.0000 and 1.0000 at three, four, five, six, seven and nine orders — six sweeps, and the widest and narrowest ratio agree to four decimal places in every one of them.

That is a stronger statement than a bound and it needs to be read as one. A quantity that is equal to the error is not a bound being attained by luck — it is a formula for the error, and the distinction matters for what a caller may do with it. A bound justifies a decision; a formula lets the decision be inverted. Given a tolerance the order can be read off directly rather than searched for, and no reduced model has to be built to find out what the next one would cost.

Removing exactly one state costs exactly 2σₙ. Measured on four systems that share nothing but the shape of the question — residues all of one sign, residues of both signs, a clustered pair of poles, a spectrum spanning three decades — the ratio of the measured H∞ error to the σ removed is 2.0000 every time.

The bound for a truncation of many states is that result applied n − r times, each step contributing twice the σ it removed, and the terms added up. So the sum can only be loose in how the removals interact, never in the single step. On a spectrum falling as fast as this one the tail sum is dominated by its own first term and the interaction has nothing to do, which is why the ratio is one.

And the sum is loose where the removals interact. On a model with two nearly equal Hankel singular values, removing both at once costs 1.52 times less than the two single steps suggest — so the looseness is in the interaction, which is where the derivation says it must be. Which property of the pair produces it is a second question, and the obvious answer is the wrong one.

What makes the single step exact #

The one-state result is the load-bearing one, so it is worth saying where it comes from rather than citing it. In balanced coordinates the state being dropped is a single direction that is equally reachable and observable, with σₙ measuring both. Its contribution to the map from input to output is a rank-one term whose largest gain over all frequencies is exactly 2σₙ — the factor of two is the round trip, in and out, and the equality is a property of the coordinates rather than of the system.

That is why the number comes back at 2.0000 on four systems with nothing in common. It is not an empirical regularity that happens to hold on well-behaved problems; it is what balancing is, and the four systems are there to show that nothing about the residues, the pole spacing or the spread of the spectrum enters it.

It also says where the equality can fail: the sum of the single steps over-counts whenever the removed directions’ contributions do not simply add.

Which pairs interact, measured #

The natural reading is that two states with nearly equal σ are two directions the balancing transformation cannot separate cleanly, so dropping them together differs from dropping them one at a time. Sixty three-state models — ten pole gaps, three choices of the third pole, and both signs for the second residue — say the spacing is not what decides it.

                        models   bound attained exactly   largest looseness
same-sign residues        30            30 of 30                1.0000
opposite-sign residues    30             1 of 30                2.8250

With same-sign residues the bound is attained exactly on every model, including ones whose second and third Hankel values are within a factor of seven. With opposite signs it is loose on twenty-nine of thirty, including ones whose Hankel values are nineteen times apart.

The two groups overlap on the axis the essay’s explanation uses: a same-sign model with σ₃/σ₂ = 0.144 attains the bound exactly, and an opposite-sign model with σ₃/σ₂ = 0.053 — three times further apart — is eleven per cent loose. So the ratio of the Hankel values does not decide it.

What does is whether the two removed modes’ contributions add or cancel at the frequency where the error peaks. With the same sign they add, the two single-step costs are simply summed, and the sum is exact. With opposite signs they partially cancel, the combined removal costs less than the two separately, and the bound over-counts by however much cancelled.

The spacing is a real second-order effect inside the opposite-sign group — the mean looseness is 1.20 for the closer half of the σ ratios and 1.86 for the wider half — so the essay’s mechanism is present and it is not the first-order one. And 1.52 is not the largest looseness this bound has anywhere: 2.83 is, in a sweep of sixty small models, and nothing in the argument suggests that is a ceiling either.

That last point matters more than the number. A bound whose looseness is a factor of two rather than a factor of 1.5 is still an excellent bound, and the hero figure’s 1.0000 at every order is still what this essay is for. What changes is the reason to expect it: the essay’s version says the bound is tight unless the σ happen to cluster, which is a condition a reader can check on the curve they already have. The measurement says it is tight unless the modes cancel, which is a condition the σ curve does not show at all.

So the practical rule is one line shorter and less reassuring than it was. The bound is exact when the removed modes reinforce and over-counts when they cancel, and the Hankel values do not say which. A reader who sees a fast-decaying σ curve and concludes the bound must be tight has concluded it from the wrong evidence — and will be right most of the time, because a fast decay means the tail sum is dominated by its first term and there is little for any interaction to do. Right for the wrong reason is the usual state of an unmeasured claim, and it is why this one held for as long as it did.

The lower bound is untouched by any of it. σᵣ₊₁ ≤ ‖H − Hᵣ‖∞ holds for every model in the sweep, with no condition on signs or spacing, because it is a statement about what any order-r model can do rather than about what this one does.

The lower bound, which is the other half #

σᵣ₊₁ is a lower bound too. No reduced model of order r — not balanced truncation’s, not anybody’s — has an H∞ error below the first Hankel singular value it neglected. Measured here at r = 4, the error is 2.10 times σ₅.

Put the two together and the picture is unusual for this site:

σᵣ₊₁  ≤  ‖H − Hᵣ‖∞  ≤  2 Σ(k > r) σₖ,

with both ends computed from one pair of Lyapunov solves and neither of them a heuristic. The question what order is worth truncating at? has an answer that can be read off a curve before any reduced model is built, and the answer is within a factor of about two of optimal.

Nothing else in this collection has that shape. Choosing without knowing is about a regularisation parameter and every published rule for choosing it is a rule of thumb; four knobs and one floor measures four such rules against each other and finds them landing within a few per cent because none of them knows anything the others do not. Here the number is not a rule of thumb. It is a bound, it is proved, and it is tight.

What it is a bound on, which is not a residual #

Worth stating plainly because every other error on this site is a residual or a forward error of a vector. ‖H − Hᵣ‖∞ is the largest, over all frequencies, of the difference between two functions — equivalently the largest possible ratio of output energy to input energy that the reduced model gets wrong. It is a statement about every input the model could ever be given.

That is a much stronger object than a residual. A residual says this computed answer nearly solves this problem. This says this smaller model behaves like the larger one under every excitation, to within this number. And it is available before the smaller model exists.

The price is visible in the same sentence. Getting it requires two Lyapunov equations, each of which is a Schur factorisation of an n × n matrix and therefore O(n³) — which is exactly what a model of a hundred thousand states cannot afford. Every large-scale reduction in use is the other kind, which asks for r numbers, costs r solves, and bounds nothing at all. That trade is the field’s central one and this essay is one half of it.

The decay, which is the assumption underneath #

Everything above assumes the σ fall. If they did not, 2Σσ would be of the same size as ‖H‖∞ at every order and the bound would say nothing.

They do fall, on this model and on nearly every model anybody discretises, and the reason is not obvious. A Gramian is defined by a matrix equation whose right-hand side has rank one; there is no low-rank structure anywhere in the statement. That the answer nonetheless has eigenvalues falling off a cliff is a theorem about rational approximation, and it has a rate with a closed form. Why a Gramian can be truncated at all is that essay, and it is the one that says which models this field works on before any of them is tried.

What the bound is not, and the three ways that matters #

Three readings of an a-priori bound are natural and all three are wrong here, so they are worth separating.

It is not a bound on the state. Nothing above says the reduced model’s internal variables resemble the original’s, and they do not: the balancing transformation mixes every state into every other, so xᵣ is a combination of all n components and has no interpretation as a subset of them. A reader wanting to know what happened to the temperature at grid point seven will not find it, and the bound says nothing about it. What is bounded is the map from u to y, which is what the model was for.

It is not a bound on the poles. The reduced model’s poles are not a subset of the original’s, and there is no statement that each of the r kept poles is close to one of the n originals. On the model here the reduced poles at order four sit at 9.9, 39.3, 87.8 and 154.6, which happen to be close to the four slowest original poles, and that is a property of a system whose modes are cleanly separated rather than a theorem. A model with clustered poles produces reduced poles that are near neither the cluster’s members nor its centre, and the bound continues to hold at exactly the same value.

It is not a bound anybody has to trust. This is the useful one. The reduced model is r × r; the difference system is (n + r) × (n + r); measuring ‖H − Hᵣ‖∞ on it is a search over frequency and costs nothing at these sizes. So the bound can be checked on any problem small enough to build one, which is how the equality above was found rather than assumed. On a problem too large to check, the bound is what remains — and the check has been done on enough small ones to say what kind of number it is.

The third point is the one worth carrying into a working habit. A bound that is only ever quoted is a bound nobody has calibrated. This one is a factor of one from the truth on the models drawn here and about a factor of two from the best possible reduced model of the same order, and both of those statements come from measurements that a user with a small test problem can repeat in a minute.

Where the two Lyapunov solves actually go #

One practical note, since the cost is the whole reason this method has a rival.

Both Gramians come from the same Schur factorisation used twice — A and Aᵀ have the same Schur vectors up to a transpose — so a careful implementation factorises once and solves two triangular systems. That halves the O(n³) and changes nothing about which side of the availability line it sits on: half of a hundred-thousand-cubed is still a hundred-thousand-cubed.

What it does change is the region where the comparison in interpolating at the model’s own poles is worth making. Up to a few thousand states the bounded method is simply available, and the argument for the cheap one is about convenience rather than necessity. The line is not sharp and it moves with the hardware, which is why this field is careful to price both in operations rather than in seconds.