Least squares

The nearest-neighbours note ended with a model that answers every query from a handful of nearby points and learns nothing about the relationship as a whole; in many dimensions the nearby points disappear. This note takes the opposite bet. It assumes that the target is, up to noise, a weighted sum of the features, one weight per feature, and fits the weights to all the data at once by minimizing the squared errors. The fit has a closed form, a geometric meaning and a probabilistic meaning, and each of the three explains a different way it can fail: on collinear features, on badly scaled numbers and on outliers. The failures set up the rest of the section.

Notation. The training set is $(x_1, y_1), \dots, (x_n, y_n)$ with $x_i \in \mathbb R^d$ and $y_i \in \mathbb R$. The design matrix $X \in \mathbb R^{n \times d}$ has the $x_i^\top$ as rows and $y \in \mathbb R^n$ collects the targets. When the model has an intercept, the first feature is the constant 1, so $d$ counts the intercept. A weight vector $w \in \mathbb R^d$ predicts $\hat y_i = \langle w, x_i\rangle$, all at once $Xw$. The residuals are $r = y - Xw$, the residual sum of squares is $\mathrm{RSS}(w) = \lVert y - Xw\rVert^2$, and the mean squared error is $\mathrm{MSE} = \mathrm{RSS}/n$. $\hat w$ is the minimizer and $\hat y = X\hat w$.

A linear model#

The model is

$$a(x) = \langle w, x\rangle = w_1 x^{(1)} + \dots + w_d x^{(d)} .$$

Each weight says how much the prediction changes when one feature grows by one unit and the others stay fixed. The model is linear in $w$, not necessarily in the original inputs: the features can be powers of a number, logarithms, products of two inputs or the one-hot columns of a category from the features note. Everything below holds for any such fixed set of features, because it only uses the matrix $X$.

The loss is the squared error, averaged over the training set:

$$Q(w) = \frac1n \sum_{i=1}^n \big(\langle w, x_i\rangle - y_i\big)^2 = \frac1n \lVert Xw - y\rVert^2 .$$

Why squares rather than absolute values or another penalty? Three answers follow in this note: the minimizer has a formula, it is an orthogonal projection, and it is the maximum-likelihood estimate when the noise is Gaussian. The last section shows the price.

The normal equations#

Expand the residual sum of squares: $\lVert y - Xw\rVert^2 = y^\top y - 2 w^\top X^\top y + w^\top X^\top X w$. The gradient of the linear term is $-2X^\top y$ and that of the quadratic form, with $X^\top X$ symmetric, is $2X^\top X w$, so

$$\nabla_w \lVert y - Xw\rVert^2 = 2X^\top (Xw - y) .$$

The Hessian is $2X^\top X$, and $v^\top X^\top X v = \lVert Xv\rVert^2 \ge 0$ for every $v$: the function is convex, and every point where the gradient vanishes is a global minimum. Setting the gradient to zero gives the normal equations

$$X^\top X\, \hat w = X^\top y .$$

If the columns of $X$ are linearly independent (which needs $n \ge d$), then $\lVert Xv\rVert^2 > 0$ for every $v \ne 0$, the matrix $X^\top X$ is positive definite and invertible, and the minimizer is unique:

$$\hat w = (X^\top X)^{-1} X^\top y .$$

If the columns are dependent, for example when one feature is the sum of two others or when $d > n$, the normal equations have infinitely many solutions, all with the same fitted values; the section on stability returns to that case.

Worked example. Four points: $x = (0, 1, 2, 3)$ and $y = (1, 3, 2, 5)$, model $\hat y = w_1 + w_2 x$. The design matrix has columns $(1, 1, 1, 1)$ and $(0, 1, 2, 3)$, so

$$X^\top X = \begin{pmatrix} 4 & 6 \\ 6 & 14 \end{pmatrix}, \qquad X^\top y = \begin{pmatrix} 11 \\ 22 \end{pmatrix} .$$

The determinant is $4 \cdot 14 - 6 \cdot 6 = 20$, and $\hat w = \frac{1}{20}\begin{pmatrix} 14 & -6 \\ -6 & 4\end{pmatrix}\begin{pmatrix} 11 \\ 22\end{pmatrix} = (1.1, 1.1)$. The fitted line is $\hat y = 1.1 + 1.1x$, the residuals are $(-0.1, 0.8, -1.3, 0.6)$ and $\mathrm{RSS} = 0.01 + 0.64 + 1.69 + 0.36 = 2.7$.

Solving the normal equations costs $O(nd^2)$ to form $X^\top X$ and $O(d^3)$ to solve the $d \times d$ system. For ten features that is nothing; for a million sparse text features it is impossible, which is where the next note starts.

The geometry: a projection#

The normal equations can be rewritten as $X^\top (y - X\hat w) = 0$: the residual vector is orthogonal to every column of $X$. Read the $n$ targets as one point $y \in \mathbb R^n$. As $w$ varies, the predictions $Xw$ sweep the column space $\operatorname{col}(X)$, a $d$-dimensional subspace of $\mathbb R^n$. Least squares picks the point of that subspace closest to $y$, and the closest point is the foot of the perpendicular.

Figure

Least squares as a projection

The orthogonality gives a proof of optimality without calculus. For any $w$, split $y - Xw = (y - X\hat w) + X(\hat w - w)$. The first term is orthogonal to the column space and the second lies in it, so by Pythagoras

$$\lVert y - Xw\rVert^2 = \lVert y - X\hat w\rVert^2 + \lVert X(\hat w - w)\rVert^2 \ \ge\ \lVert y - X\hat w\rVert^2 ,$$

with equality exactly when $Xw = X\hat w$. In the three-point example of the figure the extra term is $0.5$ for $w = (1, 1)$.

The map from $y$ to $\hat y$ is linear: $\hat y = Hy$ with the hat matrix $H = X(X^\top X)^{-1}X^\top$. It is the orthogonal projection onto $\operatorname{col}(X)$, so it is symmetric and idempotent ($H^2 = H$: projecting twice changes nothing), and its trace is $\operatorname{tr}\big((X^\top X)^{-1}X^\top X\big) = d$, the dimension of the subspace. The diagonal entry $h_{ii}$, the leverage of point $i$, is how much $\hat y_i$ moves when $y_i$ moves by one. In the four-point example the leverages are $(0.7, 0.3, 0.3, 0.7)$, summing to $d = 2$: the end points pull on the line more than the middle ones.

Two consequences of orthogonality are used constantly. When the model has an intercept, the constant column is one of the columns, so $\sum_i r_i = 0$: the residuals average to zero and the mean of $\hat y$ equals the mean of $y$. In the worked example, $-0.1 + 0.8 - 1.3 + 0.6 = 0$, and also $0 \cdot (-0.1) + 1 \cdot 0.8 + 2 \cdot (-1.3) + 3 \cdot 0.6 = 0$. And $y = \hat y + r$ with $\hat y \perp r$, so the total variation splits into an explained and an unexplained part, which is what $R^2$ below measures.

The probability: Gaussian noise and the conditional mean#

Suppose the data come from $y_i = \langle w, x_i\rangle + \varepsilon_i$ with $\varepsilon_i$ independent and $N(0, \sigma^2)$. The negative log-likelihood is $\mathrm{RSS}(w)/(2\sigma^2)$ plus a constant, as derived in the maximum likelihood note, so minimizing squared error is maximizing the likelihood. Without assuming a model at all, squared error has another justification: its population minimizer over all functions is the conditional mean $E[Y \mid X = x]$ (the best predictor). Least squares restricted to linear functions estimates the linear function closest to $E[Y \mid X]$ in mean square, which is $E[Y \mid X]$ itself when that is linear.

The distribution of the coefficients. Here the design $X$ is treated as fixed and only the noise is random: the statements below are about repeated draws of $\varepsilon$ for the same inputs. Substitute $y = Xw + \varepsilon$ into the formula:

$$\hat w = (X^\top X)^{-1}X^\top (Xw + \varepsilon) = w + (X^\top X)^{-1}X^\top \varepsilon .$$

The estimate is the true weight plus a fixed linear map of the noise. Three facts follow, each under its own assumption:

• Unbiased if $E[\varepsilon] = 0$: $E[\hat w] = w$.
• Covariance $\sigma^2 (X^\top X)^{-1}$ if the noise is uncorrelated with equal variances, $\operatorname{Var}\varepsilon = \sigma^2 I$: with $A = (X^\top X)^{-1}X^\top$, $\operatorname{Var}\hat w = A\,\sigma^2 I\,A^\top = \sigma^2 (X^\top X)^{-1}$.
• Gaussian if the noise is Gaussian, since an affine map of a Gaussian vector is Gaussian: $\hat w \sim N\big(w, \sigma^2 (X^\top X)^{-1}\big)$. For other independent noise with finite variance, the distribution is approximately Gaussian for large $n$ by the central limit theorem, provided no single point dominates the fit: the largest leverage $\max_i h_{ii}$ (defined above) must tend to 0.

The noise variance is unknown and is estimated from the residuals. Since $r = (I - H)y = (I - H)\varepsilon$ (because $(I - H)X = 0$), $\mathrm{RSS} = \varepsilon^\top (I - H)\varepsilon$ and $E[\mathrm{RSS}] = \sigma^2 \operatorname{tr}(I - H) = \sigma^2 (n - d)$. So

$$\hat\sigma^2 = \frac{\mathrm{RSS}}{n - d}$$

is unbiased: $d$ degrees of freedom are spent on the fit, as the $n - 1$ of the sample variance spends one on the mean. The square roots of the diagonal of $\hat\sigma^2 (X^\top X)^{-1}$ are the standard errors of the coefficients.

Worked example. In the four-point example $\hat\sigma^2 = 2.7/(4 - 2) = 1.35$ and $(X^\top X)^{-1} = \frac{1}{20}\begin{pmatrix} 14 & -6 \\ -6 & 4\end{pmatrix}$, so the slope $1.1$ has standard error $\sqrt{1.35 \cdot 4/20} = 0.520$ and the intercept $\sqrt{1.35 \cdot 14/20} = 0.972$. With four points the slope is barely two standard errors from zero. A simulation confirms the formula: 30 equally spaced points on $[0, 1]$, true $w = (1, 2)$, $\sigma = 0.5$ and 20,000 noise draws give coefficient means $(1.001, 2.000)$ and standard deviations $(0.177, 0.303)$, against $(0.178, 0.306)$ from $\sigma\sqrt{\operatorname{diag}(X^\top X)^{-1}}$ (the Code section).

The Gauss–Markov theorem adds that, with uncorrelated equal-variance noise, no other estimator that is linear in $y$ and unbiased has a smaller variance for any linear combination of the coefficients (Hastie, Tibshirani and Friedman, 2009, section 3.2.2). The qualifier matters: estimators that accept a little bias can have a much smaller variance and a smaller mean squared error, the trade-off of the estimators note. The regularization note builds such estimators.

When the solution is unstable#

The covariance $\sigma^2 (X^\top X)^{-1}$ is large in the directions where $X^\top X$ has small eigenvalues, that is, where some combination of features is nearly constant across the data. Two features that almost repeat each other produce exactly that.

Collinearity. Take 100 points with $x^{(2)} = x^{(1)} + 0.01\,\xi$, where $\xi$ is standard normal noise, and $y = x^{(1)} + x^{(2)} + \varepsilon$ with $\sigma = 1$. The correlation of the two features is $0.99994$ and the condition number of $X^\top X$ is $3.5 \cdot 10^4$. Three draws of the noise give the weights $(-5.00, 6.88)$, $(11.88, -9.90)$ and $(-8.83, 10.94)$, with sums $1.88$, $1.98$ and $2.11$. The data pin down the sum $w_1 + w_2$, the direction in which the features vary, and say almost nothing about the difference, the direction in which they do not. The standard error of each weight is $9.74$, against $0.10$ for a model with one of the two features alone. The predictions on data like the training data hardly change from draw to draw; the coefficients, and any interpretation of them, are noise. A weight of $-9.90$ on a feature that helps the target is the typical symptom.

Rank deficiency. When the columns are exactly dependent, or $d > n$, $X^\top X$ is singular and the normal equations have a whole affine set of solutions. The standard choice is the one with the smallest norm, $\hat w = X^+ y$ with the Moore–Penrose pseudoinverse $X^+$, which is what numpy.linalg.lstsq returns. The choice is a preference, not something the data decide, and the regularization note makes the preference explicit.

Numerical conditioning. Even when $X^\top X$ is invertible in exact arithmetic, forming it squares the condition number: $\kappa(X^\top X) = \kappa(X)^2$. Double precision carries about 16 significant digits, and a linear solve loses roughly $\log_{10}\kappa$ of them. For polynomial features $1, t, \dots, t^{10}$ on 50 points of $[0, 1]$, $\kappa(X) = 2.0 \cdot 10^7$ and $\kappa(X^\top X) = 4.1 \cdot 10^{14}$. Solving the normal equations recovers coefficients that are all equal to 1 with an error of $4 \cdot 10^{-3}$; an SVD-based solver applied to $X$ directly, as lstsq does, gets $4 \cdot 10^{-10}$ (Trefethen and Bau, 1997, lectures 11 and 19). The formula $(X^\top X)^{-1}X^\top y$ is for derivations; numerical code factorizes $X$. Centring and scaling the features, as in the features note, also lowers the condition number, and the next note shows that the same number governs the speed of gradient descent.

Measuring a regression#

The mean squared error is in squared units of the target; its root, the RMSE, is in the target’s units. The coefficient of determination compares the fit with the best constant, the mean $\bar y$:

$$R^2 = 1 - \frac{\sum_i (y_i - \hat y_i)^2}{\sum_i (y_i - \bar y)^2} .$$

On the training data of a model with an intercept, the Pythagorean split of the previous sections puts $R^2$ between 0 and 1, the share of the variation of $y$ that the fit explains. In the four-point example $\bar y = 2.75$, the total sum of squares is $8.75$, and $R^2 = 1 - 2.7/8.75 = 0.691$.

On the diabetes data bundled with scikit-learn (442 patients, 10 centred and scaled measurements, the target a measure of disease progression after a year), least squares has a training $R^2$ of $0.518$ and a 5-fold cross-validated $R^2$ of $0.489$. Nearest neighbours on standardized features reach at best $0.454$ (at $k = 20$; $0.382$, $0.444$ and $0.428$ for $k = 5$, 10 and 40). Eleven coefficients fitted to all the data beat the local vote here, because the relationship is close to additive and 442 points are few for neighbourhoods in ten dimensions.

$R^2$ has three traps:

• Training $R^2$ never decreases with features. Adding a column can only keep or enlarge the column space, so the projection’s residual never grows. With 50 training points, one useful feature and $k$ columns of pure noise, the training $R^2$ goes $0.572$, $0.649$, $0.747$, $0.921$ for $k = 0, 10, 20, 40$, while on 50 fresh points it goes $0.325$, $0.154$, $-0.158$, $-0.71$. A negative $R^2$ means worse than predicting the mean, and it only appears on data the model did not see. Report $R^2$ on validation data, as the validation note requires of every metric.
• $R^2$ depends on the spread of $x$. For $y = x + \varepsilon$ with $\sigma = 1$, the same model with the same noise has $R^2 = (1/12)/(1/12 + 1) = 0.077$ when $x$ is uniform on $[0, 1]$ and $0.893$ when it is uniform on $[0, 10]$; the RMSE is 1 in both cases. $R^2$ compares with the mean, and the mean is a worse predictor when $y$ varies more.
• A high $R^2$ does not mean a correct model. A straight line through points on a parabola can have $R^2$ near 1 and residuals that are all positive at the ends and negative in the middle. Look at the residuals.

Squared error and outliers#

Squared error estimates the conditional mean, and a mean follows its largest values. Other losses estimate other targets (other losses, other targets): the best constant under absolute error is the median. For the five values $1, 2, 4, 8, 16$ the best constant under squared error is their mean $6.2$, under absolute error their median $4$. One wrong value can therefore move a least-squares fit arbitrarily far, while it moves the median-based fit only until it crosses the other points.

Figure

One outlier and three losses

The outlier in the figure has a least-squares residual of $6.07$; squared, $36.9$, it is larger than the squared residuals of the other nine points together. Least squares tilts the line towards it, from slope $0.48$ without the point to $0.82$ with it. Two losses limit the influence of large residuals:

• Absolute error, $\rho(r) = \lvert r\rvert$, least absolute deviations. Each point pulls with the same force whatever its distance. The fitted line here is $1.27 + 0.47x$; the minimizer is not always unique, and the loss is not differentiable at 0.
• The Huber loss (Huber, 1964): $\rho_\delta(r) = r^2/2$ for $\lvert r\rvert \le \delta$ and $\delta(\lvert r\rvert - \delta/2)$ beyond. Quadratic for small residuals, so it keeps the efficiency of least squares on clean Gaussian data; linear for large ones, so an outlier’s pull is capped at $\delta$. With $\delta = 1$ the line is $0.99 + 0.53x$.

The cap is on the derivative with respect to the residual. The gradient with respect to the weights multiplies that derivative by $x_i$, so a point with an extreme input, a high-leverage point, still pulls hard. Ten points on the line $y = x$ and one gross error at $(100, 0)$ flatten all three fits alike: the slopes of least squares, Huber and least absolute deviations are all $-0.04$. These losses protect against outliers in $y$, not in $x$.

Neither has a closed form. Both can be fitted by repeatedly solving weighted least squares with weights that shrink for large residuals, which is how the figure computes them, or by the general method of the next note.

What changed#

	Nearest neighbours	Least squares
Model	the average of nearby targets	a global linear function with $d$ coefficients
Training	store the data	solve $X^\top X w = X^\top y$, $O(nd^2 + d^3)$
Prediction	search the training set	one dot product, $O(d)$
Diabetes, cross-validated $R^2$	$0.454$ at best ($k = 20$)	$0.489$
Uncertainty	none built in	$\operatorname{Var}\hat w = \sigma^2 (X^\top X)^{-1}$, standard errors
Interpretation	none	one weight per feature, if the features are not collinear
Weak point	many dimensions, scaling	nonlinear relationships, collinearity, outliers, $d$ close to $n$

Limitations#

• Linear in the chosen features. The model sees interactions and curvature only if someone puts them into $X$; the trees section and the networks section learn them instead.
• Unstable with correlated features or $d$ close to $n$. The coefficients’ variance blows up along nearly constant directions; with $d \ge n$ and $\operatorname{rank} X = n$ the training error is zero, so the fit on the training data is no evidence of generalization. Which of the many interpolating solutions is chosen (the minimum-norm one, for example) and the geometry of the data still decide whether its predictions are useful, as the double-descent caveat shows.
• Sensitive to outliers. A single gross error can turn the line, because the squared loss weighs residuals by their size; robust losses help against outliers in $y$, not against extreme inputs.
• The variance formula assumes uncorrelated, equal-variance noise. With heteroscedastic or correlated errors (time series, grouped data) the coefficients stay unbiased but the standard errors are wrong.
• The closed form costs $O(nd^2 + d^3)$ and exists only for squared error.

What you should remember#

• Least squares minimizes $\lVert y - Xw\rVert^2$; the normal equations $X^\top X\hat w = X^\top y$ have the unique solution $(X^\top X)^{-1}X^\top y$ when the columns of $X$ are independent.
• Geometrically, $\hat y$ is the orthogonal projection of $y$ onto the column space; the residual is orthogonal to every feature, and with an intercept the residuals sum to zero.
• Under Gaussian noise least squares is maximum likelihood, and $\hat w \sim N\big(w, \sigma^2 (X^\top X)^{-1}\big)$; estimate $\sigma^2$ by $\mathrm{RSS}/(n - d)$.
• Nearly collinear features make individual coefficients meaningless while the predictions stay stable; this is the opening for regularization.
• Never solve the normal equations numerically for ill-conditioned features: forming $X^\top X$ squares the condition number. Use QR or SVD (lstsq).
• $R^2$ compares with predicting the mean; on training data it never decreases as features are added, so report it on validation data.
• Squared error estimates the conditional mean and follows outliers; absolute error estimates the median, and the Huber loss interpolates between the two.

Code#

The worked examples, the sampling distribution of the coefficients, the diabetes comparison with nearest neighbours, the traps of $R^2$, collinearity, conditioning and robust fits:

import numpy as np
from scipy import optimize
from sklearn.datasets import load_diabetes
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import KFold, cross_val_score
from sklearn.neighbors import KNeighborsRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

rng = np.random.default_rng(0)
r2 = lambda y, p: 1 - ((y - p) ** 2).sum() / ((y - y.mean()) ** 2).sum()

# the worked example: four points, the normal equations, residuals orthogonal to the columns
x = np.array([0., 1, 2, 3]); y = np.array([1., 3, 2, 5])
X = np.column_stack([np.ones(4), x])
w = np.linalg.solve(X.T @ X, X.T @ y)
r = y - X @ w
H = X @ np.linalg.inv(X.T @ X) @ X.T
print(w, r.round(2), (X.T @ r).round(12) + 0, round(r @ r, 2), round(r2(y, X @ w), 3), np.diag(H).round(2))
# [1.1 1.1] [-0.1  0.8 -1.3  0.6] [0. 0.] 2.7 0.691 [0.7 0.3 0.3 0.7]
s2 = r @ r / (4 - 2)                                  # unbiased noise variance, n - d = 2
print(round(s2, 2), np.sqrt(s2 * np.diag(np.linalg.inv(X.T @ X))).round(3))   # 1.35 [0.972 0.52 ]

# the sampling distribution of w-hat: 30 points on [0, 1], w = (1, 2), sigma = 0.5, 20,000 noise draws
Xs = np.column_stack([np.ones(30), np.linspace(0, 1, 30)])
W = np.array([np.linalg.lstsq(Xs, Xs @ [1, 2] + 0.5 * rng.normal(size=30), rcond=None)[0] for _ in range(20_000)])
print(W.mean(0).round(3), W.std(0).round(3), (0.5 * np.sqrt(np.diag(np.linalg.inv(Xs.T @ Xs)))).round(3))
# [1.001 2.   ] [0.177 0.303] [0.178 0.306]: simulated means, simulated and theoretical standard deviations

# diabetes (442 patients, 10 features): the closed form is what scikit-learn computes; linear model against kNN
Xd, yd = load_diabetes(return_X_y=True)
A = np.column_stack([np.ones(len(yd)), Xd])
wd = np.linalg.lstsq(A, yd, rcond=None)[0]
lr = LinearRegression().fit(Xd, yd)
print(f'{np.abs(wd - np.r_[lr.intercept_, lr.coef_]).max():.0e}', round(r2(yd, A @ wd), 3))   # 2e-12 0.518
cv = KFold(5, shuffle=True, random_state=0)
print(round(cross_val_score(LinearRegression(), Xd, yd, cv=cv, scoring='r2').mean(), 3),
      [round(float(cross_val_score(make_pipeline(StandardScaler(), KNeighborsRegressor(k)), Xd, yd, cv=cv, scoring="r2").mean()), 3) for k in (5, 10, 20, 40)])
# 0.489 [0.382, 0.444, 0.454, 0.428]

# R² on the training data never decreases with useless features; on fresh data it collapses (50 + 50 points)
z = rng.normal(size=(100, 1)); t = z[:, 0] + rng.normal(size=100)
for k in (0, 10, 20, 40):
    Z = np.column_stack([np.ones(100), z, rng.normal(size=(100, k))])
    v = np.linalg.lstsq(Z[:50], t[:50], rcond=None)[0]
    print(k, round(r2(t[:50], Z[:50] @ v), 3), round(r2(t[50:], Z[50:] @ v), 3))
# 0 0.572 0.325 | 10 0.649 0.154 | 20 0.747 -0.158 | 40 0.921 -0.71

# two nearly identical features: the coefficients swing, their sum and the predictions do not
x1 = rng.normal(size=100); x2 = x1 + 0.01 * rng.normal(size=100); X2 = np.column_stack([x1, x2])
print(round(np.corrcoef(x1, x2)[0, 1], 5), f'{np.linalg.cond(X2.T @ X2):.1e}')   # 0.99994 3.5e+04
for _ in range(3):
    v = np.linalg.lstsq(X2, x1 + x2 + rng.normal(size=100), rcond=None)[0]
    print(v.round(2), round(v.sum(), 2))
# [-5.    6.88] 1.88 | [11.88 -9.9 ] 1.98 | [-8.83 10.94] 2.11
print(np.sqrt(np.diag(np.linalg.inv(X2.T @ X2))).round(2), round(1 / np.sqrt(x1 @ x1), 2))   # [9.74 9.73] 0.1: standard errors for sigma = 1

# conditioning: polynomial features of degree 10 on 50 points; normal equations against QR/SVD (lstsq)
V = np.vander(np.linspace(0, 1, 50), 11, increasing=True); yv = V @ np.ones(11)
print(f'{np.linalg.cond(V):.1e} {np.linalg.cond(V.T @ V):.1e}',
      f'{np.abs(np.linalg.solve(V.T @ V, V.T @ yv) - 1).max():.0e} {np.abs(np.linalg.lstsq(V, yv, rcond=None)[0] - 1).max():.0e}')
# 2.0e+07 4.1e+14 4e-03 4e-10

# one outlier: least squares, Huber (delta = 1) and least absolute deviations
xo = np.arange(10.); yo = np.array([1.3, 1.1, 2.2, 3.0, 2.7, 3.6, 3.5, 4.9, 13.0, 5.3]); Xo = np.column_stack([np.ones(10), xo])
huber = lambda r: np.where(np.abs(r) <= 1, 0.5 * r ** 2, np.abs(r) - 0.5)
fits = {'ls': np.linalg.lstsq(Xo, yo, rcond=None)[0],
        'ls without x=8': np.linalg.lstsq(np.delete(Xo, 8, 0), np.delete(yo, 8), rcond=None)[0],
        'huber': optimize.minimize(lambda v: huber(yo - Xo @ v).sum(), [0, 0], method='BFGS', options={'gtol': 1e-10}).x,
        'lad': optimize.linprog(np.r_[0, 0, np.ones(10)], A_ub=np.block([[Xo, -np.eye(10)], [-Xo, -np.eye(10)]]),
                                b_ub=np.r_[yo, -yo], bounds=[(None, None)] * 2 + [(0, None)] * 10).x[:2]}
for k, v in fits.items(): print(k, v.round(3))
# ls [0.373 0.819] | ls without x=8 [1.107 0.477] | huber [0.986 0.533] | lad [1.267 0.467]

# a gross error at an extreme input: ten points on y = x and one at (100, 0)
xl = np.r_[np.arange(10.), 100]; yl = np.r_[np.arange(10.), 0]; Xl = np.column_stack([np.ones(11), xl])
print(round(float(np.linalg.lstsq(Xl, yl, rcond=None)[0][1]), 2),
      round(float(optimize.minimize(lambda v: huber(yl - Xl @ v).sum(), [0, 1], method='BFGS', options={'gtol': 1e-10}).x[1]), 2),
      round(float(optimize.linprog(np.r_[0, 0, np.ones(11)], A_ub=np.block([[Xl, -np.eye(11)], [-Xl, -np.eye(11)]]),
                                   b_ub=np.r_[yl, -yl], bounds=[(None, None)] * 2 + [(0, None)] * 11).x[1]), 2))
# -0.04 -0.04 -0.04: the slopes of least squares, Huber and least absolute deviations

Next#

The normal equations need the whole matrix $X^\top X$ and a $d \times d$ solve, which is out of reach for a million features, and they exist only for squared error: the Huber and absolute losses of the last section, and the logistic loss of the classification section, have no closed form at all. The next note replaces the formula by an iteration that only needs the gradient, $2X^\top(Xw - y)$ here, and shows that its speed is set by the same condition number that made the normal equations inaccurate.

Least squares reappears across the courses: GloVe word vectors are fitted by weighted least squares on log co-occurrence counts, the least-squares GAN replaces the discriminator’s cross-entropy by squared error, and alternating least squares fits the matrix factorization of recommender systems.

References#

• E. Sokolov et al. Машинное обучение (HSE), lecture notes, lecture 2: linear models, loss functions and quality measures for regression and lecture 3: overfitting, training linear regression.
• T. Hastie, R. Tibshirani, J. Friedman. The Elements of Statistical Learning, 2nd ed. Springer, 2009. Section 3.2 (least squares, its geometry and the sampling distribution of the coefficients), Section 3.2.2 (the Gauss–Markov theorem).
• C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006. Section 3.1.1 (maximum likelihood and least squares), Section 3.1.2 (the geometry of least squares).
• M. P. Deisenroth, A. A. Faisal, C. S. Ong. Mathematics for Machine Learning. Cambridge University Press, 2020. Chapter 9 (linear regression), Section 3.8 (orthogonal projections).
• L. N. Trefethen, D. Bau III. Numerical Linear Algebra. SIAM, 1997. Lecture 11 (least squares problems), Lecture 19 (stability of least squares algorithms: the normal equations against QR and SVD).
• P. J. Huber. Robust estimation of a location parameter. Annals of Mathematical Statistics, 35(1), 1964.