Numerical Linear Algebra
where the algebra stops being the arithmetic
Reduction, and what a model is for

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 this site can predict from one number — 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.
38 min read 34 figures Identical algebra, different arithmeticWhere the algebra stops being the arithmeticMeasured, not asserted

Worth reading first: The bound that is known in advance · An index that is a pair · The road that squares the problem.

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.

Hankel singular values, and the half of them a product cannot seeσₖ divided by σ₁, for a 20-state model of McMillan degree 12. The lower curve is the square-root route — Cholesky-like factors of the two Gramians and one SVD of RᵀS — and it descends to the level at which the Gramians were computed. The upper curve eigendecomposes the product PQ, agrees for the first 6 values and then stops, flattening at 9.67·10⁻¹⁰. The dashed line is σ₁√u = 2.32·10⁻⁹, which is where a route that squares the conditioning must stop: the eigenvalues of PQ are σ², an absolute error of uσ₁² on them is a relative error of √u on σ. κ(PQ) is 1.8·10³⁶ against 1.34·10¹⁸. Every σ below the line is a term in the error bound, so this is not an academic loss.1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =6product floors at9.7·10⁻¹⁰σ₁√u2.3·10⁻⁹κ(P)κ(Q)1.8·10³⁶√ of it1.3·10¹⁸do not form the productthe σ below the line are the bound
Fig. 1 The σ by both routes. They agree, and then one of them stops.

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⁻¹⁰.

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.

Why the product route stops: two condition numbers and their productThe two Gramians of a 22-state model of McMillan degree 10, and what each route to the Hankel singular values is charged. κ(P) = 1.26·10¹⁸ and κ(Q) = 2.16·10¹⁸. The square-root route works with RᵀS, whose condition number is their geometric mean, 1.65·10¹⁸; the route that eigendecomposes PQ works at 2.72·10³⁶, which is past 1/u = 4.5·10¹⁵ — the point at which nothing small survives at all. The Gramians themselves are right: this one agrees with its closed form to 2.15·10⁻¹³. What is lost is lost in the last step, to a product nobody had to form.condition numbers, on a logarithmic scaleκ(P)1.26·10¹⁸κ(Q)2.16·10¹⁸√(κ(P)κ(Q)) — the SVD route1.65·10¹⁸κ(P)κ(Q) — the product route2.72·10³⁶1/u4.5·10¹⁵the mean, or the productand only one of them fits
Fig. 2 Two condition numbers, their mean and their product, and the line at 1/u.

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Σ_{k>r}σₖ 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Σ^{−1/2} and T = SVΣ^{−1/2}, where U and V are the singular vectors of RᵀS and Σ holds the σ being kept. The Σ^{−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.

At other settings

Hankel singular values, and the half of them a product cannot seeσₖ divided by σ₁, for a 10-state model of McMillan degree 6. The lower curve is the square-root route — Cholesky-like factors of the two Gramians and one SVD of RᵀS — and it descends to the level at which the Gramians were computed. The upper curve eigendecomposes the product PQ, agrees for the first 5 values and then stops, flattening at 1.45·10⁻⁹. The dashed line is σ₁√u = 2.39·10⁻⁹, which is where a route that squares the conditioning must stop: the eigenvalues of PQ are σ², an absolute error of uσ₁² on them is a relative error of √u on σ. κ(PQ) is 2.01·10³⁵ against 4.49·10¹⁷. Every σ below the line is a term in the error bound, so this is not an academic loss.1357910⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =5product floors at1.4·10⁻⁹σ₁√u2.4·10⁻⁹κ(P)κ(Q)2·10³⁵√ of it4.5·10¹⁷do not form the productthe σ below the line are the bound
Fig. 3 A ten-state model, where the σ run out before the floor is reached.
Hankel singular values, and the half of them a product cannot seeσₖ divided by σ₁, for a 14-state model of McMillan degree 10. The lower curve is the square-root route — Cholesky-like factors of the two Gramians and one SVD of RᵀS — and it descends to the level at which the Gramians were computed. The upper curve eigendecomposes the product PQ, agrees for the first 5 values and then stops, flattening at 9.19·10⁻⁹. The dashed line is σ₁√u = 2.23·10⁻⁹, which is where a route that squares the conditioning must stop: the eigenvalues of PQ are σ², an absolute error of uσ₁² on them is a relative error of √u on σ. κ(PQ) is 2.91·10³⁴ against 1.71·10¹⁷. Every σ below the line is a term in the error bound, so this is not an academic loss.13579111310⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =5product floors at9.2·10⁻⁹σ₁√u2.2·10⁻⁹κ(P)κ(Q)2.9·10³⁴√ of it1.7·10¹⁷do not form the productthe σ below the line are the bound
Fig. 4 Fourteen, where the two routes part at the fifth value.
Hankel singular values, and the half of them a product cannot seeσₖ divided by σ₁, for a 26-state model of McMillan degree 12. The lower curve is the square-root route — Cholesky-like factors of the two Gramians and one SVD of RᵀS — and it descends to the level at which the Gramians were computed. The upper curve eigendecomposes the product PQ, agrees for the first 6 values and then stops, flattening at 8.07·10⁻¹⁰. The dashed line is σ₁√u = 2.32·10⁻⁹, which is where a route that squares the conditioning must stop: the eigenvalues of PQ are σ², an absolute error of uσ₁² on them is a relative error of √u on σ. κ(PQ) is 1.58·10³⁸ against 1.26·10¹⁹. Every σ below the line is a term in the error bound, so this is not an academic loss.1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =6product floors at8.1·10⁻¹⁰σ₁√u2.3·10⁻⁹κ(P)κ(Q)1.6·10³⁸√ of it1.3·10¹⁹do not form the productthe σ below the line are the bound
Fig. 5 Twenty-six, and the floor has not moved, because it is about the arithmetic.
Hankel singular values, and the half of them a product cannot seeσₖ divided by σ₁, for a 34-state model of McMillan degree 12. The lower curve is the square-root route — Cholesky-like factors of the two Gramians and one SVD of RᵀS — and it descends to the level at which the Gramians were computed. The upper curve eigendecomposes the product PQ, agrees for the first 6 values and then stops, flattening at 8.06·10⁻¹⁰. The dashed line is σ₁√u = 2.32·10⁻⁹, which is where a route that squares the conditioning must stop: the eigenvalues of PQ are σ², an absolute error of uσ₁² on them is a relative error of √u on σ. κ(PQ) is 1.76·10³⁸ against 1.33·10¹⁹. Every σ below the line is a term in the error bound, so this is not an academic loss.1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =6product floors at8.1·10⁻¹⁰σ₁√u2.3·10⁻⁹κ(P)κ(Q)1.8·10³⁸√ of it1.3·10¹⁹do not form the productthe σ below the line are the bound
Fig. 6 Thirty-four, where there are more σ below the floor than above it.
Why the product route stops: two condition numbers and their productThe two Gramians of a 22-state model of McMillan degree 4, and what each route to the Hankel singular values is charged. κ(P) = 2.69·10¹⁸ and κ(Q) = 4.12·10¹⁸. The square-root route works with RᵀS, whose condition number is their geometric mean, 3.33·10¹⁸; the route that eigendecomposes PQ works at 1.11·10³⁷, which is past 1/u = 4.5·10¹⁵ — the point at which nothing small survives at all. The Gramians themselves are right: this one agrees with its closed form to 2.14·10⁻¹³. What is lost is lost in the last step, to a product nobody had to form.condition numbers, on a logarithmic scaleκ(P)2.69·10¹⁸κ(Q)4.12·10¹⁸√(κ(P)κ(Q)) — the SVD route3.33·10¹⁸κ(P)κ(Q) — the product route1.11·10³⁷1/u4.5·10¹⁵the mean, or the productand only one of them fits
Fig. 7 A model of degree four, whose Gramians are better conditioned and whose product still is not.
Why the product route stops: two condition numbers and their productThe two Gramians of a 22-state model of McMillan degree 14, and what each route to the Hankel singular values is charged. κ(P) = 2.76·10¹⁹ and κ(Q) = 1.3·10¹⁸. The square-root route works with RᵀS, whose condition number is their geometric mean, 5.99·10¹⁸; the route that eigendecomposes PQ works at 3.59·10³⁷, which is past 1/u = 4.5·10¹⁵ — the point at which nothing small survives at all. The Gramians themselves are right: this one agrees with its closed form to 2.14·10⁻¹³. What is lost is lost in the last step, to a product nobody had to form.condition numbers, on a logarithmic scaleκ(P)2.76·10¹⁹κ(Q)1.3·10¹⁸√(κ(P)κ(Q)) — the SVD route5.99·10¹⁸κ(P)κ(Q) — the product route3.59·10³⁷1/u4.5·10¹⁵the mean, or the productand only one of them fits
Fig. 8 And one of degree fourteen, where the decay is slower and the product route survives further.
Least squares by QR and by the normal equations in binary32Relative error of the computed coefficients against epsilon, on log axes. The QR route is a flat line near the bottom; the normal-equations route climbs and then stops, at the epsilon where the cross-product matrix becomes exactly singular.10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴
Fig. 9 The essay this one is shaped after, in the field it was written for.
QᵀQ from classical Gram–Schmidt and from Householder on the 10×10 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.A = H10 · κ = 1.6·10¹³ · both factorisations reconstruct A to 6.4·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.00000000000001.00000000000001.00000000000001.00000000000001.0000-0.001-0.001-0.001-0.001000001.000-0.053-0.042-0.035-0.0280000-0.001-0.0531.0001.0001.0001.0000000-0.001-0.0421.0001.0001.0001.0000000-0.001-0.0351.0001.0001.0001.0000000-0.001-0.0281.0001.0001.0001.000classical Gram–Schmidt1.00000000000001.00000000000001.00000000000001.00000000000001.00000000000001.00000000000001.00000000000001.00000000000001.00000000000001.000Householderclassical ‖QᵀQ − I‖3.5Householder ‖QᵀQ − I‖1.9·10⁻¹⁵largest off-diagonal 1 against 4·10⁻¹⁶length is not angle
Fig. 10 Two procedures that are the same algebra and not the same arithmetic.
Classical and modified Gram–Schmidt: the same subtraction, in a different orderTwo panels of pseudocode differing in one argument, with the resulting pairwise dot products of the computed columns listed beneath each.for each previous column i, subtract the projection of column j onto qᵢclassicalr[i][j] = qᵢ · a[j] ↑ the ORIGINAL columnv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 8 and 9: 1columns 9 and 10: 1columns 7 and 8: 1modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 2 and 10: 1.8·10⁻⁴columns 1 and 10: 3.2·10⁻⁵columns 4 and 10: 7.9·10⁻⁶The two R factors agree to 2.9·10⁻⁶ relative. The two Q factors do not.the 10×10 Hilbert matrixone word, eight orders
Fig. 11 And the measurement that tells them apart.
Along a sequence of exactly rank-two tensors: the step's condition number, and what a fit from a random start achievesEvery A_n on this sequence *is* a rank-two tensor and its two rank-one terms are written down, so nothing here is about existence. The rising curve is the condition number of the r × r system each alternating sweep solves, which on this sequence has a closed form in n whose asymptote is 2n² — the marks are measured and the dashed line is that closed form, agreeing to 1.1·10⁻¹¹. The lower marks are what a three-hundred-sweep fit from a random start returns: 2.1·10⁻⁴ at n = 2 rising to 0.0118 at n = 1024. Started at the answer instead, the same code stays within 9.2·10⁻⁹ of it at every n — not the rounding level, because the drift from an exact start is itself about κ times the unit roundoff, but nine orders below what a random start reaches. That is the control that says the failure is the conditioning and not the implementation. A tensor away from the boundary conditions its step at 6.03.110¹10²10³10⁻¹²10⁻⁸10⁻⁴110⁴ncondition number, and residual reachedmarks above: κ of the step · dashes: its closed formmiddle: a fit from a random startbelow: the same fit started at the answera decomposition that is ill-conditionedκ at n = 10242.1·10⁶its closed form2.1·10⁶cosine of the terms1from a random start0.012from the answer9.2·10⁻⁹the answer existsand cannot be found
Fig. 12 What a decomposition inherits from the matrix it is taken of.
The error a reduced model has, and the bound computed before it existed‖H − Hᵣ‖∞ for balanced truncation of a 20-state model, against 2Σ_{k>r}σₖ — a number available from two Lyapunov solves before any reduced model is formed — and against σᵣ₊₁, which is a lower bound no reduced model of that order can beat. The measured error and the upper bound are the same curve: the ratio is 1.0000 to 1.0000 across the sweep, spread 1.0000. The bound is attained rather than approached, because on a spectrum falling this fast the tail sum is its own first term and removing one state costs exactly twice it. The measured error sits at 2.03 times the lower bound, so the whole question of what order to truncate at is answered by a curve nobody has to compute a model to see.123456710⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹order r kepterror and bounds2Σσ, and the errorσᵣ₊₁a bound that is an equalityorders7bound ÷ error, worst1spread over the sweep1error ÷ σᵣ₊₁, worst2.1computed before the modeland attained by it
Fig. 13 The bound the σ below the floor are terms of.
Removing one state costs exactly twice the Hankel singular value it removed‖H − Hₙ₋₁‖∞ divided by σₙ, for 5 systems that share nothing but the shape of the question: residues of one sign and of both, a clustered pair, a spectrum spanning three decades. Every one comes back at 2.0000. That is where the bound for a many-state truncation comes from — it is this equality applied once per removed state and the terms added up — and it is why the sum can only be loose in how the removals interact, never in the single step. The measurement is the whole reason to trust a bound that was computed from two Lyapunov solves and no reduced model.‖H − Hₙ₋₁‖∞ ÷ σₙpositive residues2.0000mixed residues2.0000clustered poles2.0000four poles2.0000wide spectrum2.0000‖ΔH‖ 3.89·10⁻⁴ σ 1.94·10⁻⁴‖ΔH‖ 3.92·10⁻⁴ σ 1.96·10⁻⁴‖ΔH‖ 0.00502 σ 0.00251‖ΔH‖ 8.41·10⁻⁴ σ 4.21·10⁻⁴‖ΔH‖ 4.6·10⁻⁴ σ 2.3·10⁻⁴2four systems, four ratios, one number: 2.0000one state, one equalityand the sum is that, repeated
Fig. 14 Where that bound comes from, one state at a time.
The growth factor of a 10×10 elimination, against the condition number of the matrixThree curves against κ. Cholesky's growth factor is exactly 1 at every condition number drawn — the elimination never produces an entry larger than the matrix already had. Partial pivoting on the same matrices reaches the same growth, and reaches it by making up to 5 row interchanges, each of which destroys the symmetry that was the reason to use a symmetric factorisation. The bound the general theory allows, 2^9 = 512, is drawn above them both.10¹10³10⁵10⁷10⁹10¹¹110¹10²10³10⁴condition number of the matrixgrowth factorbound 2^9partial pivotingCholeskyno pivot to gain fromCholesky growth, every κ1Cholesky interchanges0partial pivoting, worst5the bound, 2^9512both eliminations reach the same growthand only one of them had to swap to get there
Fig. 15 The factorisation that is not the right square root here, and why.
How often Cholesky still calls a 10×10 matrix positive definite, against λ_min/λ_max in units of the format's own roundoffThree curves, one per precision, of the share of 24 seeded matrices on which the factorisation succeeds. Measured in units of each format's unit roundoff the three lie almost on top of one another, with the edge — the smallest ratio at which every seed succeeds — at 1.0, 1.0, 1.8 times u. The absolute thresholds are 6·10⁻⁸, 2.3·10⁻¹⁰ and 2·10⁻¹⁶: nine orders of magnitude apart, and the same number in the format's own units.10⁻³10⁻²10⁻¹110¹10²10³10⁴00.250.50.751λ_min / λ_max, in units of the format's own ushare of seeds that succeed24 bits32 bits53 bitsone threshold, three formats24-bit edge, absolute6·10⁻⁸32-bit edge, absolute2.3·10⁻¹⁰53-bit edge, absolute2·10⁻¹⁶in units of u, at 53 bits1.8a yes-or-no question with a precision in itand a coin flip three decades below the edge
Fig. 16 The edge a Cholesky of a numerically singular matrix falls off.
Every column's 1-norm on a 12×12 matrix built to stop Hager's walk one column short12 bars, one per column, each the sum of the absolute values in that column — the quantity the estimator maximises, whose largest is the matrix's 1-norm. The walk visits 1 of them and returns 12.00; column 2 has 1-norm 114.00 and is never visited, because its entries alternate in sign and cancel against every sign vector the walk stands on.the estimator maximises this quantity over the columns it visitscolumn 1 ‹visited›12column 2 ‹the answer›114column 311.4column 411.4column 511.4column 611.4column 711.4column 811.4column 911.4column 1011.4column 1111.4column 1211.4estimate 12.0a walk that stopped earlythe estimate returned12the true 1-norm114columns visited1products with the matrix5the walk's own stopping test firedand every column it could see was smaller
Fig. 17 An estimate that is right until it is not, one field over.
What Hager's estimator returns, as a share of the truth, against the size of the matrix built to defeat itThe estimate over the true 1-norm against n, both axes logarithmic. It falls from 0.526 at n = 4 to 0.0658 at n = 32, along the line 1/t where t is the construction's own multiplier — which is bounded only by the size of the matrix, so the ratio has no floor. A line at one marks a correct estimate.10¹10⁻²10⁻¹1size of the matrixestimate ÷ true 1-norma correct estimatewhat it returns1 / tno floorratio at n = 40.53ratio at n = 320.066products, either size5the estimate is always a lower boundwhich is the direction that flatters the matrix
Fig. 18 And the case built to fool it.
What each way of eliminating a constraint inherits, over five decades of κ(A)The same system solved twice, at 8 unknowns and 3 constraints with κ(H) = 1. The range-space method forms S = AH⁻¹Aᵀ and inherits κ(S), which rises from 1 to 10·10⁹ — the square of κ(A), for the reason the normal equations square it. The null-space method solves with the reduced Hessian ZᵀHZ, whose condition number is 1 at the start of the sweep and 1 at the end: it does not contain κ(A) at all. The two forward errors, measured against a solution computed in BigInt rationals, follow their own condition numbers: 7.341·10⁻⁷ against 2.025·10⁻¹² at the far end. Both methods are correct and one of them is usable.01234510⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹log₁₀ κ(A)condition number, and relative errorκ(S)κ(ZᵀHZ)range-space errornull-space erroragainst a BigInt answerκ(S) at κ(A) = 10⁵10·10⁹κ(ZᵀHZ), all stops1range-space forward error7.3·10⁻⁷null-space forward error2·10⁻¹²both are the same algebraand only one squares
Fig. 19 Two eliminations of one matrix, differing in what they form.
Two condition numbers of one 8×8 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 9.83 to 1.9·10¹⁰ while the componentwise condition number stays at 6.98 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.4·10¹⁰ and 1.2·10¹⁰: a matrix whose sensitivity no scaling repairs.0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 20 The free decision that governs the visible one, in another field.
What two lines of scaling are worth: the same quadratic in nine systems of unitsThe forward error against the closed form for an overdamped chain of 8 masses, before and after Fan–Lin–Van Dooren scaling. The unscaled curve runs 7.49·10⁻¹⁴ to 0.001258 — every digit gone by the far end — and the scaled one runs 1.26·10⁻¹³ to 8.23·10⁻¹⁴, flat to within a factor of 2.01. The scaling is γ = √(‖K‖/‖M‖) and δ = 2/(‖K‖ + γ‖C‖), computed from three norms and nothing else, and the map back is λ = γμ with no rounding in the statement. Flatness is the half that matters: after scaling every stop of the sweep IS the same problem, so there is nothing left for the change of units to do.0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.0013scaled, worst1.7·10⁻¹³orders recovered10scaled coefficient spread4.5the answer was never the problemthe units were
Fig. 21 And the clearest instance of it on this site.
How much a perturbation of the right-hand side is amplified, κ = 10⁴The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10000worst found 76206×6, 200 directionsmedian reaches 0.29 of κ
Fig. 22 The amplifier, and what squaring it costs.
Forward error of a 6×6 Hilbert solve at eight precisionsA bar for each significand width from 12 to 53 bits showing the relative error in the computed solution, with the condition number times the unit roundoff marked as a prediction.κ = 1.5·10⁷ · the exact answer is (1, 2, …, 6)12 bits3.316 bits0.8620 bits4.324 bits0.09330 bits7.6·10⁻⁴36 bits6.2·10⁻⁶43 bits2·10⁻⁷53 bits4.5·10⁻¹¹dashed: κ · unit roundoffone matrix, eight arithmeticsmeasured against a known answer
Fig. 23 The knob this site is built around, for comparison with a floor that cannot be moved by it.
Relative error of each singular value of a bidiagonal matrix graded over 11 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 1.934·10⁻¹¹. One-sided Jacobi holds every value to 7.1·10⁻¹⁶ relative and the zero-shift bidiagonal sweep to 4.41·10⁻¹⁶; the eigenvalues of BᵀB are wrong by 14.5 — not a rounding error but the whole number — and the damage is entirely at the small end.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly1.9·10⁻¹¹worst, one-sided Jacobi7.1·10⁻¹⁶worst, zero-shift QR4.4·10⁻¹⁶worst, eigenvalues of BᵀB15a relative error is a ratioand the denominator is the answer
Fig. 24 A spectrum whose small end is what a squared condition number destroys.
Singular values of a rank-4 matrix with noise of relative size 10⁻⁴Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 820decisivetrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 25 Reading a rank off a curve that has stopped descending.
Eigenvalues of a singular 6×6 pencil after a perturbation of 10^-11, over 8 seedsA and B share a null vector, so det(A − λB) is identically zero and every λ is an eigenvalue: the question has no answer. Perturbing by 10^-11 makes the pencil regular, and a solver then returns 6 eigenvalues with residuals no larger than 1.89·10⁻¹² — small by any standard a caller would apply, and 2.2·10¹⁶ times smaller than the disagreement between seeds. Each row here is one seed. No two rows agree, the values are spread over 40800, and nothing in the output distinguishes this from a well-posed problem. The only thing that says so is the exact characteristic polynomial, which has no coefficients at all.-40-27-14-11225380123456789computed eigenvalueseedevery mark has a residual below 10⁻⁸a small residual, and no answercoefficients of det(A − λB)0worst residual over all seeds1.9·10⁻¹²spread of the answers4.1·10⁴seeds drawn8the residual is small at every markand none of the marks means anything
Fig. 26 Small singular values under perturbation, which is what the floor is.
Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 30-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 2.48·10⁻⁸ in 10 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.86 and 3834, whose ratio is 388.8 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates30κ of the spectrum389λ11 ÷ λ₁2.5·10⁻⁸the bound there5.2·10⁻⁴the cliff everything rests onand the reason for it
Fig. 27 Why the σ fall at all, and how far they would fall if they could be seen.
The two 3×3 matrices of a Sylvester equation, and the 9×9 matrix it meansA and B are 3×3 with 5 nonzeros each. The coefficient matrix of the map X ↦ AX + XB is 9×9 with 21 nonzeros, drawn at the same scale so the ratio is visible rather than quoted. Bartels and Stewart's algorithm never forms it: two Schur reductions and a back-substitution give the same X to 5.5·10⁻¹⁷, and that X satisfies AX + XB = C to 8.8·10⁻¹⁷.A3×3B3×3I ⊗ A + Bᵀ ⊗ I9×9one equation, two objectsentries in A and B18entries in the coefficient matrix81two routes, relative gap5.5·10⁻¹⁷‖AX + XB − C‖/‖C‖8.8·10⁻¹⁷the small squares are the problemand the large one is the notation
Fig. 28 The equation both Gramians solve.
sep(A, B) and the smallest |λᵢ(A) + μⱼ(B)| as the departure from normality grows, n = 6A and B are upper triangular with diagonals 1, 2, …, 6 and μ above them, so every eigenvalue sum is at least 2 at every point on the axis and the flat line is exact rather than nearly flat. sep — the smallest ‖AX + XB‖_F over ‖X‖_F = 1, and the reciprocal of the amplification a perturbation of C receives — starts equal to the gap at μ = 0, where the matrices are normal, and falls to 9.47·10⁻⁴ at μ = 8. The number a reader is invited to consult is the one that does not move.01234567810⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹μ, the entry above the diagonalsep, and the eigenvalue gapmin |λᵢ + μⱼ|sep(A, B)the spectra never moveeigenvalue gap, throughout2sep at μ = 02sep at μ = 89.5·10⁻⁴amplification there1056solvability is the eigenvaluesand conditioning is not
Fig. 29 And the number that decides how accurately it can be solved.
The transfer function of a 24-state model, by a solve and by a formula|H(iω)| for the second-difference operator on 24 points with one actuator and one sensor. The curve is C(iωI − A)⁻¹B, one factorisation of a complex 24 × 24 matrix per frequency, and the crosses are Σ(cᵀvₖ)(vₖᵀb)/(iω − λₖ) with the eigenvalues and eigenvectors written down in closed form. The two share the coefficient matrices and no code. Below ω = 100 they agree to 1.05·10⁻¹⁴ relatively; above it the modal sum is computing a number far below its own largest term, because the residues sum to cᵀb and the sensor is not at the actuator, so the relative comparison stops meaning what it looks like. Measured against the size of the terms being added, the two agree to 3.85·10⁻¹⁵ at every frequency. A state dimension of 24 and a function of one variable: the gap between those is what the field is about.10⁻¹110¹10²10³10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹frequency ω|H(iω)|one function, two routesstates24relative, ω ≤ 10010⁻¹⁴against the terms3.9·10⁻¹⁵‖H‖∞0.0924 statesone function of one variable
Fig. 30 The model all of this is about.
Two rank decisions for one integer: the McMillan degree of a 40-state modelSingular values, each divided by the largest of its own set. The Loewner matrix is built from 32 samples of H on the imaginary axis and never touches A: its 6th singular value stands 9.46·10⁴ above the next, so the rank decision is not a close call. The Hankel matrix of Markov parameters CAᵏB is built from A, B and C — strictly more information — and has no cliff, because ‖A‖ is (n+1)² and its entries span 3·10³⁷. It returns 2 where the answer is 6. The definition is not the algorithm, and the route with less information is the one that works.1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular value ÷ the largestdegree 6from samples of Hfrom A, B and Cone integer, two routesstates40Loewner rank6its gap9.5·10⁴Markov rank2its dynamic range3·10³⁷the rank of a divided differenceand the rank that cannot be seen
Fig. 31 Its degree, which bounds how many σ there are to lose.
Three iterations to the orthogonal polar factor, κ = 10^1Newton's iteration, X ← (X + X⁻ᵀ)/2, halves its error per step while it is far away and only becomes quadratic near the end: it is at 1.9·10⁻⁶ after six steps. Higham's scaling costs two norms and no extra factorisation and reaches thirteen digits in 5. Newton–Schulz, X ← X(3I − XᵀX)/2, uses no inverse at all — two matrix products a step and nothing that reads an entry — and needs 11 steps to get to the same place.1471013161910⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹steprelative error in the orthogonal factorNewtonNewton, scaledNewton–Schulzone fixed point, three costsscaled Newton, steps5Newton–Schulz, steps11Newton at step 61.9·10⁻⁶scaled Newton at step 67.5·10⁻¹⁶a Newton step needs an inverseand a Schulz step needs two products
Fig. 32 Another square root of a matrix, computed without forming what it is a root of.
‖A − Q‖_F for four orthogonal matrices, on an 8×8 matrix with κ = 10The polar factor is 1.8554 from A. QR with its column signs fixed is 2.1265 — 15 per cent further. QR as Householder returns it, with 7 of 8 columns negated, is 3.8226, which is further than the best of two hundred orthogonal matrices drawn at random. The strip beneath the bars is those two hundred draws, whose best is 2.6928; the polar factor is to the left of all of them, which is the minimisation being checked rather than assumed.how far is it from A to an orthogonal matrix?smaller is nearer · the polar factor minimises this in every unitarily invariant normpolar factor U1.8554QR, signs fixed2.1265QR as returned3.8226200 drawn at randomκ = 10polar factor1.9QR, signs fixed2.1QR as returned3.8best of 200 random2.7‖A − QR‖ is the same either wayand ‖A − Q‖ is not
Fig. 33 And the object that iteration converges to.
Iterative refinement from a 11-bit factorisation, κ = 10⁴A semi-logarithmic plot of forward error against refinement step. One curve falls steeply to the level of a double-precision solve; the other is nearly flat.012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in11-bitresidual indoubleone argument apartκ·u of the factorisation4.9double residual, final0.0049same-precision, final1.430×30, κ = 10⁴, same factors in both runsidentical cost
Fig. 34 Buying digits back after losing them, which is not available here.

Where this goes next

123456710⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹order r kepterror and bounds2Σσ, and the errorσᵣ₊₁a bound that is an equalityorders7bound ÷ error, worst1spread over the sweep1error ÷ σᵣ₊₁, worst2.1computed before the modeland attained by it The bound that is known in advance Almost every error on this site 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. 10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴ The road that squares the problem The normal equations are the first method every course teaches and the method no library uses. Forming AᵀA squares the condition number, and below ε = √u it does not degrade — it produces a matrix that is exactly singular, from data that was perfectly usable. -12-10-8-6-4-2010⁻¹⁷10⁻¹³10⁻⁹10⁻⁵10⁻¹10³10⁷10¹¹10¹⁵log₁₀ μcondition number, and relative errorκ₂, augmentedcomponentwise, condensedcomponentwise, augmentederror, augmentedone matrix, two numbersκ₂ at μ = 10⁻¹²3·10¹³componentwise, same matrix13their ratio2.3·10¹²measured error9.4·10⁻¹⁶every library prints the top lineand the error obeys the third Two condition numbers of one matrix κ₂ is a worst case over perturbations of a given norm, and a normwise perturbation may put its whole budget on the smallest entry. The componentwise number is a worst case over perturbations proportional to the entries, which is what a backward-stable factorisation actually makes. On one matrix they are 3·10¹³ and 13.3, and the error obeys the second.

What links here

Computed from the collection, not written here: the essays that point at this one.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Balanced truncationCondition numberGramianHankel singular valuesLyapunov equationNormal equationsSingular valuesSquare root methodUnit roundoff