Numerical Linear Algebra
where the algebra stops being the arithmetic
Least squares, and the road not to take

A constraint is a weight at infinity

Stack an equality constraint on top of a least-squares problem with a large weight and the answer approaches the constrained one like 1/τ². The limit is takeable to any accuracy — and how far it can be taken is a property of the solver, not of the problem. One of them stops at the square root of the precision, and one of them does not stop.
38 min read 29 figures Identical algebra, different arithmeticExact ground truthWhere the algebra stops being the arithmetic

Worth reading first: The projection and the right angle · The road that squares the problem · The zero that is not a missing entry.

Least squares with an equality constraint —

minimise ‖Ax − b‖ subject to Bx = d

— is what a data fit looks like when part of the answer is known: a curve required to pass through a point, a set of coefficients required to sum to one, a flow required to balance at a node. Its optimality conditions are a saddle-point system with AᵀA in the (1, 1) block, so everything the constraint field says applies to it.

It also has a property that nothing else in this collection has. A constraint is the limit of a weight, and the limit can be watched.

Minimise ‖Ax − b‖ subject to Bx = d, solved as a weighted least-squares problem, three waysStack the constraint on top of the objective with a weight τ and solve the ordinary least-squares problem that results. In exact arithmetic the answer approaches the constrained one like 1/τ² — measured here as exactly four orders of error per two decades of τ, against a solution computed in BigInt rationals from the problem's own optimality conditions. What stops the limit is the solver and not the problem. The normal equations on the weighted problem form entries of size τ², so the constraint block is lost once τ² passes 1/u: the last weight at which they are within an order of the right answer is 10⁴ here, against 1/√u = 6.71·10⁷, and by τ = 10¹⁴ they are wrong by 14.3. Householder QR has no such ceiling and is at 4.8·10⁻¹⁵ at the same weight. Classical Gram–Schmidt is worse than either, at 4.9·10¹⁰.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷10¹⁰log₁₀ τ — the weight on the constraintrelative error against the exact answerτ = 1/√uGram–Schmidtnormal equationsHouseholder QRthe ceiling is the method'sHouseholder, τ = 10¹⁴4.8·10⁻¹⁵normal equations14Gram–Schmidt4.9·10¹⁰1/√u6.7·10⁷a constraint is a weight at infinityand the solver decides how far infinity is
Fig. 1 Three solvers on one weighted problem. The slope they share is the algebra; the three places they stop are the arithmetic.

The construction, and why anybody would use it

Stack the constraint on top of the objective with a weight τ and solve the ordinary unconstrained problem

minimise ‖ [τB; A]x − [τd; b] ‖

Nothing about it is constrained. It is a least-squares problem of the sort the field’s first essay is about, with more rows than before, and any least-squares routine will take it.

As τ grows the constraint is enforced more and more nearly, and the solution approaches the constrained one. Measured against a solution computed in BigInt rationals from the problem’s own optimality conditions, the error is 4.69·10⁻⁵ at τ = 10², 4.69·10⁻⁹ at 10⁴ and 4.75·10⁻¹³ at 10⁶ — four orders per two decades of τ, which is 1/τ² exactly, at every stop.

The reason to do it at all is that it needs no special routine. A code that already has a least-squares solver acquires a constrained one for the price of stacking two arrays, and that is why the method is in use rather than a curiosity.

Minimise ‖Ax − b‖ subject to Bx = d, solved as a weighted least-squares problem, three waysStack the constraint on top of the objective with a weight τ and solve the ordinary least-squares problem that results. In exact arithmetic the answer approaches the constrained one like 1/τ² — measured here as exactly four orders of error per two decades of τ, against a solution computed in BigInt rationals from the problem's own optimality conditions. What stops the limit is the solver and not the problem. The normal equations on the weighted problem form entries of size τ², so the constraint block is lost once τ² passes 1/u: the last weight at which they are within an order of the right answer is 10⁴ here, against 1/√u = 6.71·10⁷, and by τ = 10¹⁴ they are wrong by 35.3. Householder QR has no such ceiling and is at 1.35·10⁻¹⁵ at the same weight. Classical Gram–Schmidt is worse than either, at 0.74.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷10¹⁰log₁₀ τ — the weight on the constraintrelative error against the exact answerτ = 1/√uGram–Schmidtnormal equationsHouseholder QRthe ceiling is the method'sHouseholder, τ = 10¹⁴1.4·10⁻¹⁵normal equations35Gram–Schmidt0.741/√u6.7·10⁷a constraint is a weight at infinityand the solver decides how far infinity is
Fig. 2 With a single constraint, where the same slope appears and the ceilings do not move.

Why the limit is quadratic and not linear

The exponent is worth deriving, because “converges” is not a measurement and “τ⁻²” is.

At finite τ the weighted problem’s solution x(τ) satisfies the normal equations (τ²BᵀB + AᵀA)x = τ²Bᵀd + Aᵀb. Write x(τ) = x* + e where x* is the constrained solution, and expand. The leading term of the residual comes from the fact that x* satisfies Bx* = d exactly but does not make Aᵀ(Ax* − b) zero — it makes it equal to −Bᵀλ for the multiplier λ. Balancing that against τ²BᵀBe gives e = O(τ⁻²) in the constrained directions and, because the objective already determines the rest, O(τ⁻²) overall.

The feasibility itself converges at the same rate: ‖Bx(τ) − d‖/‖d‖ is measured at 2.94·10⁻⁵, 2.94·10⁻⁹ and 2.93·10⁻¹³ at τ = 10², 10⁴ and 10⁶, which is 1/τ² to three digits.

That the two quantities share an exponent is not automatic — a method can approach feasibility faster than it approaches optimality or the other way round — and it is worth having measured, because the randomised field’s version of this page finds a case where the exponent is halved and the two come apart.

Minimise ‖Ax − b‖ subject to Bx = d, solved as a weighted least-squares problem, three waysStack the constraint on top of the objective with a weight τ and solve the ordinary least-squares problem that results. In exact arithmetic the answer approaches the constrained one like 1/τ² — measured here as exactly four orders of error per two decades of τ, against a solution computed in BigInt rationals from the problem's own optimality conditions. What stops the limit is the solver and not the problem. The normal equations on the weighted problem form entries of size τ², so the constraint block is lost once τ² passes 1/u: the last weight at which they are within an order of the right answer is 10⁴ here, against 1/√u = 6.71·10⁷, and by τ = 10¹⁴ they are wrong by 2.62. Householder QR has no such ceiling and is at 6.55·10⁻¹⁶ at the same weight. Classical Gram–Schmidt is worse than either, at 1.24·10¹¹.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷10¹⁰log₁₀ τ — the weight on the constraintrelative error against the exact answerτ = 1/√uGram–Schmidtnormal equationsHouseholder QRthe ceiling is the method'sHouseholder, τ = 10¹⁴6.6·10⁻¹⁶normal equations2.6Gram–Schmidt1.2·10¹¹1/√u6.7·10⁷a constraint is a weight at infinityand the solver decides how far infinity is
Fig. 3 With more data rows, where the constant in front of τ⁻² moves and the exponent does not.

Where the limit stops, and it is not where the problem says

In exact arithmetic the limit is takeable to any accuracy: 1/τ² is a statement with no rounding in it. What stops it is the solver, and the three solvers stop in three different places.

The normal equations on the weighted problem form (τB)ᵀ(τB) + AᵀA, whose constraint part has entries of size τ². Once τ² reaches 1/u the objective’s contribution is lost in the rounding of the constraint’s, and the matrix becomes numerically the constraint block alone — which is singular, because B has fewer rows than columns. The last weight at which they are within an order of the right answer is 10⁴, against 1/√u = 6.7·10⁷; by τ = 10¹⁴ the answer is wrong by a factor of fourteen.

Householder QR has no such ceiling. It never forms a Gram matrix, so nothing is ever squared, and its error falls with the 1/τ² slope until it reaches the rounding level at τ = 10⁶ and stays there: 4.8·10⁻¹⁵ at τ = 10¹⁴.

Classical Gram–Schmidt is worse than either, for the reason the essay that compared the two Gram–Schmidts gives: its loss of orthogonality grows like κ², and here κ grows like τ.

So the ceiling belongs to the method. The same problem, the same weights, and one route is at 10⁻¹⁴ where another has no digits at all.

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. 4 The same square, in the field’s own essay about it.

The square root of the precision, again

τ ≈ 1/√u is a threshold this site has met three times now and it is the same threshold each time.

The normal equations lose the solution when κ(A) passes 1/√u, because they solve with κ(A)². The definiteness test fails at a distance from singularity of a few ulps for a related reason. And here, a weight passes 1/√u and its square passes 1/u, at which point adding an O(1) quantity to a τ² one changes nothing.

The rule of thumb it produces is worth stating because it is often quoted without the “for the normal equations” attached: the method of weighting can only impose a constraint to about half the machine’s digits. That is true of a solver that forms a Gram matrix and false of one that does not, and the figure separates the two on one pair of axes.

‖QᵀQ − I‖ of the implied Q, against the condition numberLoss of orthogonality against κ, both axes logarithmic, for three factorisations of the same 256×8 matrix. The column sweep and the reduction tree run from 2.7·10⁻¹⁴ to 1.5·10⁻¹⁰. Cholesky QR — the one with a single reduction — runs from 2.3·10⁻¹³ to 1.8·10⁻⁵, a fitted slope of 1.97 against the others' near-flat one.10²10³10⁴10⁵10⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶condition number κ‖QᵀQ − I‖Cholesky QRsweeptreethe price of one roundfitted slope, Cholesky QR2‖QᵀQ − I‖ at κ = 10·10⁵1.8·10⁻⁵the sweep's, at the same κ1.5·10⁻¹⁰one reduction instead of nand the condition number squared
Fig. 5 The same threshold in the orthogonality field, where the quantity crossing it is ‖QᵀQ − I‖.

And the order of the rows is worth eleven orders

Householder QR on the stacked matrix has no ceiling, and it has a condition attached that is easier to violate than to notice: the heavily weighted rows have to come first.

Same matrix, same right-hand side, same reflectors, same operation count, and the only difference is whether the constraint rows are above the data rows or below them:

τ constraint first constraint last 10⁴ 4.69·10⁻⁹ 4.69·10⁻⁹ 10⁸ 1.12·10⁻¹⁴ 5.19·10⁻⁹ 10¹⁴ 4.80·10⁻¹⁵ 2.08·10⁻³

Up to τ = 10⁴ they agree. At 10¹⁴ they are a factor of 4.3·10¹¹ apart, and the wrong ordering has lost ten digits.

The mechanism is in the first reflector. A Householder reflector is built from a whole column and annihilates all of it below the first entry; a column whose large entries are at the bottom has them added to the small ones on the way to computing the column’s norm, and the small ones are lost before they are ever used. With the large rows first, the first reflector eliminates them immediately and the small rows enter a problem that has already been scaled.

This is Powell and Reid’s observation, and it is the sharpest instance on this site of the theme it belongs to: identical algebra, identical arithmetic count, different sequencing, ten digits.

The same Householder QR on the same rows, in two ordersOne weighted least-squares problem, one factorisation, and the only difference is whether the 2 heavily weighted constraint rows are placed above the 12 data rows or below them. The algebra is identical, the operation count is identical, the reflectors are the same reflectors. Up to τ = 10⁴ the two answers agree. At τ = 10¹⁴ they are 4.8·10⁻¹⁵ and 0.00208 — a factor of 4.324·10¹¹, from the sequence of the rows and nothing else. The reason is that the first reflector is built from the first column, and a column whose large entries are at the bottom has them annihilated last, after the small ones have already been added to them and lost.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ τrelative error against the exact answerconstraint rows lastconstraint rows firstthe same arithmetic, resequencedconstraint first, τ = 10¹⁴4.8·10⁻¹⁵constraint last, same τ0.0021the ratio4.3·10¹¹multiplications, either order0same reflectors, same countand ten digits between them
Fig. 6 The two orderings, on one problem, with one solver.

What a row ordering is not

Three readings that survive the paragraph above and should not.

It is not pivoting. A least-squares QR with column pivoting reorders the unknowns and is about rank; this reorders the equations and is about scale. A code can do both, and doing the second does not require the first.

It is not equilibration. Scaling the rows to have equal norm would destroy the weighting, which is the whole point of the construction: the constraint rows are supposed to be large. Row scaling and row ordering are different operations and only the second is available here.

And it is not a property of Householder QR being fragile. The measurement shows the opposite: with the right ordering it is at the rounding level at a weight of 10¹⁴, which is a weight at which the alternative methods have nothing. The ordering is a condition on using a good method well.

The backward error of a solve along a sequence, with the pivot order chosen fresh, kept, and kept after equilibrationA 10×10 conflict grid, 100 unknowns, with three of its rows scaled down through 9 decades as the sequence runs — a row that was ordinary becoming a row that is small in the matrix's own units. The sparsity pattern is identical at every member, so a symbolic phase computed once stays valid throughout. Choosing a fresh order each time holds the backward error at the working precision. Keeping the first member's order costs 2 replaced pivots and a backward error of 4.82·10⁻⁹. Keeping the same order after dividing each row by its largest entry costs nothing at all: 6.55·10⁻¹⁷, with no pivot replaced anywhere in the run. The ringed points on the upper curve are the members at which another pivot fell below the floor and was replaced.02468101210⁻¹⁸10⁻¹⁶10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸member of the sequencebackward error of the solvefresh order · equilibrated kept orderringed: another pivot replaced hereone order, two sets of unitsfresh order, last member9.5·10⁻¹⁷kept order, last member4.8·10⁻⁹equilibrated, last member6.5·10⁻¹⁷pivots replaced, kept2pivots replaced, equilibrated0the perturbation floor √u1.1·10⁻⁸a reused pivot order is safe or notdepending on what the rows are measured in
Fig. 7 The other operation on rows, from the elimination field, and what it does to a pivot choice.

Which raises the question of why not solve the saddle-point system

The weighted problem is an approximation and the optimality conditions are exact, so a reader who has read the constraint field will ask why anybody weights at all.

Three answers, and none of them is that weighting is better.

The routine is already there, and adding rows to an array is a smaller change to a program than adding a solver. That is not a numerical argument and it is the reason the method is common.

The weighted problem is unconstrained, so it can be handed to an iterative least-squares method — LSQR, or a sketched solve — where the constrained version would need a projection at every step.

And the weight is a knob. Some problems want the constraint enforced nearly rather than exactly, because it came from a measurement rather than from a law, and then τ is a modelling parameter with a meaning rather than a numerical device. That is the case where the regularisation field’s whole question reappears: how much of the data to believe, expressed as a weight.

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) = 100. The range-space method forms S = AH⁻¹Aᵀ and inherits κ(S), which rises from 12.99 to 3.981·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 21.13 at the start of the sweep and 21.13 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: 5.314·10⁻⁶ against 5.788·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⁵4·10¹⁰κ(ZᵀHZ), all stops21range-space forward error5.3·10⁻⁶null-space forward error5.8·10⁻¹²both are the same algebraand only one squares
Fig. 8 The exact route, priced in the constraint field.

The other way the same trade appears

There is a mirror image of this page in the constraint field, and putting the two side by side is worth a paragraph because they look like opposite advice.

The sixth essay there takes an exact constrained system and perturbs it — adds δ to one block and −γ to the other — in order to make the factorisation well behaved, and then measures how much the perturbation costs. This page takes an exact constrained system and approximates it by a weight, and measures how much the approximation costs.

They are the same move in opposite directions. There, the constraint is exact and is loosened to buy a factorisation; here, the constraint is a weight and is tightened to buy exactness. And both end with the same number attached: the error is proportional to the perturbation, the constant is computable, and the practical range is bounded at both ends by the arithmetic.

The difference in the two answers is instructive. There the perturbation is removable, by refining against the unperturbed matrix. Here it is not: refining against the constrained problem would mean having a constrained solver, which is what the whole construction was avoiding. So the weighting method’s error has to be lived with, and how large it is at the largest usable τ is the number that decides whether the method is good enough.

What the regularisation costs, and what 0 steps of refinement take backSolving [[H + δI, Aᵀ], [A, −δI]] instead of K gives the exact answer to a different problem, so its error is proportional to δ: measured at 1489·δ across six decades, which is a slope of one and not a trend. Refining against the unregularised matrix — the residual formed with K and the correction solved with the regularised factorisation — removes that term entirely, because the perturbation was never in the residual. It works while δ is below σ_min(K) = 6.797·10⁻⁴, marked on the axis, and stops working above it: the iteration's contraction factor is δ/σ_min and a fixed point needs that under one. So the trade-off curve every regularisation essay on this site has drawn — a term falling in δ against a term rising in it — has, here, a left branch that can simply be removed.-14-12-10-8-6-4-210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ δrelative error against the exact answerδ = σ_min(K)no refinement0 stepsa left branch that can be removederror ÷ δ, unrefined1489σ_min(K)6.8·10⁻⁴refined at δ = 10⁻⁶0.0014refined at δ = 10⁻²0.93the perturbation is known exactlybecause the code chose it
Fig. 9 The mirror image, from the constraint field: a perturbation of known size, and an error proportional to it.

What the exact answer is, and how it is available

The comparison needs a constrained solution that was not computed by any of the three methods, and this page uses the optimality conditions in exact arithmetic.

Form AᵀA and Aᵀb in floating point, assemble the saddle-point system [[AᵀA, Bᵀ], [B, 0]] with the stored doubles, convert every entry to the dyadic rational it already is, and solve in BigInt. That gives the exact solution of the problem the machine is holding, which is the right reference: a comparison against the problem somebody meant would be measuring the rounding of AᵀA rather than the behaviour of the three solvers.

There is a pleasing circularity in it. The reference is computed by forming a Gram matrix and solving a saddle-point system — the two moves this page is about the numerical cost of — and it is exact because BigInt does not round. The moves are not wrong; their arithmetic is.

The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 10 The discipline in its older form, on a matrix whose entries were rationals before storage.

A weight is not always a device

The construction has been presented as a numerical trick, and in a good many problems it is the model rather than a trick, which changes what the measurement is for.

A constraint that came from a conservation law is exact and should be imposed exactly: mass does not nearly balance. A constraint that came from a calibration measurement is not exact, and imposing it exactly means asserting that one observation is infinitely more reliable than the others. In the second case τ is a statement about relative confidence and the right value of it is finite, small, and chosen from what the measurements are worth — which is the whole subject of another field here.

The numerical content of this page then reads differently. At a modelling τ of, say, 10³, none of the three ceilings is anywhere near: all three solvers agree, the ordering does not matter, and the 1/τ² statement is about a bias that was deliberately introduced and is not an error at all. The measurements above are about the regime where τ is being used to approximate infinity, and that regime is entered on purpose by somebody who wanted an equality and had a least-squares routine.

Separating the two is the difference between reading the figure as a warning and reading it as a specification.

The L-curve, and where four rules put λThe norm of the solution against the norm of its residual, on logarithmic axes, as λ sweeps eight decades. The curve has a corner: to the left of it the noise is being amplified and to the right the signal is being thrown away. Four points are marked — the three rules that use only the data, and the oracle, which requires the exact answer and is not a method.10⁻²10⁻¹110¹10²10³10⁴‖Ax − b‖‖x‖the oraclediscrepancyL-curvegeneralisedscored against a truth none hasoracle, relative error0.11discrepancy principle, as a multiple1.1L-curve corner, as a multiple2.3generalised cross-validation, as a multiple1the oracle needs the exact answer and is not a methodit is the reference the others are scored on
Fig. 11 Choosing a weight from the data rather than from the arithmetic, in the field that owns the question.

The refusals

The claim the page has to close is the one that reads as a fact about norms: the order of the rows of a least-squares problem does not affect its answer, since a permutation of the rows is a permutation of the residual and the norm is the same. It is exactly true in the algebra. The assertion is fed the pair at τ = 10¹⁴ and required to reject it.

The second refusal covers the reading that a large enough weight enforces the constraint whatever the solver: fed the normal equations at τ = 10¹², where the answer is off by a factor of two, and required to fail. The weight is exact and the arithmetic is not, and past a point more weight buys less.

The third is a control on the reference itself: the exact answer is fed the claim that it misses its own constraint, and required to refuse. An equality constraint is an equality, and a page whose whole subject is how nearly a constraint is satisfied has to know that its yardstick satisfies it exactly.

The same Householder QR on the same rows, in two ordersOne weighted least-squares problem, one factorisation, and the only difference is whether the 2 heavily weighted constraint rows are placed above the 30 data rows or below them. The algebra is identical, the operation count is identical, the reflectors are the same reflectors. Up to τ = 10⁴ the two answers agree. At τ = 10¹⁴ they are 8.97·10⁻¹⁶ and 0.00171 — a factor of 1.902·10¹², from the sequence of the rows and nothing else. The reason is that the first reflector is built from the first column, and a column whose large entries are at the bottom has them annihilated last, after the small ones have already been added to them and lost.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ τrelative error against the exact answerconstraint rows lastconstraint rows firstthe same arithmetic, resequencedconstraint first, τ = 10¹⁴9·10⁻¹⁶constraint last, same τ0.0017the ratio1.9·10¹²multiplications, either order0same reflectors, same countand ten digits between them
Fig. 12 With thirty data rows, where the gap between the two orderings is wider — more small rows for the large ones to be added to.

What to do with all of this

The page has three separate findings and they combine into one rule, which is worth stating plainly because the individual measurements do not.

If a constraint has to be imposed exactly, solve the saddle-point system. The weighting method cannot reach exact — its error is 1/τ² for a finite τ, and τ is bounded by the arithmetic — so “exactly” is not on offer.

If a constraint has to be imposed nearly, weight it, and then two things follow. Use a solver that does not form a Gram matrix, because the one that does stops at half the machine’s digits and stops without saying so. And put the weighted rows first, because a Householder QR that sees them last loses ten of the digits it would otherwise keep.

Neither of the second pair is expensive and neither is discoverable from the answer: the residual is small either way, the factorisation completes either way, and the only quantity that would say otherwise is a comparison against something the code does not have. That is the recurring shape of this whole collection, and it is why the measurement is the essay.

Minimise ‖Ax − b‖ subject to Bx = d, solved as a weighted least-squares problem, three waysStack the constraint on top of the objective with a weight τ and solve the ordinary least-squares problem that results. In exact arithmetic the answer approaches the constrained one like 1/τ² — measured here as exactly four orders of error per two decades of τ, against a solution computed in BigInt rationals from the problem's own optimality conditions. What stops the limit is the solver and not the problem. The normal equations on the weighted problem form entries of size τ², so the constraint block is lost once τ² passes 1/u: the last weight at which they are within an order of the right answer is 10⁴ here, against 1/√u = 6.71·10⁷, and by τ = 10¹⁴ they are wrong by 14.3. Householder QR has no such ceiling and is at 4.8·10⁻¹⁵ at the same weight. Classical Gram–Schmidt is worse than either, at 4.9·10¹⁰.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷10¹⁰log₁₀ τ — the weight on the constraintrelative error against the exact answerτ = 1/√uGram–Schmidtnormal equationsHouseholder QRthe ceiling is the method'sHouseholder, τ = 10¹⁴4.8·10⁻¹⁵normal equations14Gram–Schmidt4.9·10¹⁰1/√u6.7·10⁷a constraint is a weight at infinityand the solver decides how far infinity is
Fig. 13 The three ceilings again, which is what the rule above is a summary of.

At other settings

Minimise ‖Ax − b‖ subject to Bx = d, solved as a weighted least-squares problem, three waysStack the constraint on top of the objective with a weight τ and solve the ordinary least-squares problem that results. In exact arithmetic the answer approaches the constrained one like 1/τ² — measured here as exactly four orders of error per two decades of τ, against a solution computed in BigInt rationals from the problem's own optimality conditions. What stops the limit is the solver and not the problem. The normal equations on the weighted problem form entries of size τ², so the constraint block is lost once τ² passes 1/u: the last weight at which they are within an order of the right answer is 10⁴ here, against 1/√u = 6.71·10⁷, and by τ = 10¹⁴ they are wrong by 1.45. Householder QR has no such ceiling and is at 3.29·10⁻¹⁶ at the same weight. Classical Gram–Schmidt is worse than either, at 5.09·10¹¹.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷10¹⁰log₁₀ τ — the weight on the constraintrelative error against the exact answerτ = 1/√uGram–Schmidtnormal equationsHouseholder QRthe ceiling is the method'sHouseholder, τ = 10¹⁴3.3·10⁻¹⁶normal equations1.4Gram–Schmidt5.1·10¹¹1/√u6.7·10⁷a constraint is a weight at infinityand the solver decides how far infinity is
Fig. 14 Three constraints, where the Householder line settles higher because more of the answer is fixed by the constraint.
Minimise ‖Ax − b‖ subject to Bx = d, solved as a weighted least-squares problem, three waysStack the constraint on top of the objective with a weight τ and solve the ordinary least-squares problem that results. In exact arithmetic the answer approaches the constrained one like 1/τ² — measured here as exactly four orders of error per two decades of τ, against a solution computed in BigInt rationals from the problem's own optimality conditions. What stops the limit is the solver and not the problem. The normal equations on the weighted problem form entries of size τ², so the constraint block is lost once τ² passes 1/u: the last weight at which they are within an order of the right answer is 10⁴ here, against 1/√u = 6.71·10⁷, and by τ = 10¹⁴ they are wrong by 0.575. Householder QR has no such ceiling and is at 6.04·10⁻¹⁶ at the same weight. Classical Gram–Schmidt is worse than either, at 2.69·10¹¹.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷10¹⁰log₁₀ τ — the weight on the constraintrelative error against the exact answerτ = 1/√uGram–Schmidtnormal equationsHouseholder QRthe ceiling is the method'sHouseholder, τ = 10¹⁴6·10⁻¹⁶normal equations0.57Gram–Schmidt2.7·10¹¹1/√u6.7·10⁷a constraint is a weight at infinityand the solver decides how far infinity is
Fig. 15 And four.
The same Householder QR on the same rows, in two ordersOne weighted least-squares problem, one factorisation, and the only difference is whether the 2 heavily weighted constraint rows are placed above the 8 data rows or below them. The algebra is identical, the operation count is identical, the reflectors are the same reflectors. Up to τ = 10⁴ the two answers agree. At τ = 10¹⁴ they are 1.13·10⁻¹⁵ and 0.0125 — a factor of 1.11·10¹³, from the sequence of the rows and nothing else. The reason is that the first reflector is built from the first column, and a column whose large entries are at the bottom has them annihilated last, after the small ones have already been added to them and lost.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ τrelative error against the exact answerconstraint rows lastconstraint rows firstthe same arithmetic, resequencedconstraint first, τ = 10¹⁴1.1·10⁻¹⁵constraint last, same τ0.013the ratio1.1·10¹³multiplications, either order0same reflectors, same countand ten digits between them
Fig. 16 Eight data rows, where the two orderings are closest.
The same Householder QR on the same rows, in two ordersOne weighted least-squares problem, one factorisation, and the only difference is whether the 2 heavily weighted constraint rows are placed above the 20 data rows or below them. The algebra is identical, the operation count is identical, the reflectors are the same reflectors. Up to τ = 10⁴ the two answers agree. At τ = 10¹⁴ they are 6.55·10⁻¹⁶ and 0.00784 — a factor of 1.197·10¹³, from the sequence of the rows and nothing else. The reason is that the first reflector is built from the first column, and a column whose large entries are at the bottom has them annihilated last, after the small ones have already been added to them and lost.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ τrelative error against the exact answerconstraint rows lastconstraint rows firstthe same arithmetic, resequencedconstraint first, τ = 10¹⁴6.6·10⁻¹⁶constraint last, same τ0.0078the ratio1.2·10¹³multiplications, either order0same reflectors, same countand ten digits between them
Fig. 17 And twenty, midway.
How nearly the constraint is satisfied, when a sketch is allowed to see it and when it is notSketching reduces a least-squares problem to 8 rows and keeps its minimiser within (1 + ε) of the original's. That is a statement about a norm, and an equality constraint is a statement that a quantity is zero, which no multiplicative distortion preserves. Kept out of the sketch and imposed exactly, the constraint holds to 1.68·10⁻¹⁶ at every weight. Written as a weight and sketched along with everything else, it holds at 1.747·10⁻⁸ at τ = 10⁸ — and the slope is the finding: the violation falls as 1/τ, one power, where the unsketched weighted problem's falls as 1/τ². The sketch costs exactly half the decades the weight was buying, and buying them back means squaring τ, which the previous figure says the arithmetic stops allowing.2468101210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²log₁₀ τ‖Bx − d‖ / ‖d‖sketched with the objectiveweighted, not sketchedkept out of the sketchone power instead of twokept out — feasibility1.7·10⁻¹⁶sketched at τ = 10⁸1.7·10⁻⁸unsketched at τ = 10⁸1.7·10⁻¹⁶objective ÷ optimum1.2a sketch preserves a normand a constraint is not one
Fig. 18 What happens to the weighted problem when a sketch is allowed to see it, from the randomised field.
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 7 and 8: 1columns 6 and 8: 0.13columns 6 and 7: 0.13modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 8: 4.4·10⁻⁷columns 2 and 8: 2.7·10⁻⁷columns 3 and 8: 2.4·10⁻⁸The two R factors agree to 1.2·10⁻⁶ relative. The two Q factors do not.the 8×8 Hilbert matrixone word, eight orders
Fig. 19 The other sequencing result on this site, in the field that owns it.
The least-squares solution as a projection, with the right angle measuredThe column space drawn edge-on as a plane, the data vector above it, and the perpendicular dropped to the plane, with the residual marked at a right angle to it.everything Ax can reachb = (1.1, 0.4, 1.5)Ax, the closest reachable pointr = b − Ax‖Aᵀr‖ / (‖A‖‖r‖)1.7·10⁻¹⁶‖b‖² − ‖Ax‖² − ‖r‖²1.3·10⁻¹⁵‖r‖1.3200 random nearby points of the plane were tried; none is closer.a 3×2 system, Householder QRperpendicularity is checked
Fig. 20 The unconstrained problem the weighted one is.
How far the coefficients can move without changing the fit, degree 9Relative increase in the residual against relative change in the coefficients, along the least determined direction. The residual does not move measurably until the coefficients have changed by more than a factor of one.10⁻⁴10⁻³10⁻²10⁻¹110¹10²10³10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹relative change in the coefficients, along the worst directionrelative increase in the residualcoefficients doubled39% change, fit unmoved in the sixth digit308×: the third digit movesκ(A) = 3.6·10⁶. Exact arithmetic would pick one point on this floor. It would not raise it.24 points, degree 9, monomial basisthe data leaves them free
Fig. 21 And the case where the objective does not determine the answer at all.
The 12 eigenvalues of a saddle-point matrix with 10 unknowns and 2 constraints, inside their closed-form bracketsK = [[H, Aᵀ], [A, 0]] with κ(H) = 100 and κ(A) = 100. The shaded bands are the Rusten–Winther brackets, computed from four numbers — the extreme eigenvalues of H and the extreme singular values of A — before the matrix was assembled: negative eigenvalues in [-0.995, -10·10⁻⁵] and positive ones in [0.01, 1.62]. The marks are the computed spectrum. There are exactly 10 above zero and 2 below, which is Sylvester's law of inertia and not a property of this matrix: K is congruent to blkdiag(H, −AH⁻¹Aᵀ), both blocks are definite, and congruence preserves signs. The two groups are separated by a gap containing zero, at a ratio of 4.48 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue2 negative10 positivecounted before it was formedpositive10negative2at zero0innermost ratio4.5the zero block is a theoremand so is the count either side of it
Fig. 22 The optimality conditions, as a matrix.
The pivots Cholesky computes on a saddle-point matrix, and the one it stops at — row 11 of 14Every pivot Cholesky forms on K = [[H, Aᵀ], [A, 0]], in order, with the 10 taken inside H drawn as positive and the 11-th, which is not, drawn as the failure. It is not a near miss and it is not a rounding accident: the leading n × n block is positive definite, so the first 10 pivots are positive at any conditioning, and the next one is 0 − a₁ᵀH⁻¹a₁ where a₁ᵀ is the first row of A. That equals −(AH⁻¹Aᵀ)₁₁ = -3.7312, and the pivot the routine actually reports is -3.7312 — the same number to 15 digits. The failure has an address and a value, and both were available before the factorisation started.02468101210⁻³10⁻²10⁻¹110¹10²pivot, in order|pivot|negative, at row 11|S₁₁| from AH⁻¹Aᵀa failure with an addressrows before it10the pivot-3.7−(AH⁻¹Aᵀ)₁₁-3.7disagreement3.6·10⁻¹⁶it does not fail somewhereit fails at the constraint
Fig. 23 Which nothing in this field’s usual toolkit can factorise.
Leverage and the deleted residual for 40 observations of a 6-column fitThe upper panel is the diagonal of the hat matrix, one bar per observation, with the average p/m = 0.150 drawn through it; the leverages sum to 6.000000000, which is exactly the number of columns. The lower panel is the leave-one-out residual, computed in closed form as eᵢ/(1 − hᵢ) and, separately, by refitting the model 40 times without each observation; the two agree to 1.3·10⁻¹². The first observation carries a leverage of 0.5000 by construction, and 1 − h is the number a hyperbolic downdate takes the square root of.the diagonal of the hat matrix, hᵢ = aᵢᵀ(AᵀA)⁻¹aᵢaverage p/m = 0.15010the leave-one-out residual: eᵢ/(1 − hᵢ), and forty refitsbars: closed form · dots: refitted without that pointone number, two fieldsΣ hᵢ, exactly p6largest leverage0.5closed form against refits1.3·10⁻¹²1 − h of the first row0.5y appears in the residualand nowhere in the leverage
Fig. 24 A quantity from the least-squares field that a constraint changes and a weight only nearly changes.
QᵀQ from classical Gram–Schmidt and from Householder on the 8×8 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 = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle
Fig. 25 And the number that decides whether a QR of the stacked matrix was a QR at all.
Minimise ‖Ax − b‖ subject to Bx = d, solved as a weighted least-squares problem, three waysStack the constraint on top of the objective with a weight τ and solve the ordinary least-squares problem that results. In exact arithmetic the answer approaches the constrained one like 1/τ² — measured here as exactly four orders of error per two decades of τ, against a solution computed in BigInt rationals from the problem's own optimality conditions. What stops the limit is the solver and not the problem. The normal equations on the weighted problem form entries of size τ², so the constraint block is lost once τ² passes 1/u: the last weight at which they are within an order of the right answer is 10⁴ here, against 1/√u = 6.71·10⁷, and by τ = 10¹⁴ they are wrong by 9.57. Householder QR has no such ceiling and is at 8.97·10⁻¹⁶ at the same weight. Classical Gram–Schmidt is worse than either, at 1.69·10¹¹.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹10⁴10⁷10¹⁰log₁₀ τ — the weight on the constraintrelative error against the exact answerτ = 1/√uGram–Schmidtnormal equationsHouseholder QRthe ceiling is the method'sHouseholder, τ = 10¹⁴9·10⁻¹⁶normal equations9.6Gram–Schmidt1.7·10¹¹1/√u6.7·10⁷a constraint is a weight at infinityand the solver decides how far infinity is
Fig. 26 Thirty data rows, where the constant in front of 1/τ² moves and the ceilings do not.
The same Householder QR on the same rows, in two ordersOne weighted least-squares problem, one factorisation, and the only difference is whether the 3 heavily weighted constraint rows are placed above the 12 data rows or below them. The algebra is identical, the operation count is identical, the reflectors are the same reflectors. Up to τ = 10⁴ the two answers agree. At τ = 10¹⁴ they are 3.29·10⁻¹⁶ and 0.00146 — a factor of 4.431·10¹², from the sequence of the rows and nothing else. The reason is that the first reflector is built from the first column, and a column whose large entries are at the bottom has them annihilated last, after the small ones have already been added to them and lost.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ τrelative error against the exact answerconstraint rows lastconstraint rows firstthe same arithmetic, resequencedconstraint first, τ = 10¹⁴3.3·10⁻¹⁶constraint last, same τ0.0015the ratio4.4·10¹²multiplications, either order0same reflectors, same countand ten digits between them
Fig. 27 Three constraint rows to misplace rather than two.
How nearly the constraint is satisfied, when a sketch is allowed to see it and when it is notSketching reduces a least-squares problem to 12 rows and keeps its minimiser within (1 + ε) of the original's. That is a statement about a norm, and an equality constraint is a statement that a quantity is zero, which no multiplicative distortion preserves. Kept out of the sketch and imposed exactly, the constraint holds to 1.19·10⁻¹⁶ at every weight. Written as a weight and sketched along with everything else, it holds at 1.587·10⁻⁸ at τ = 10⁸ — and the slope is the finding: the violation falls as 1/τ, one power, where the unsketched weighted problem's falls as 1/τ². The sketch costs exactly half the decades the weight was buying, and buying them back means squaring τ, which the previous figure says the arithmetic stops allowing.2468101210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²log₁₀ τ‖Bx − d‖ / ‖d‖sketched with the objectiveweighted, not sketchedkept out of the sketchone power instead of twokept out — feasibility1.2·10⁻¹⁶sketched at τ = 10⁸1.6·10⁻⁸unsketched at τ = 10⁸1.7·10⁻¹⁶objective ÷ optimum1.2a sketch preserves a normand a constraint is not one
Fig. 28 What a sketch does to the same construction, from the randomised field.
The 14 eigenvalues of a saddle-point matrix with 10 unknowns and 4 constraints, inside their closed-form bracketsK = [[H, Aᵀ], [A, 0]] with κ(H) = 100 and κ(A) = 100. The shaded bands are the Rusten–Winther brackets, computed from four numbers — the extreme eigenvalues of H and the extreme singular values of A — before the matrix was assembled: negative eigenvalues in [-0.995, -10·10⁻⁵] and positive ones in [0.01, 1.62]. The marks are the computed spectrum. There are exactly 10 above zero and 4 below, which is Sylvester's law of inertia and not a property of this matrix: K is congruent to blkdiag(H, −AH⁻¹Aᵀ), both blocks are definite, and congruence preserves signs. The two groups are separated by a gap containing zero, at a ratio of 39.3 between the innermost positive eigenvalue and the innermost negative one.-1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue4 negative10 positivecounted before it was formedpositive10negative4at zero0innermost ratio39the zero block is a theoremand so is the count either side of it
Fig. 29 And the exact system the weight is approximating.

Where this goes next

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. -1-0.582271-0.1645420.2531860.6709151.088641.506370eigenvalue4 negative10 positivecounted before it was formedpositive10negative4at zero0innermost ratio39the zero block is a theoremand so is the count either side of it The zero that is not a missing entry A constrained minimisation produces a matrix with a zero block, and the zero is a theorem rather than a sparsity pattern. No pivot order makes it positive definite, no precision changes that, and Cholesky does not fail somewhere on it — it fails at the first constraint row, on a number the problem already contained. 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 7 and 8: 1columns 6 and 8: 0.13columns 6 and 7: 0.13modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 8: 4.4·10⁻⁷columns 2 and 8: 2.7·10⁻⁷columns 3 and 8: 2.4·10⁻⁸The two R factors agree to 1.2·10⁻⁶ relative. The two Q factors do not.the 8×8 Hilbert matrixone word, eight orders Two Gram–Schmidts One argument changes. Classical Gram–Schmidt projects the original column onto each previous direction; modified projects what is left of it. In exact arithmetic the coefficients are identical. In floating point they differ by eight orders of magnitude in the thing that matters.

What links here

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

Reads more easily once this is understood

Essays that name this one as worth reading first.

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.

Condition numberEquality constrained least squaresExact ground truthHouseholderLeast squaresMethod of weightingNormal equationsQR factorisationRow orderingSaddle-point systems