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

The projection and the right angle

The least-squares solution is the one whose residual is perpendicular to everything the columns can reach. That is not a mnemonic — it is an equation, Aᵀr = 0, and the computed answer satisfies it to 10⁻¹⁶.
15 min read 6 figures Measured, not assertedTwo routes to a numberWhere the algebra stops being the arithmetic

There are more equations than unknowns and no exact solution. A has m rows and n columns with m > n, and b is not in the column space, so Ax = b has no answer. The least-squares problem asks for the x that comes closest: the one minimising ‖b − Ax‖.

The answer has a geometric characterisation that is worth more than the formula. Ax ranges over the column space of A — a plane through the origin in m-dimensional space — and the closest point of a plane to a point off it is the foot of the perpendicular. So the residual r = b − Ax is perpendicular to the column space, which means perpendicular to every column of A, which is the equation Aᵀr = 0.

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. 1 The column space drawn edge-on, the data vector above it, and the perpendicular dropped to it. The badge reports the three quantities that make this a measurement rather than a diagram: ‖Aᵀr‖ scaled by the sizes of A and r, the Pythagorean discrepancy, and the residual norm itself.

The perpendicularity is measured

Aᵀr = 0 is the defining property, so it is the thing to check, and the check is not “the picture shows a right angle”.

For the three-by-two problem in the figure, ‖Aᵀr‖ divided by ‖A‖‖r‖ — a scale-free version of the condition — comes out below 10⁻¹⁶. The computed answer satisfies the defining equation of the problem to the last bit the arithmetic has.

That number is worth having on the figure because it is the thing that could be wrong. A least-squares routine with a sign error, a transposition, or a mishandled permutation still returns a vector, and the vector still looks plausible; what it does not do is make the residual perpendicular. One matrix–vector product distinguishes them.

Pythagoras as the second route

Perpendicularity implies a decomposition of the length of b:

‖b‖² = ‖Ax‖² + ‖r‖²

which is the theorem of Pythagoras applied to the triangle in the figure. It follows from Aᵀr = 0 by expanding, so it is not independent as mathematics — but as arithmetic it is a completely different computation, involving three norms rather than a matrix–vector product, and it must agree to 10⁻¹⁴.

The reason to compute both is the reason two routes to a number runs through this whole site. A single check confirms that the code is self-consistent. Two checks sharing no arithmetic confirm something about the world. If a scaling error had crept into the residual, the perpendicularity check would still pass — the error would scale out — and the Pythagorean one would not.

And nothing closer exists

The third check is the one that turns the argument from a derivation into a measurement.

The claim is that x minimises ‖b − Ax‖. So: perturb x in two hundred random directions, of a size comparable to x itself, and measure the residual at each. Every one must be at least as large as the one at the computed answer.

All two hundred are. That is not a proof — a random search never is — but it is the kind of evidence that catches the errors a derivation cannot, and it costs a few hundred matrix–vector products at build time. A sign flip that produced a maximum rather than a minimum would pass the perpendicularity check, pass Pythagoras, and fail this immediately.

Why the projection is the whole story

Three consequences follow from the geometry, and each answers a question people bring to least squares.

The fitted values are a projection, so they are what the model can see. Ax̂ is the part of b that the columns can express; r is the part they cannot. That is the honest statement of what a fit is, and it makes “the model explains the data” a statement about the angle between b and a subspace.

The residual is orthogonal to every predictor. This is why a residual plot against any column of A is expected to show no trend: there is no trend, exactly, by construction. A trend visible in such a plot means the fit is not the least-squares fit, or that the plot is against something not in the column space — which is how omitted-variable problems are detected.

The projection is unique even when x is not. If the columns are dependent, many x give the same Ax̂, and the projection is still a single well-defined point. That is the distinction between the fit being determined and the coefficients being determined, and it is the entire subject of the valley with no bottom — where the columns are not dependent, only nearly so, and the coefficients are correspondingly nearly undetermined.

The two roads from here

Aᵀr = 0 expands to Aᵀ(b − Ax) = 0, that is, AᵀAx = Aᵀb — the normal equations. They are correct, they are the shortest derivation, and forming AᵀA is the single worst thing that can be done to this problem.

The reason: the singular values of AᵀA are the squares of those of A, so κ(AᵀA) = κ(A)². A fit with κ(A) = 10⁶ — which is an ordinary polynomial fit — becomes a linear system with κ = 10¹², and half the digits available are spent before the solve begins.

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. 2 The two roads, measured. Above the axis, the QR route holds accuracy across the whole range. Below it, the normal-equations route degrades and then breaks outright at the ε where 1 + ε² rounds to 1 and the cross-product matrix becomes exactly singular.

The other road keeps the geometry. Factor A = QR with Q orthogonal, and the projection onto the column space is QQᵀ, so the residual condition becomes Rx = Qᵀb — a triangular system with κ® = κ(A), not κ(A)². The condition number is untouched, because every step was orthogonal and orthogonal steps have κ = 1.

The road that squares the problem works that through with numbers, including the exact ε at which the normal equations stop being a method at all.

What orthogonality buys here, concretely

This is where orthogonal is a number becomes practical rather than a matter of principle.

The QR route computes x from Rx = Qᵀb, and that derivation is valid only if QᵀQ = I. Feed it a Q from classical Gram–Schmidt on an ill-conditioned matrix — where ‖QᵀQ − I‖ can be of order one — and the method silently reverts to something no better than the normal equations it was chosen to avoid.

The chain is: orthogonal Q → the projection QQᵀ is a projection → Rx = Qᵀb solves the right problem → the conditioning is κ(A) rather than κ(A)². Break the first link and everything after it is a different algorithm wearing the same name. That is why the loss-of-orthogonality measurement matters downstream and not only as a curiosity about factorisations.

Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 3 The first link, measured. Householder’s Q is orthogonal to 10⁻¹⁵ whatever the matrix; classical Gram–Schmidt’s is not, at any condition number above about 10⁵. A least-squares routine built on the second one inherits a κ² behaviour it was designed to avoid.

What is asserted here

Four measurements, all on the figure and all checked when it is generated.

Aᵀr = 0, scaled, below 10⁻¹⁵. The defining property.

‖b‖² − ‖Ax‖² − ‖r‖² = 0 to 10⁻¹⁴. The second route, sharing no arithmetic with the first.

No nearby point is closer. Two hundred random perturbations, all giving a residual at least as large.

The factorisation carries its residual, because this figure factorises a matrix and the site’s rule applies without exception: no decomposition is drawn without the number that says the factors factorise what they were given.

The four together are the difference between a diagram of the projection theorem and a measurement of it. The diagram would be identical if the code were wrong.

Where this goes

The geometry above is exact and it is the easy part. Two things complicate it, and they are the remaining essays in this field.

The conditioning of the problem — how well the coefficients are determined — is a separate question from whether the projection is computed correctly, and it is what the road that squares the problem and the valley with no bottom are about. The projection can be computed to sixteen digits while the coefficients producing it are determined to two, and both statements are about the same fit.

And the basis — which columns were chosen to span the space — is a modelling decision that changes κ(A) by eight orders of magnitude without changing the fitted curve at all. The column space is what determines the projection; the particular columns used to describe it determine whether the coefficients mean anything.

The pseudoinverse, and what it does when the answer is not unique

Everything above assumed A has full column rank, so the least-squares solution is unique. When it does not — or nearly does not, which is the case that actually arises — the projection is still a single well-defined point and the coefficients producing it are not.

The standard resolution is the pseudoinverse, A⁺ = VΣ⁺Uᵀ, where Σ⁺ inverts the nonzero singular values and leaves the zeros alone. Among all the x that achieve the minimum residual, A⁺b is the one of smallest norm — a second criterion imposed to break the tie, and a reasonable one, but a choice rather than a consequence.

What makes the pseudoinverse genuinely useful is that “zero” becomes a decision, and the decision is the same one rank is a decision is about. Setting 1/σ to zero for σ below a tolerance is exactly truncating the directions the data does not determine, and it turns an unstable computation into a stable answer to a slightly different question.

The alternative — inverting a singular value of 10⁻¹⁶ and getting 10¹⁶ — is what a solver without a tolerance does, and the result is a coefficient vector of enormous norm producing the same projection. The fit is identical; the coefficients are noise multiplied by 10¹⁶.

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 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 4 Where the tolerance goes. Ten singular values, four of them signal and six of them noise, with the cutoff marked. Everything to the right of it is a direction the pseudoinverse will discard, and the gap says how confident that decision is.

Weighted least squares, in one paragraph

The geometry generalises without changing, which is worth knowing because it covers most of applied statistics.

If the observations have different reliabilities, the right objective is not ‖b − Ax‖ but a weighted norm, and the perpendicularity condition becomes AᵀW(b − Ax) = 0 for the weight matrix W. That is the same projection theorem in a different inner product: the residual is orthogonal to the column space with respect to W, and everything above holds with the word “perpendicular” reinterpreted.

Numerically the sensible route is to scale the rows by √W and run ordinary least squares on the scaled problem, which keeps the QR machinery and its conditioning. Forming AᵀWA is the same mistake as forming AᵀA, one weight matrix later.

What the residual is for

A closing distinction, because the word does three jobs on this site and this is the one place they meet.

The least-squares residual ‖b − Ax‖ at the optimum is not an error — it is the answer. It says how far the data is from the model’s reach, and a large one means the model is wrong rather than the computation. This is the quantity everything in this field is about.

The factorisation residual ‖A − QR‖/‖A‖ says whether the factors factorise, and is the site’s rule.

The solve residual ‖b − Ax̂‖ for a square system says whether the algorithm did its job and nothing about the answer, which is a small residual is not a small error.

Three quantities, one word. The least-squares one is the only one of the three that is supposed to be large.

Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 5 And the reason to keep them apart. The lower line is a factorisation and solve residual — small, flat, and about the algorithm. The upper line is the error in the answer. A least-squares residual would belong to neither: it is a property of the data, and it does not move when the arithmetic does.
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. 6 The link on which all of this depends, printed. QᵀQ from two factorisations of the same matrix: the projection QQᵀ is a projection only for the lower one. Every claim above about the geometry being computed to 10⁻¹⁶ is a claim about which Q was used.

Two further threads run out of this field. Orthogonal is a number is the measurement that guards the link above, and an answer that is known is what makes it possible to say that a fitted coefficient is wrong rather than merely different.

Why the geometry is worth insisting on

There is a version of this material that never mentions a plane: define the objective, differentiate, set the gradient to zero, obtain AᵀAx = Aᵀb, solve. It is shorter, it is correct, and it loses two things that turn out to matter.

It loses the reason the answer exists and is unique. The projection of a point onto a subspace exists because a closed subspace is closed, and it is unique because the subspace is convex — facts that survive when A is rank deficient and the algebraic derivation quietly does not, since AᵀA is then singular and the manipulation divides by nothing.

And it loses the diagnostic. Aᵀr = 0 is a property of the computed answer that can be checked in one matrix–vector product, and the calculus route offers nothing comparable: a gradient that was set to zero symbolically leaves no residue to inspect. Every figure in this field prints a number that comes from the geometry, and there is no corresponding number available from the derivation.

The general form of that observation is worth carrying past least squares. A characterisation — the answer is the thing with this property — gives a test. A construction — the answer is what this procedure produces — gives only a procedure, and a procedure cannot be verified except by comparing it against another procedure. Every check on this site is a characterisation used as a test: the factorisation must reconstruct, the eigenvector must satisfy its own equation, the residual must be perpendicular. None of those is available from the derivation that produces the object.

What the two hundred random points are for

The third check in the badge deserves a defence, because a random search is the least rigorous thing on this site and it is here deliberately.

The perpendicularity condition and Pythagoras both follow from the same derivation, so as evidence they are correlated: an error in the interpretation of the derivation — a sign, an orientation, a minimisation that has become a maximisation — would satisfy both. What distinguishes a minimum from a maximum is not an equation; it is the behaviour of the objective nearby, and the only way to observe behaviour nearby is to look nearby.

Two hundred points is not a proof and it is not offered as one. It is a test that can fail, that no derivation-derived check can fail in the same way, and that costs a few hundred matrix–vector products at build time. The site’s habit is to prefer a test that can fail over a test that is elegant, and where both are available to run both.

That is also why the count is stated on the figure rather than left in the code. A reader who thinks two hundred is too few knows how many were tried; a reader who thinks the check is worthless knows exactly what it does. Neither is possible if the caption says only that the answer is optimal.

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. 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 The valley with no bottom A degree-nine fit's coefficients can be moved by a third of their own size before the residual changes in the sixth significant figure. The arithmetic did not lose those digits. The data never contained them. 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 Orthogonal is a number "Q is orthogonal" is a claim about a measurable quantity, ‖QᵀQ − I‖, and on the eight-by-eight Hilbert matrix two standard algorithms return 10⁻¹⁵ and 1 for it. The one that returns 1 still reconstructs the matrix perfectly, which is why nothing warns you.