Orthogonality, Projections, and Least Squares

Many machine-learning problems ask for the best approximation inside a restricted family. Linear regression is one. The prediction vector has to lie in the column space of the design matrix. Orthogonality supplies the certificate. At the least-squares solution the residual is perpendicular to every direction the model can produce, so the model cannot cut squared error by moving inside its own subspace. That single fact explains the normal equations, the role of rank, and why good software never forms an inverse explicitly.

The geometry is old. The arithmetic is not settled. Legendre put the normal equations in print in 1805, in an appendix to his book on the orbits of comets. Two centuries later the National Institute of Standards and Technology still publishes least-squares datasets with answers certified to 16 significant digits, because widely used software gets them wrong.

One mathematical objective, four methods, four different numerical behaviours. The difference is measured, not argued. The National Institute of Standards and Technology publishes Statistical Reference Datasets whose least-squares answers are certified to 16 significant digits. Filip is one of them: 82 observations fitted by a tenth-degree polynomial, 11 parameters, classified “Higher Level of Difficulty.” Its certified values include β0 = −1467.48961422980 (SD 298.084530995537), a residual standard deviation of 0.334801051324544E-02 and R² = 0.996727416185620.

Someone checked. A 2009 audit ran four freely used packages against those certified digits. On Filip, R's lm() returned no usable answer at all, and neither did Ox. Octave's regress scored a log-relative-error of 7 digits on its least accurate coefficient and 1 digit on the residual standard deviation. NIST states plainly why the datasets exist: “Even with the availability of reliable code for linear least squares fitting, problems persist. Failure to use the best algorithms and to implement them most effectively is often the cause.” The card that calls normal equations reasonable only in well-behaved small problems is describing that table.

The iterative column is not generic either. It has a name and a date: LSQR, published by Paige and Saunders in 1982. It is built on the Golub–Kahan bidiagonalization, the same reduction that sits under the dense SVD driver in the third column. It touches A only through the products Ax and Aᵀx, which is what lets it run on matrices too large to factor. Its authors' claim for it is precise: “It is analytically equivalent to the standard method of conjugate gradients, but possesses more favorable numerical properties.” They report it as the most reliable of the conjugate-gradient-type methods they tested when A is ill-conditioned. scipy.sparse.linalg.lsqr implements it, and accepts A as a linear operator producing Ax and Aᵀx.

Minimizing ‖Xβ − y‖² gives the gradient 2Xᵀ(Xβ − y). Setting it to zero yields XᵀXβ = Xᵀy. If X has full column rank, XᵀX is invertible and the solution is unique. When columns are dependent, several parameter vectors produce the same projection.

Production software does not leave that ambiguity open. It settles it by a documented rule. LAPACK describes its driver in one line: “DGELSD computes the minimum-norm solution to a real linear least squares problem: minimize 2-norm(| b - A*x |) using the singular value decomposition (SVD) of A. A is an M-by-N matrix which may be rank-deficient.” Among the many minimisers it returns the one of smallest 2-norm. It also reports an integer effective rank, because “The effective rank of A is determined by treating as zero those singular values which are less than RCOND times the largest singular value.” That is why numpy.linalg.lstsq documents that “If there are multiple minimizing solutions, the one with the smallest 2-norm ||x|| is returned”. It is also why scipy.linalg.lstsq defaults to lapack_driver='gelsd' and hands back the effective rank of a alongside the coefficients. The alternatives are 'gelsy' and 'gelss'.

So the fitted values can be unique when the coefficients are not. And the coefficient vector you receive is a choice, made by a cutoff you are allowed to set. Telling those two apart matters whenever parameters are interpreted rather than used only for prediction.

A small calculation makes the orthogonality condition visible. The device that makes it work — centring the feature — is not a classroom convenience. It is the remedy a standards body prescribes when a real fit falls apart.

NIST's other “Higher Level of Difficulty” real-world least-squares case is Longley: 16 observations, 7 parameters, certified at β0 = −3482258.63459582 with a residual standard deviation of 304.854073561965. The data comes from a 1967 paper by J. W. Longley appraising least squares programs for the electronic computer from the point of view of the user. R's datasets manual describes it as “A macroeconomic data set which provides a well-known example for a highly collinear regression” — 7 variables, 1947–1962, n=16. Measured against the certified digits, Octave's regress scored a log-relative-error of 0 on that residual standard deviation. No correct digits at all. R's lm scored 14.

Use these checks before interpreting a fitted linear model. Several of them have documented answers rather than judgement calls.

Step 2 has a prescribed first move. Near-dependence among columns is exactly what makes Longley hard, and NIST's advice is to centre the predictors — subtract each predictor's mean — and rerun before concluding anything else. Step 5 has a return value: scipy.linalg.lstsq reports the effective rank that GELSD derived from its RCOND cutoff on the singular values, so a rank below the column count is a fact the library will state if you read it. Steps 4 and 6 have Filip and Longley standing behind them. On those two certified datasets Octave's regress returned a residual standard deviation with 1 correct digit and 0 correct digits. R's lm() returned nothing usable on Filip at all. Every one of those runs still produced output shaped like a fitted model, with coefficients and a residual line and no warning that the digits were gone.

A vector splits into a component inside a subspace and a residual orthogonal to it. The target y sits outside col(X). The projection ŷ is the closest reachable point under Euclidean distance. The residual r = y − ŷ carries what the subspace cannot reach. And Xᵀr = 0 is the certificate that no move inside the subspace improves the fit.

Rank is where the picture meets the software. When col(X) is spanned by fewer directions than X has columns, many β land on the same ŷ. DGELSD picks among them by an explicit rule: the minimum 2-norm solution, with singular values below RCOND times the largest treated as zero.