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.

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 σ are close together. 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. That is the entire looseness of the bound, located: it lives in the interaction between removals of comparable size, and nowhere else.

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 explains the one place the equality does not survive. Two states with nearly equal σ are two directions the balancing transformation cannot separate cleanly, so dropping them together is not the same operation as dropping them one at a time. The measured 1.52 is the size of that effect on a deliberately clustered pair, and it is the largest looseness this bound has anywhere.

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.

At other settings #