Linear Regression and Least Squares: Reading the Normal Equations as an Orthogonal Projection

0. Key points

• Training a linear regression model means searching for a $\boldsymbol{w}$ that satisfies $\boldsymbol{y} \approx X\boldsymbol{w}$ as well as possible for an observed vector $\boldsymbol{y}$. The intercept $b$ can be absorbed into the weights by appending a column of all ones to the design matrix $X$.
• Minimising the sum of squared errors $L(\boldsymbol{w}) = \lVert \boldsymbol{y} - X\boldsymbol{w}\rVert^2$ is exactly equivalent to solving the normal equations $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$, and no differentiation is needed to see this (Theorem 3.3). That learning collapses into a system of linear equations is the whole point of the model.
• The prediction $X\hat{\boldsymbol{w}}$ produced by a solution $\hat{\boldsymbol{w}}$ of the normal equations is nothing other than the orthogonal projection of $\boldsymbol{y}$ onto the column space $\operatorname{Im} X$ (Theorem 4.2). The residual is orthogonal to every explanatory variable.
• The normal equations always have a solution, whatever $X$ and $\boldsymbol{y}$ may be. The solution is unique precisely when $\operatorname{rank} X = p$ (the columns are linearly independent), and then $\hat{\boldsymbol{w}} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}$. Even when the rank drops, the prediction $X\hat{\boldsymbol{w}}$ remains uniquely determined.
• In an implementation, never form $(X^{\mathsf{T}}X)^{-1}$ explicitly. Since $\kappa_2(X^{\mathsf{T}}X) = \kappa_2(X)^2$ squares the condition number, going through a QR decomposition or a singular value decomposition is safer.

1. Motivation

1.1. The problem of having too much information

As one learns in school, there is exactly one line through two distinct points in the plane. So which line “passes through” 100 points contaminated by measurement error? The answer is: none. There are only two unknowns, the slope and the intercept, but there are 100 equations. Such an overdetermined system — a system of linear equations with more equations than unknowns — has no solution in the ordinary sense.

This situation was entirely routine in the astronomy and geodesy of the early nineteenth century. An orbit is described by six elements, yet observations are made dozens of times. Practitioners of the day resorted to ad hoc procedures: selecting a few observations and writing down equations, or grouping observations and averaging them. What was wanted was a procedure that uses all the observations impartially and still returns a unique answer.

That procedure is the method of least squares. On 1 January 1801 Giuseppe Piazzi discovered the asteroid Ceres, but after 41 days of observation it moved towards the Sun and was lost. Gauss, then 24, computed an orbit from the scanty data and predicted where it would reappear; at the end of the same year astronomers recovered Ceres exactly where he had said it would be. The episode made the power of least squares widely known. The method itself was published in 1805 by Legendre, as “méthode des moindres carrés” in an appendix to his Nouvelles méthodes pour la détermination des orbites des comètes, whereupon Gauss claimed in his Theoria Motus of 1809 that he had been using it since 1795 — the origin of a famous priority dispute.

1.2. Why squares?

The sum of squares is not the only way to measure error. The sum of absolute values $\sum_i |y_i - \hat{y}_i|$ or the maximum error $\max_i |y_i - \hat{y}_i|$ would do just as well. There are nevertheless clear reasons why the sum of squares became standard.

1. It is differentiable. The absolute value is not differentiable at the origin, whereas the square is smooth everywhere and, moreover, a quadratic in $\boldsymbol{w}$. The condition for a stationary point of a quadratic is a linear equation.
2. It connects directly with geometry. The sum of squares is the square of the Euclidean norm, that is, the distance coming from an inner product. Hence “minimise the error” means the same as “drop a perpendicular onto a subspace”, and the Pythagorean theorem applies verbatim (Theorem 4.2).
3. It has a probabilistic meaning. If the observation noise is assumed to be independent and identically normally distributed, maximising the log-likelihood coincides exactly with minimising the sum of squared errors. This was Gauss’s own justification; it is treated in The Role of Probability Theory and Bayesian Statistics at maximum likelihood under Gaussian noise is least squares（Proposition 3.3）［確率論とベイズ統計］.

There is a weakness on the other side. Squaring emphasises large errors still further, so a single outlier can move the solution a long way. Remedies for this weakness (robust regression, regularisation) are touched on in the later sections and in the exercises.

1.3. Where this sits in machine learning

As we saw in Why Does Machine Learning Need Mathematics?, machine learning algorithms become suddenly transparent once written in the language of mathematics. Among supervised learning models, linear regression is essentially the only one whose learning problem admits a closed-form solution. For logistic regression or neural networks the loss ceases to be a quadratic in $\boldsymbol{w}$, the closed form disappears, and one falls back on iterative methods such as Gradient Descent. It is precisely because a closed form exists here that it is worth understanding completely what is going on. The triple that appears here — orthogonal projection, rank condition, condition number — carries over unchanged to Principal Component Analysis and to the iteratively reweighted least squares of Logistic Regression.

2. Preliminaries: arranging the data in a matrix

The data consist of $n$ observations $(\boldsymbol{x}_i, y_i)$ for $i = 1, \dots, n$. Here $\boldsymbol{x}_i \in \mathbb{R}^d$ is the vector of explanatory variables (features) and $y_i \in \mathbb{R}$ is the response. The relationship we wish to predict has the form

$$\hat{y} = w_1 x_1 + w_2 x_2 + \cdots + w_d x_d + b$$

where $b$ is the intercept (bias term). To handle these $d+1$ unknowns together we set up the following notation.

Definition 2.1（Design matrix and the linear regression model）

Given $n$ observations $(\boldsymbol{x}_i, y_i)$ with $\boldsymbol{x}_i = (x_{i1}, \dots, x_{id})^{\mathsf{T}} \in \mathbb{R}^d$ and $y_i \in \mathbb{R}$, put $p = d + 1$ and set

$$X = \begin{pmatrix} 1 & x_{11} & \cdots & x_{1d} \\ 1 & x_{21} & \cdots & x_{2d} \\ \vdots & \vdots & & \vdots \\ 1 & x_{n1} & \cdots & x_{nd} \end{pmatrix} \in \mathbb{R}^{n \times p}, \qquad \boldsymbol{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} \in \mathbb{R}^{n}.$$

The matrix $X$ is called the design matrix. Its first column is the vector $\boldsymbol{1}$ all of whose entries are $1$. Writing the parameters as $\boldsymbol{w} = (b, w_1, \dots, w_d)^{\mathsf{T}} \in \mathbb{R}^{p}$, the $n$ predicted values are collected into $X\boldsymbol{w} \in \mathbb{R}^n$. This is the linear regression model.

By placing $\boldsymbol{1}$ in the first column, the intercept becomes “the weight attached to the constant feature $1$”. From now on the intercept needs no special treatment. The $i$-th row of $X$ is the $i$-th data point, and the $j$-th column is the $j$-th feature recorded across all data. We shall use the correspondence rows are data, columns are features throughout.

Let us make the meaning of “linear” precise. The “linear” in linear regression means linear in the parameters $\boldsymbol{w}$, not in the explanatory variables. This distinction widens the range of applications enormously.

Definition 2.2（Extension by a feature map）

Fix an arbitrary map $\varphi \colon \mathbb{R}^d \to \mathbb{R}^{p}$ and let $\Phi \in \mathbb{R}^{n \times p}$ be the matrix obtained by replacing the $i$-th row of the design matrix with $\varphi(\boldsymbol{x}_i)^{\mathsf{T}}$. The model $\hat{y} = \boldsymbol{w}^{\mathsf{T}}\varphi(\boldsymbol{x})$ is called the linear model with feature map $\varphi$. However nonlinear $\varphi$ may be, the model is linear in $\boldsymbol{w}$.

For instance, taking $d = 1$ and $\varphi(x) = (1, x, x^2, x^3)^{\mathsf{T}}$ gives a cubic polynomial fit, yet the model is still linear in the unknown $\boldsymbol{w}$. Everything we build below therefore applies directly to polynomial regression, trigonometric expansions and radial basis function models. A concrete numerical instance appears in Example 6.2.

Note

Throughout, $\mathbb{R}^n$ carries the standard inner product $\langle \boldsymbol{u}, \boldsymbol{v}\rangle = \boldsymbol{u}^{\mathsf{T}}\boldsymbol{v} = \sum_{i} u_i v_i$ and the norm $\lVert \boldsymbol{u}\rVert = \sqrt{\langle \boldsymbol{u}, \boldsymbol{u}\rangle}$. For a matrix $A$ we write $\operatorname{Im} A = \{A\boldsymbol{v} : \boldsymbol{v}\}$ (the column space, or image) and $\operatorname{Ker} A = \{\boldsymbol{v} : A\boldsymbol{v} = \boldsymbol{0}\}$ (the kernel). From linear algebra we take as known the rank–nullity theorem $\operatorname{rank} A + \dim \operatorname{Ker} A = (\text{number of columns})$, the equality of row rank and column rank $\operatorname{rank} A = \operatorname{rank} A^{\mathsf{T}}$, and the orthogonal decomposition $\mathbb{R}^n = W \oplus W^{\perp}$ for a subspace $W \subseteq \mathbb{R}^n$. For the first two see the rank–nullity theorem（Theorem 8.3）［Matrices and Linear Systems］ and the equality of row and column rank（Theorem 8.2）［Matrices and Linear Systems］ in Matrices and Linear Systems; for the last, see the orthogonal decomposition theorem（Theorem 7.1）［内積空間とグラム・シュミット直交化］ in Inner Product Spaces and Gram–Schmidt Orthogonalisation.

3. The least squares problem and the normal equations

3.1. Formulating the problem

When $n > p$, there is in general no $\boldsymbol{w}$ satisfying $X\boldsymbol{w} = \boldsymbol{y}$ exactly. We therefore give up on exactness and choose the $\boldsymbol{w}$ that makes the discrepancy as small as possible.

Definition 3.1（The least squares problem）

For $X \in \mathbb{R}^{n \times p}$ and $\boldsymbol{y} \in \mathbb{R}^n$, define the residual vector $\boldsymbol{r}(\boldsymbol{w}) = \boldsymbol{y} - X\boldsymbol{w}$ and the sum of squared errors (residual sum of squares)

$$L(\boldsymbol{w}) = \lVert \boldsymbol{y} - X\boldsymbol{w}\rVert^2 = \sum_{i=1}^{n} \bigl(y_i - (X\boldsymbol{w})_i\bigr)^2 .$$

A vector $\hat{\boldsymbol{w}}$ minimising $L$ over all of $\mathbb{R}^p$, that is, one satisfying

$$\hat{\boldsymbol{w}} \in \operatorname*{arg\,min}_{\boldsymbol{w} \in \mathbb{R}^{p}} \lVert \boldsymbol{y} - X\boldsymbol{w}\rVert^2 ,$$

is called a least squares solution.

$L$ is a quadratic in each component of $\boldsymbol{w}$ and is bounded below (always $L \ge 0$). As we shall see, the minimum is always attained. The minimiser need not be unique, however. Capturing this “exists but need not be unique” structure precisely is the goal of this section and the next.

3.2. The gradient of a quadratic form

We first prove, independently, the differentiation formula we shall need. The partial derivatives here are those of multivariable calculus; if necessary, review the definition of a partial derivative（Definition 3.1）［多変数関数の微分と偏微分］ in Differentiation of Multivariable Functions and Partial Derivatives.

Lemma 3.2（Gradient of a quadratic form）

Let $A \in \mathbb{R}^{p \times p}$ be symmetric ($A^{\mathsf{T}} = A$), let $\boldsymbol{b} \in \mathbb{R}^p$ and $c \in \mathbb{R}$, and define $f \colon \mathbb{R}^p \to \mathbb{R}$ by

$$f(\boldsymbol{w}) = \boldsymbol{w}^{\mathsf{T}} A \boldsymbol{w} - 2\boldsymbol{b}^{\mathsf{T}}\boldsymbol{w} + c .$$

Then $f$ is of class $C^{\infty}$, and its gradient and Hessian are

$$\nabla f(\boldsymbol{w}) = 2A\boldsymbol{w} - 2\boldsymbol{b}, \qquad \nabla^2 f(\boldsymbol{w}) = 2A .$$

Proof（Lemma 3.2）

In components, $f(\boldsymbol{w}) = \sum_{i=1}^{p}\sum_{j=1}^{p} a_{ij} w_i w_j - 2\sum_{i=1}^{p} b_i w_i + c$, a polynomial in $w_1, \dots, w_p$, hence of class $C^{\infty}$. Differentiate with respect to the $k$-th component. The product $w_i w_j$ involves $w_k$ when $i = k$ and when $j = k$; the term $a_{kk}w_k^2$ with $i = j = k$ differentiates to $2a_{kk}w_k$, which agrees with the sum of the two counts. Therefore

$$\frac{\partial f}{\partial w_k}(\boldsymbol{w}) = \sum_{j=1}^{p} a_{kj} w_j + \sum_{i=1}^{p} a_{ik} w_i - 2 b_k = (A\boldsymbol{w})_k + (A^{\mathsf{T}}\boldsymbol{w})_k - 2b_k .$$

By the hypothesis $A^{\mathsf{T}} = A$ the right-hand side equals $2(A\boldsymbol{w})_k - 2b_k$, which gives the first formula. Differentiating once more with respect to $w_l$ gives $\partial^2 f / \partial w_l \partial w_k = 2a_{kl}$, so the Hessian is $2A$.

∎

3.3. The normal equations

We come to the central theorem. Note that it can be proved without any differentiation at all. Because the proof needs nothing but an algebraic identity, it simultaneously establishes that we have a minimum, and not merely a stationary point.

Theorem 3.3（The normal equations）

Let $X \in \mathbb{R}^{n\times p}$ and $\boldsymbol{y} \in \mathbb{R}^n$, and set $L(\boldsymbol{w}) = \lVert \boldsymbol{y} - X\boldsymbol{w}\rVert^2$. For $\hat{\boldsymbol{w}} \in \mathbb{R}^p$ the following two conditions are equivalent.

1. $\hat{\boldsymbol{w}}$ minimises $L$ on $\mathbb{R}^p$, that is, $L(\hat{\boldsymbol{w}}) \le L(\boldsymbol{w})$ for every $\boldsymbol{w} \in \mathbb{R}^p$.
2. $\hat{\boldsymbol{w}}$ satisfies the normal equations $X^{\mathsf{T}}X\hat{\boldsymbol{w}} = X^{\mathsf{T}}\boldsymbol{y}$.

Moreover, in that case the identity

$$L(\boldsymbol{w}) = L(\hat{\boldsymbol{w}}) + \lVert X(\boldsymbol{w} - \hat{\boldsymbol{w}})\rVert^2$$

holds for every $\boldsymbol{w} \in \mathbb{R}^p$, and the set of all least squares solutions coincides with $\hat{\boldsymbol{w}} + \operatorname{Ker} X = \{\hat{\boldsymbol{w}} + \boldsymbol{v} : X\boldsymbol{v} = \boldsymbol{0}\}$.

Proof（Theorem 3.3）

For arbitrary $\hat{\boldsymbol{w}}, \boldsymbol{v} \in \mathbb{R}^p$, put $\hat{\boldsymbol{r}} = \boldsymbol{y} - X\hat{\boldsymbol{w}}$ and expand $L(\hat{\boldsymbol{w}} + \boldsymbol{v})$. Since the squared norm is an inner product,

$$\begin{aligned} L(\hat{\boldsymbol{w}} + \boldsymbol{v}) &= \lVert \hat{\boldsymbol{r}} - X\boldsymbol{v}\rVert^2 \\ &= \langle \hat{\boldsymbol{r}} - X\boldsymbol{v},\ \hat{\boldsymbol{r}} - X\boldsymbol{v}\rangle \\ &= \lVert \hat{\boldsymbol{r}}\rVert^2 - 2\langle \hat{\boldsymbol{r}},\ X\boldsymbol{v}\rangle + \lVert X\boldsymbol{v}\rVert^2 . \end{aligned}$$

By the transpose rule $\langle \boldsymbol{u}, A\boldsymbol{v}\rangle = \langle A^{\mathsf{T}}\boldsymbol{u}, \boldsymbol{v}\rangle$, the middle term can be written

$$\langle \hat{\boldsymbol{r}},\ X\boldsymbol{v}\rangle = \langle X^{\mathsf{T}}\hat{\boldsymbol{r}},\ \boldsymbol{v}\rangle = \boldsymbol{v}^{\mathsf{T}}\bigl(X^{\mathsf{T}}\boldsymbol{y} - X^{\mathsf{T}}X\hat{\boldsymbol{w}}\bigr) .$$

So, setting $\boldsymbol{g} = X^{\mathsf{T}}\boldsymbol{y} - X^{\mathsf{T}}X\hat{\boldsymbol{w}} = X^{\mathsf{T}}\hat{\boldsymbol{r}}$, we obtain for every $\boldsymbol{v}$

$$L(\hat{\boldsymbol{w}} + \boldsymbol{v}) = L(\hat{\boldsymbol{w}}) - 2\boldsymbol{v}^{\mathsf{T}}\boldsymbol{g} + \lVert X\boldsymbol{v}\rVert^2 .$$

Both implications follow from this single identity.

(2) $\Rightarrow$ (1). If the normal equations hold then $\boldsymbol{g} = \boldsymbol{0}$, so the identity becomes $L(\hat{\boldsymbol{w}} + \boldsymbol{v}) = L(\hat{\boldsymbol{w}}) + \lVert X\boldsymbol{v}\rVert^2$. Since $\lVert X\boldsymbol{v}\rVert^2 \ge 0$, we get $L(\hat{\boldsymbol{w}} + \boldsymbol{v}) \ge L(\hat{\boldsymbol{w}})$ for every $\boldsymbol{v}$; that is, $\hat{\boldsymbol{w}}$ is a minimiser. Substituting $\boldsymbol{w} = \hat{\boldsymbol{w}} + \boldsymbol{v}$ gives exactly the asserted identity. Equality holds if and only if $\lVert X\boldsymbol{v}\rVert = 0$, that is, $X\boldsymbol{v} = \boldsymbol{0}$, so the set of minimisers is $\hat{\boldsymbol{w}} + \operatorname{Ker} X$.

(1) $\Rightarrow$ (2). We prove the contrapositive. Suppose $\boldsymbol{g} \ne \boldsymbol{0}$ and substitute $\boldsymbol{v} = t\boldsymbol{g}$ with $t > 0$:

$$L(\hat{\boldsymbol{w}} + t\boldsymbol{g}) - L(\hat{\boldsymbol{w}}) = -2t\lVert \boldsymbol{g}\rVert^2 + t^2 \lVert X\boldsymbol{g}\rVert^2 = t\bigl(t\lVert X\boldsymbol{g}\rVert^2 - 2\lVert \boldsymbol{g}\rVert^2\bigr) .$$

If $X\boldsymbol{g} = \boldsymbol{0}$ the right-hand side equals $-2t\lVert\boldsymbol{g}\rVert^2 < 0$ for every $t > 0$. If $X\boldsymbol{g} \ne \boldsymbol{0}$, taking for example $t = \lVert \boldsymbol{g}\rVert^2 / \lVert X\boldsymbol{g}\rVert^2 > 0$ makes the right-hand side $-\lVert\boldsymbol{g}\rVert^4/\lVert X\boldsymbol{g}\rVert^2 < 0$. In either case $L(\hat{\boldsymbol{w}} + t\boldsymbol{g}) < L(\hat{\boldsymbol{w}})$, so $\hat{\boldsymbol{w}}$ is not a minimiser. Hence a minimiser must satisfy $\boldsymbol{g} = \boldsymbol{0}$, that is, the normal equations.

∎

The normal equations $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$ form a system of $p$ linear equations. However large $n$ may be — even hundreds of millions of data points — the size of the system to be solved is governed solely by the number of features $p$. What actually happens behind the word “learning” is the solution of a $p \times p$ linear system.

Remark 3.4（Derivation by partial differentiation）

Let us also record the derivation using calculus. Expanding,

$$L(\boldsymbol{w}) = \lVert \boldsymbol{y}\rVert^2 - 2\boldsymbol{y}^{\mathsf{T}}X\boldsymbol{w} + \boldsymbol{w}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{w} ,$$

so applying Lemma 3.2 with $A = X^{\mathsf{T}}X$ (symmetric, since $(X^{\mathsf{T}}X)^{\mathsf{T}} = X^{\mathsf{T}}X$), $\boldsymbol{b} = X^{\mathsf{T}}\boldsymbol{y}$ and $c = \lVert\boldsymbol{y}\rVert^2$ gives

$$\nabla L(\boldsymbol{w}) = 2X^{\mathsf{T}}X\boldsymbol{w} - 2X^{\mathsf{T}}\boldsymbol{y} .$$

The condition $\nabla L(\hat{\boldsymbol{w}}) = \boldsymbol{0}$ is the normal equations. Being a stationary point does not in general imply being a minimum. Here, however, the Hessian is $2X^{\mathsf{T}}X$, and $\boldsymbol{v}^{\mathsf{T}}(X^{\mathsf{T}}X)\boldsymbol{v} = \lVert X\boldsymbol{v}\rVert^2 \ge 0$ for every $\boldsymbol{v}$, so it is positive semidefinite; $L$ is therefore convex and a stationary point is a global minimum. The proof of Theorem 3.3 can be viewed as replacing this convexity argument by a single identity.

4. Geometry: the least squares solution is an orthogonal projection

4.1. The residual is orthogonal to every feature

Rewrite the normal equations as $X^{\mathsf{T}}(\boldsymbol{y} - X\hat{\boldsymbol{w}}) = \boldsymbol{0}$, that is, $X^{\mathsf{T}}\hat{\boldsymbol{r}} = \boldsymbol{0}$. The $j$-th component of $X^{\mathsf{T}}\hat{\boldsymbol{r}}$ is the inner product of the $j$-th column of $X$ with $\hat{\boldsymbol{r}}$, so the equation reads:

the residual vector is orthogonal to every column of $X$, that is, to every feature.

In statistical terms, “the residual contains nothing further that could be explained by a linear combination of the available features”. Least squares is a procedure that stops exactly when everything the explanatory variables can account for has been used up. Let us turn this reading into a theorem.

Lemma 4.1（Orthogonal complement of the column space）

For any $X \in \mathbb{R}^{n\times p}$ we have $(\operatorname{Im} X)^{\perp} = \operatorname{Ker} X^{\mathsf{T}}$.

Proof（Lemma 4.1）

That $\boldsymbol{u} \in (\operatorname{Im}X)^{\perp}$ means $\langle \boldsymbol{u}, X\boldsymbol{v}\rangle = 0$ for every $\boldsymbol{v}\in\mathbb{R}^p$. By the transpose rule $\langle \boldsymbol{u}, X\boldsymbol{v}\rangle = \langle X^{\mathsf{T}}\boldsymbol{u}, \boldsymbol{v}\rangle$, this is equivalent to $\langle X^{\mathsf{T}}\boldsymbol{u}, \boldsymbol{v}\rangle = 0$ for every $\boldsymbol{v}$. Taking in particular $\boldsymbol{v} = X^{\mathsf{T}}\boldsymbol{u}$ gives $\lVert X^{\mathsf{T}}\boldsymbol{u}\rVert^2 = 0$, hence $X^{\mathsf{T}}\boldsymbol{u} = \boldsymbol{0}$. Conversely, if $X^{\mathsf{T}}\boldsymbol{u} = \boldsymbol{0}$ then the inner product vanishes for every $\boldsymbol{v}$. The two conditions are therefore equivalent, and the sets coincide.

∎

Theorem 4.2（Existence of least squares solutions and their projection form）

Let $X \in \mathbb{R}^{n\times p}$ and $\boldsymbol{y}\in\mathbb{R}^n$ be arbitrary. Then the following hold.

1. The normal equations $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$ have at least one solution. In particular a least squares solution always exists.
2. If $\hat{\boldsymbol{w}}$ is any least squares solution, then $X\hat{\boldsymbol{w}}$ equals the orthogonal projection of $\boldsymbol{y}$ onto the subspace $\operatorname{Im}X$. In particular $X\hat{\boldsymbol{w}}$ is uniquely determined, independently of which least squares solution is chosen.
3. $\hat{\boldsymbol{y}} = X\hat{\boldsymbol{w}}$ is the unique point of $\operatorname{Im}X$ closest to $\boldsymbol{y}$: if $\boldsymbol{z} \in \operatorname{Im}X$ and $\boldsymbol{z} \ne \hat{\boldsymbol{y}}$, then $\lVert \boldsymbol{y}-\boldsymbol{z}\rVert > \lVert \boldsymbol{y}-\hat{\boldsymbol{y}}\rVert$.

Proof（Theorem 4.2）

(1). Since $W = \operatorname{Im}X$ is a subspace of $\mathbb{R}^n$, the orthogonal decomposition $\mathbb{R}^n = W \oplus W^{\perp}$ gives a unique decomposition $\boldsymbol{y} = \hat{\boldsymbol{y}} + \boldsymbol{s}$ with $\hat{\boldsymbol{y}} \in W$ and $\boldsymbol{s}\in W^{\perp}$. As $\hat{\boldsymbol{y}} \in \operatorname{Im}X$, there is a $\hat{\boldsymbol{w}}$ with $\hat{\boldsymbol{y}} = X\hat{\boldsymbol{w}}$. By Lemma 4.1, $\boldsymbol{s} \in W^{\perp} = \operatorname{Ker}X^{\mathsf{T}}$, that is, $X^{\mathsf{T}}\boldsymbol{s} = \boldsymbol{0}$. Hence

$$X^{\mathsf{T}}X\hat{\boldsymbol{w}} = X^{\mathsf{T}}\hat{\boldsymbol{y}} = X^{\mathsf{T}}(\boldsymbol{y} - \boldsymbol{s}) = X^{\mathsf{T}}\boldsymbol{y} ,$$

so this $\hat{\boldsymbol{w}}$ solves the normal equations. By Theorem 3.3 it is a least squares solution.

(2). Let $\hat{\boldsymbol{w}}$ be any least squares solution. By Theorem 3.3, $X^{\mathsf{T}}(\boldsymbol{y} - X\hat{\boldsymbol{w}}) = \boldsymbol{0}$, hence by Lemma 4.1 we have $\boldsymbol{y} - X\hat{\boldsymbol{w}} \in W^{\perp}$. On the other hand $X\hat{\boldsymbol{w}} \in W$, so $\boldsymbol{y} = X\hat{\boldsymbol{w}} + (\boldsymbol{y}-X\hat{\boldsymbol{w}})$ is a decomposition along $W \oplus W^{\perp}$. Such a decomposition is unique, so $X\hat{\boldsymbol{w}}$ equals the $\hat{\boldsymbol{y}}$ constructed in (1), namely the orthogonal projection of $\boldsymbol{y}$ onto $W$, and does not depend on the choice of least squares solution.

(3). For $\boldsymbol{z}\in W$ we have $\hat{\boldsymbol{y}} - \boldsymbol{z} \in W$, while $\boldsymbol{y}-\hat{\boldsymbol{y}} \in W^{\perp}$, so the two are orthogonal. Hence by the Pythagorean theorem（Remark 4.4）［内積空間とグラム・シュミット直交化］

$$\lVert \boldsymbol{y}-\boldsymbol{z}\rVert^2 = \lVert (\boldsymbol{y}-\hat{\boldsymbol{y}}) + (\hat{\boldsymbol{y}}-\boldsymbol{z})\rVert^2 = \lVert \boldsymbol{y}-\hat{\boldsymbol{y}}\rVert^2 + \lVert \hat{\boldsymbol{y}}-\boldsymbol{z}\rVert^2 .$$

If $\boldsymbol{z}\ne\hat{\boldsymbol{y}}$ the second term is positive, giving the strict inequality.

∎

The geometry of least squares. The observed vector y generally lies outside the column space Im X; the fitted vector is the foot of the perpendicular dropped from y onto Im X (the orthogonal projection). The residual r is orthogonal to the whole column space.

4.2. The hat matrix and the coefficient of determination

When $\operatorname{rank}X = p$ the projection can be written explicitly as a matrix.

Proposition 4.3（The hat matrix）

Suppose $X \in \mathbb{R}^{n\times p}$ satisfies $\operatorname{rank}X = p$ (so that $X^{\mathsf{T}}X$ is invertible by Corollary 5.2). Put

$$P = X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}} \in \mathbb{R}^{n\times n} .$$

Then the following hold.

1. $P^{\mathsf{T}} = P$ and $P^2 = P$ (symmetric idempotent).
2. For every $\boldsymbol{y}\in\mathbb{R}^n$, $P\boldsymbol{y}$ is the orthogonal projection of $\boldsymbol{y}$ onto $\operatorname{Im}X$, and $P\boldsymbol{y} = X\hat{\boldsymbol{w}}$.
3. $I_n - P$ is the orthogonal projection onto $(\operatorname{Im}X)^{\perp} = \operatorname{Ker}X^{\mathsf{T}}$, and the residual is $\hat{\boldsymbol{r}} = (I_n - P)\boldsymbol{y}$.
4. $\operatorname{tr} P = p$.

Proof（Proposition 4.3）

1. $(X^{\mathsf{T}}X)^{-1}$ is the inverse of a symmetric matrix, hence symmetric (substitute $A^{\mathsf{T}}=A$ into $(A^{-1})^{\mathsf{T}} = (A^{\mathsf{T}})^{-1}$). Therefore $P^{\mathsf{T}} = X\bigl((X^{\mathsf{T}}X)^{-1}\bigr)^{\mathsf{T}}X^{\mathsf{T}} = P$. Also

$$P^2 = X(X^{\mathsf{T}}X)^{-1}\underbrace{X^{\mathsf{T}}X(X^{\mathsf{T}}X)^{-1}}_{= I_p}X^{\mathsf{T}} = X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}} = P .$$

2. $\hat{\boldsymbol{w}} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}$ solves the normal equations (multiply both sides by $X^{\mathsf{T}}X$ to check), so $P\boldsymbol{y} = X\hat{\boldsymbol{w}}$, and by part (2) of Theorem 4.2 this is the orthogonal projection.

3. We have $\hat{\boldsymbol{r}} = \boldsymbol{y}-P\boldsymbol{y} = (I_n-P)\boldsymbol{y}$, and the proof of Theorem 4.2 shows $\hat{\boldsymbol{r}} \in (\operatorname{Im}X)^{\perp}$. Since $(I_n-P)^{\mathsf{T}} = I_n - P$ and $(I_n-P)^2 = I_n - 2P + P^2 = I_n - P$, this matrix is also symmetric idempotent; and for $\boldsymbol{u}\in(\operatorname{Im}X)^{\perp}$ we have $X^{\mathsf{T}}\boldsymbol{u} = \boldsymbol{0}$ by Lemma 4.1, hence $P\boldsymbol{u} = \boldsymbol{0}$ and $(I_n-P)\boldsymbol{u} = \boldsymbol{u}$, so $I_n-P$ is the identity on $(\operatorname{Im}X)^{\perp}$.

4. Using the cyclic property of the trace, $\operatorname{tr}(AB) = \operatorname{tr}(BA)$, with $A = X$ and $B = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}$,

$$\operatorname{tr}P = \operatorname{tr}\bigl((X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}X\bigr) = \operatorname{tr}I_p = p .$$∎

Because $\hat{\boldsymbol{y}} = P\boldsymbol{y}$ “puts a hat on $\boldsymbol{y}$”, $P$ is called the hat matrix. In statistics the diagonal entries $P_{ii}$ measure “how strongly the $i$-th observation pulls its own fitted value” (the leverage), and $\operatorname{tr}P = p$ can be read as ”$p$ degrees of freedom have been spent on $n$ observations”.

Corollary 4.4（Decomposition of the sum of squares and the coefficient of determination）

Assume that the columns of $X$ include the all-ones vector $\boldsymbol{1}$. Put $\hat{\boldsymbol{y}} = X\hat{\boldsymbol{w}}$, $\hat{\boldsymbol{r}} = \boldsymbol{y}-\hat{\boldsymbol{y}}$ and $\bar{y} = \frac{1}{n}\sum_i y_i$. Then the following hold.

1. $\sum_{i=1}^n \hat{r}_i = 0$, and consequently $\frac{1}{n}\sum_i \hat{y}_i = \bar{y}$.
2. The decomposition of the sum of squares

$$\underbrace{\sum_{i=1}^n (y_i-\bar{y})^2}_{\text{total } S_{\mathrm{tot}}} = \underbrace{\sum_{i=1}^n (\hat{y}_i-\bar{y})^2}_{\text{explained}} + \underbrace{\sum_{i=1}^n \hat{r}_i^{\,2}}_{\text{residual } S_{\mathrm{res}}}$$

holds. In particular, when $S_{\mathrm{tot}} \ne 0$ the coefficient of determination $R^2 = 1 - S_{\mathrm{res}}/S_{\mathrm{tot}}$ satisfies $0 \le R^2 \le 1$.

Proof（Corollary 4.4）

1. By Theorem 3.3 we have $X^{\mathsf{T}}\hat{\boldsymbol{r}} = \boldsymbol{0}$, and the component corresponding to $\boldsymbol{1}$ is $\langle \boldsymbol{1}, \hat{\boldsymbol{r}}\rangle = \sum_i \hat{r}_i = 0$. Dividing by $n$ gives $\bar{y} - \frac{1}{n}\sum_i\hat{y}_i = 0$.

2. Decompose $\boldsymbol{y}-\bar{y}\boldsymbol{1} = (\hat{\boldsymbol{y}}-\bar{y}\boldsymbol{1}) + \hat{\boldsymbol{r}}$. Since $\boldsymbol{1} \in \operatorname{Im}X$ and $\hat{\boldsymbol{y}}\in\operatorname{Im}X$, we have $\hat{\boldsymbol{y}}-\bar{y}\boldsymbol{1}\in\operatorname{Im}X$, while $\hat{\boldsymbol{r}} \in (\operatorname{Im}X)^{\perp}$ by Theorem 4.2; the two are therefore orthogonal. The asserted identity follows from the Pythagorean theorem. Finally $R^2 = (\text{explained})/S_{\mathrm{tot}}$, and both terms are non-negative with sum $S_{\mathrm{tot}}$, so $0\le R^2\le 1$.

∎

Example 4.5（Fitting a line to three points, seen as a projection）

Take $n = 3$, $d = 1$ and the data $(x_i, y_i) = (0,1), (1,1), (2,4)$. The design matrix and observation vector are

$$X = \begin{pmatrix} 1 & 0 \\ 1 & 1 \\ 1 & 2\end{pmatrix},\qquad \boldsymbol{y} = \begin{pmatrix} 1 \\ 1 \\ 4\end{pmatrix} .$$

Here $\operatorname{Im}X$ is a two-dimensional plane in $\mathbb{R}^3$, and $\boldsymbol{y}$ does not lie on it (the three points are not collinear). Computing the ingredients of the normal equations,

$$X^{\mathsf{T}}X = \begin{pmatrix} 3 & 3 \\ 3 & 5\end{pmatrix},\qquad X^{\mathsf{T}}\boldsymbol{y} = \begin{pmatrix} 1+1+4 \\ 0\cdot 1 + 1\cdot 1 + 2\cdot 4\end{pmatrix} = \begin{pmatrix} 6 \\ 9\end{pmatrix} .$$

Since $\det(X^{\mathsf{T}}X) = 15-9 = 6 \ne 0$, the inverse exists and

$$\hat{\boldsymbol{w}} = \frac{1}{6}\begin{pmatrix} 5 & -3 \\ -3 & 3\end{pmatrix}\begin{pmatrix} 6 \\ 9\end{pmatrix} = \frac{1}{6}\begin{pmatrix} 30-27 \\ -18+27\end{pmatrix} = \begin{pmatrix} 0.5 \\ 1.5\end{pmatrix} ,$$

that is, $\hat{y} = 0.5 + 1.5x$. The fitted values and residuals are

$$\hat{\boldsymbol{y}} = X\hat{\boldsymbol{w}} = (0.5,\ 2.0,\ 3.5)^{\mathsf{T}},\qquad \hat{\boldsymbol{r}} = \boldsymbol{y}-\hat{\boldsymbol{y}} = (0.5,\ -1.0,\ 0.5)^{\mathsf{T}} .$$

Let us verify orthogonality: $\langle \boldsymbol{1}, \hat{\boldsymbol{r}}\rangle = 0.5-1.0+0.5 = 0$ and $\langle (0,1,2)^{\mathsf{T}}, \hat{\boldsymbol{r}}\rangle = 0\cdot 0.5 + 1\cdot(-1.0) + 2\cdot 0.5 = 0$, so the residual is indeed orthogonal to both columns of $X$. The Pythagorean theorem also checks out:

$$\lVert\boldsymbol{y}\rVert^2 = 1+1+16 = 18,\quad \lVert\hat{\boldsymbol{y}}\rVert^2 = 0.25+4+12.25 = 16.5,\quad \lVert\hat{\boldsymbol{r}}\rVert^2 = 0.25+1+0.25 = 1.5$$

gives $18 = 16.5 + 1.5$. The hat matrix can be written down explicitly:

$$P = X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}} = \frac{1}{6}\begin{pmatrix} 5 & 2 & -1 \\ 2 & 2 & 2 \\ -1 & 2 & 5\end{pmatrix} .$$

Indeed $P\boldsymbol{y} = \frac{1}{6}(5+2-4,\ 2+2+8,\ -1+2+20)^{\mathsf{T}} = (0.5, 2.0, 3.5)^{\mathsf{T}} = \hat{\boldsymbol{y}}$. Symmetry is visible at a glance, and $\operatorname{tr}P = (5+2+5)/6 = 2 = p$ agrees with part 4 of Proposition 4.3.

5. When is the solution unique? The rank condition

By Theorem 4.2 the prediction $\hat{\boldsymbol{y}} = X\hat{\boldsymbol{w}}$ is always unique. Whether the coefficient vector $\hat{\boldsymbol{w}}$ itself is unique is another matter. The last assertion of Theorem 3.3 says the solution set is $\hat{\boldsymbol{w}} + \operatorname{Ker}X$, so uniqueness is equivalent to $\operatorname{Ker}X = \{\boldsymbol{0}\}$. Let us translate this into the language of $X^{\mathsf{T}}X$.

Lemma 5.1（Kernel of the Gram matrix）

For any $X\in\mathbb{R}^{n\times p}$,

$$\operatorname{Ker}(X^{\mathsf{T}}X) = \operatorname{Ker}X, \qquad \operatorname{rank}(X^{\mathsf{T}}X) = \operatorname{rank}X .$$

Proof（Lemma 5.1）

If $X\boldsymbol{v} = \boldsymbol{0}$, multiplying on the left by $X^{\mathsf{T}}$ gives $X^{\mathsf{T}}X\boldsymbol{v} = \boldsymbol{0}$, so $\operatorname{Ker}X \subseteq \operatorname{Ker}(X^{\mathsf{T}}X)$. Conversely, suppose $X^{\mathsf{T}}X\boldsymbol{v} = \boldsymbol{0}$ and multiply on the left by $\boldsymbol{v}^{\mathsf{T}}$:

$$0 = \boldsymbol{v}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{v} = (X\boldsymbol{v})^{\mathsf{T}}(X\boldsymbol{v}) = \lVert X\boldsymbol{v}\rVert^2 .$$

Only the zero vector has zero norm, so $X\boldsymbol{v} = \boldsymbol{0}$, giving $\operatorname{Ker}(X^{\mathsf{T}}X)\subseteq\operatorname{Ker}X$. The two kernels therefore coincide. As for the ranks, both $X$ and $X^{\mathsf{T}}X$ have $p$ columns, so applying the rank–nullity theorem

$$\operatorname{rank}A = p - \dim\operatorname{Ker}A$$

to each and using the equality of the kernel dimensions gives the equality of the ranks.

∎

Corollary 5.2（Uniqueness of the least squares solution）

Let $X\in\mathbb{R}^{n\times p}$ and $\boldsymbol{y}\in\mathbb{R}^n$. The following three conditions are equivalent.

1. The $p$ columns of $X$ are linearly independent, that is, $\operatorname{rank}X = p$.
2. $X^{\mathsf{T}}X$ is invertible.
3. There is exactly one least squares solution.

In that case the least squares solution is

$$\hat{\boldsymbol{w}} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y} .$$

Furthermore, $n \ge p$ is necessary for these conditions to hold.

Proof（Corollary 5.2）

1 $\Leftrightarrow$ 2. $X^{\mathsf{T}}X$ is a $p\times p$ square matrix, so invertibility is equivalent to $\operatorname{rank}(X^{\mathsf{T}}X) = p$. By Lemma 5.1 this is equivalent to $\operatorname{rank}X = p$, which is exactly linear independence of the columns.

1 $\Leftrightarrow$ 3. By Theorem 4.2 a least squares solution always exists, and by Theorem 3.3 the solution set is $\hat{\boldsymbol{w}} + \operatorname{Ker}X$. Hence there is exactly one solution if and only if $\operatorname{Ker}X = \{\boldsymbol{0}\}$, which by the rank–nullity theorem is equivalent to $\operatorname{rank}X = p$.

The formula. Under condition 2, setting $\boldsymbol{w} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}\boldsymbol{y}$ gives $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$, so this solves the normal equations and is therefore a least squares solution by Theorem 3.3. Uniqueness has just been shown.

Necessity of $n\ge p$. Since $\operatorname{rank}X \le \min(n,p)$, the equality $\operatorname{rank}X = p$ forces $p \le n$.

∎

When there are fewer data points than features ($n < p$), the solution is in this sense never unique. The situation of “more parameters than data”, an everyday occurrence in deep learning, already shows its face at the level of linear regression.

Example 5.3（When features are duplicated (multicollinearity)）

Take the same data $(x_i, y_i) = (0,1), (1,1), (2,4)$ as in Example 4.5 and add the pointless feature $x' = 2x$, “twice $x$”. The design matrix becomes

$$X = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 1 & 2 \\ 1 & 2 & 4\end{pmatrix} ,$$

whose third column is exactly twice the second, so $\operatorname{rank}X = 2 < 3 = p$. The kernel is

$$\operatorname{Ker}X = \{\, t(0, 2, -1)^{\mathsf{T}} : t\in\mathbb{R}\,\}$$

(indeed $0\cdot\boldsymbol{1} + 2\boldsymbol{x} - 1\cdot(2\boldsymbol{x}) = \boldsymbol{0}$). Appending a third coordinate $0$ to the $(0.5, 1.5)$ found in Example 4.5 gives the least squares solution $\hat{\boldsymbol{w}}_0 = (0.5, 1.5, 0)^{\mathsf{T}}$, so by Theorem 3.3 the solution set is the whole line

$$\{\,(0.5,\ 1.5 - 2t,\ t)^{\mathsf{T}} : t\in\mathbb{R}\,\} .$$

For example $t = 0.75$ gives $(0.5, 0, 0.75)$ and $t = -1$ gives $(0.5, 3.5, -1)$; the coefficients look completely different. Whether “the effect of $x$ is $1.5$” or “the effect of $x$ is $0$ and that of $x'$ is $0.75$” cannot be decided from the data. The prediction, however, is

$$X(0.5,\ 1.5-2t,\ t)^{\mathsf{T}} = (0.5,\ 0.5 + (1.5-2t) + 2t,\ 0.5 + 2(1.5-2t)+4t)^{\mathsf{T}} = (0.5,\ 2,\ 3.5)^{\mathsf{T}} ,$$

constant in $t$, exactly as part (2) of Theorem 4.2 predicts.

Caution

In practice the troublesome case is not exact linear dependence but near linear dependence. When $\operatorname{rank}X = p$ but the smallest eigenvalue of $X^{\mathsf{T}}X$ is nearly $0$, the solution $\hat{\boldsymbol{w}}$ is unique yet swings wildly under tiny changes in $\boldsymbol{y}$. Including height both in centimetres and in inches as features, or including all levels of a dummy-coded categorical variable together with the intercept column (the so-called dummy variable trap), creates exactly this state. The remedies are to drop the redundant columns or to regularise, as in the ridge regression of exercise Exercise 8.3.

6. Working the computations by hand

6.1. The closed formula for simple regression

Let us write out the case $d = 1$ (a single explanatory variable) in full generality. Here $p = 2$ and

$$X^{\mathsf{T}}X = \begin{pmatrix} n & \sum_i x_i \\ \sum_i x_i & \sum_i x_i^2\end{pmatrix}, \qquad X^{\mathsf{T}}\boldsymbol{y} = \begin{pmatrix} \sum_i y_i \\ \sum_i x_i y_i\end{pmatrix} .$$

Put $\bar{x} = \frac1n\sum_i x_i$, $\bar{y} = \frac1n\sum_i y_i$, $S_{xx} = \sum_i (x_i-\bar{x})^2$ and $S_{xy} = \sum_i (x_i-\bar{x})(y_i-\bar{y})$. Expanding,

$$\det(X^{\mathsf{T}}X) = n\sum_i x_i^2 - \Bigl(\sum_i x_i\Bigr)^2 = n\Bigl(\sum_i x_i^2 - n\bar{x}^2\Bigr) = nS_{xx} ,$$

so if $S_{xx}\ne 0$ (the $x_i$ are not all equal) then Corollary 5.2 applies and

$$\hat{\boldsymbol{w}} = \frac{1}{nS_{xx}}\begin{pmatrix} \sum_i x_i^2 & -\sum_i x_i \\ -\sum_i x_i & n\end{pmatrix}\begin{pmatrix} n\bar{y} \\ \sum_i x_i y_i\end{pmatrix} .$$

The second component is $\bigl(n\sum_i x_iy_i - n\bar{y}\sum_i x_i\bigr)/(nS_{xx}) = (\sum_i x_iy_i - n\bar{x}\bar{y})/S_{xx} = S_{xy}/S_{xx}$. Computing the first component in the same way gives $(\bar{y}\sum_i x_i^2 - \bar{x}\sum_i x_iy_i)/S_{xx}$, which agrees with $\bar{y}S_{xx} - \bar{x}S_{xy} = \bar{y}\sum_i x_i^2 - n\bar{x}^2\bar{y} - \bar{x}\sum_i x_iy_i + n\bar{x}^2\bar{y} = \bar{y}\sum_i x_i^2 - \bar{x}\sum_i x_iy_i$. In summary we recover the formula familiar from school:

$$\hat{w}_1 = \frac{S_{xy}}{S_{xx}}, \qquad \hat{b} = \bar{y} - \hat{w}_1\bar{x} .$$

The second equation says that the regression line always passes through the centroid $(\bar{x},\bar{y})$. It is also a restatement of part 1 of Corollary 4.4 (that the residuals sum to $0$).

Example 6.1（Carrying a five-point simple regression through to the end）

Take the data $(x_i,y_i) = (1,2), (2,3), (3,5), (4,4), (5,6)$. With $n = 5$,

$$\sum_i x_i = 15,\quad \sum_i y_i = 20,\quad \sum_i x_i^2 = 55,\quad \sum_i x_iy_i = 2+6+15+16+30 = 69 ,$$

so $\bar{x} = 3$ and $\bar{y} = 4$. The normal equations read

$$\begin{pmatrix} 5 & 15 \\ 15 & 55\end{pmatrix}\begin{pmatrix} b \\ w_1\end{pmatrix} = \begin{pmatrix} 20 \\ 69\end{pmatrix} .$$

Since $\det = 5\cdot 55 - 15^2 = 275-225 = 50 \ne 0$,

$$\begin{pmatrix} \hat{b} \\ \hat{w}_1\end{pmatrix} = \frac{1}{50}\begin{pmatrix} 55 & -15 \\ -15 & 5\end{pmatrix}\begin{pmatrix} 20 \\ 69\end{pmatrix} = \frac{1}{50}\begin{pmatrix} 1100 - 1035 \\ -300 + 345\end{pmatrix} = \frac{1}{50}\begin{pmatrix} 65 \\ 45\end{pmatrix} = \begin{pmatrix} 1.3 \\ 0.9\end{pmatrix} ,$$

that is, $\hat{y} = 1.3 + 0.9x$. Let us check with the formula as well: $S_{xx} = 4+1+0+1+4 = 10$ and $S_{xy} = (-2)(-2)+(-1)(-1)+0\cdot 1+1\cdot 0+2\cdot 2 = 9$, so $\hat{w}_1 = 9/10 = 0.9$ and $\hat{b} = 4 - 0.9\cdot 3 = 1.3$, in agreement.

The fitted values are $\hat{\boldsymbol{y}} = (2.2,\ 3.1,\ 4.0,\ 4.9,\ 5.8)^{\mathsf{T}}$ and the residuals are

$$\hat{\boldsymbol{r}} = (-0.2,\ -0.1,\ 1.0,\ -0.9,\ 0.2)^{\mathsf{T}} .$$

Verifying the orthogonality conditions, $\sum_i \hat{r}_i = -0.2-0.1+1.0-0.9+0.2 = 0$ and $\sum_i x_i\hat{r}_i = -0.2-0.2+3.0-3.6+1.0 = 0$, both zero. The sums of squares are $S_{\mathrm{res}} = 0.04+0.01+1.00+0.81+0.04 = 1.90$ and $S_{\mathrm{tot}} = 4+1+1+0+4 = 10$, so by Corollary 4.4 the explained variation is $10-1.9 = 8.1$, which matches $S_{xy}^2/S_{xx} = 81/10 = 8.1$. The coefficient of determination is $R^2 = 1 - 1.9/10 = 0.81$.

The same computation in NumPy looks as follows. Note that it solves a linear system rather than forming an inverse.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([2.0, 3.0, 5.0, 4.0, 6.0])
X = np.column_stack([np.ones_like(x), x])        # design matrix (5, 2)

# solve the normal equations X^T X w = X^T y directly
w_ne = np.linalg.solve(X.T @ X, X.T @ y)         # [1.3, 0.9]

# numerically preferable (uses an SVD internally)
w_ls, *_ = np.linalg.lstsq(X, y, rcond=None)     # [1.3, 0.9]

r = y - X @ w_ls
print(X.T @ r)          # very close to [0. 0.] (up to rounding error)
print(1 - r @ r / ((y - y.mean()) ** 2).sum())   # 0.81

6.2. Polynomial fitting

As Definition 2.2 says, merely swapping the feature map turns the procedure into curve fitting. None of the machinery changes.

Example 6.2（Fitting a quadratic polynomial）

Take the data $(x_i, y_i) = (-2,6), (-1,2), (0,2), (1,2), (2,5)$ and $\varphi(x) = (1, x, x^2)^{\mathsf{T}}$, that is, fit $\hat{y} = w_0 + w_1x + w_2x^2$. Here $p = 3$, and the required sums are

$$\sum_i x_i = 0,\quad \sum_i x_i^2 = 10,\quad \sum_i x_i^3 = 0,\quad \sum_i x_i^4 = 34$$

(the odd powers vanish because the $x_i$ are symmetric about $0$), together with

$$\sum_i y_i = 17,\quad \sum_i x_iy_i = -12-2+0+2+10 = -2,\quad \sum_i x_i^2y_i = 24+2+0+2+20 = 48 .$$

The normal equations are therefore

$$\begin{pmatrix} 5 & 0 & 10 \\ 0 & 10 & 0 \\ 10 & 0 & 34 \end{pmatrix} \begin{pmatrix} w_0 \\ w_1 \\ w_2\end{pmatrix} = \begin{pmatrix} 17 \\ -2 \\ 48\end{pmatrix} .$$

The second row immediately gives $\hat{w}_1 = -0.2$. The first and third rows form a system in $w_0$ and $w_2$ alone; doubling the first row and subtracting it from the third gives

$$(34 - 20)w_2 = 48 - 34, \qquad \text{that is}\quad 14w_2 = 14 ,$$

so $\hat{w}_2 = 1$, and then the first row gives $\hat{w}_0 = (17-10)/5 = 1.4$. The fitted curve is

$$\hat{y} = 1.4 - 0.2x + x^2 .$$

The fitted values are $\hat{\boldsymbol{y}} = (5.8,\ 2.6,\ 1.4,\ 2.2,\ 5.0)^{\mathsf{T}}$ and the residuals are $\hat{\boldsymbol{r}} = (0.2,\ -0.6,\ 0.6,\ -0.2,\ 0)^{\mathsf{T}}$. Checking the orthogonality conditions again, all three hold:

$$\begin{aligned} \sum_i \hat{r}_i &= 0.2-0.6+0.6-0.2+0 = 0,\\ \sum_i x_i\hat{r}_i &= -0.4+0.6+0-0.2+0 = 0,\\ \sum_i x_i^2\hat{r}_i &= 0.8-0.6+0-0.2+0 = 0. \end{aligned}$$

The residual sum of squares is $S_{\mathrm{res}} = 0.04+0.36+0.36+0.04+0 = 0.8$, and since $\bar{y} = 3.4$ the total is $S_{\mathrm{tot}} = 6.76+1.96+1.96+1.96+2.56 = 15.2$, so $R^2 = 1-0.8/15.2 \approx 0.947$. We fitted a curve, yet all we solved was a $3\times 3$ linear system. That is the power of being linear in the parameters.

7. Numerical considerations

7.1. The overall picture

Let us gather the results so far onto a single page.

flowchart TD
A["data X, y"] --> B["minimise the sum of squared errors L"]
B --> C["normal equations X^T X w = X^T y"]
C --> D{"is rank X = p?"}
D -- "yes" --> E["unique solution; implement via Cholesky or QR"]
D -- "no" --> F["solution set is w0 + Ker X; choose the minimum-norm one"]
E --> G["the prediction Xw is the orthogonal projection of y onto Im X"]
F --> G

From the least squares problem to the solution. The branch is the rank condition, and the vector of predictions is uniquely determined along either branch.

When $\operatorname{rank}X = p$, the matrix $X^{\mathsf{T}}X$ is symmetric and positive definite. Indeed, by Lemma 5.1, if $\boldsymbol{v}\ne\boldsymbol{0}$ then $X\boldsymbol{v}\ne\boldsymbol{0}$, so $\boldsymbol{v}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{v} = \lVert X\boldsymbol{v}\rVert^2 > 0$. A positive definite symmetric matrix admits a Cholesky decomposition $X^{\mathsf{T}}X = LL^{\mathsf{T}}$ with $L$ lower triangular, so the implementation reduces to: form $X^{\mathsf{T}}X$ and $X^{\mathsf{T}}\boldsymbol{y}$, compute the Cholesky factor, and solve two triangular systems.

Tip

It is tempting to write np.linalg.inv(X.T @ X) @ X.T @ y, but do not. Forming an inverse explicitly is wasteful both in operation count and in accuracy. The standard practice is to solve a linear system (np.linalg.solve) or to use the QR or singular value decomposition discussed next (np.linalg.lstsq).

7.2. The condition number gets squared

The normal equations, clean as they are in theory, have a numerical drawback. Measuring the condition number of a matrix $A$ by $\kappa_2(A) = \sigma_{\max}(A)/\sigma_{\min}(A)$ (the ratio of singular values), one has, for $X$ of full rank,

$$\kappa_2(X^{\mathsf{T}}X) = \kappa_2(X)^2 .$$

The reason is that with the singular value decomposition $X = U\Sigma V^{\mathsf{T}}$ we get $X^{\mathsf{T}}X = V\Sigma^{\mathsf{T}}\Sigma V^{\mathsf{T}}$, so the eigenvalues of $X^{\mathsf{T}}X$ are the squares of the singular values of $X$ (for the diagonalisation of symmetric matrices see orthogonal diagonalisation of real symmetric matrices（Corollary 4.3）［スペクトル定理］ in The Spectral Theorem). The condition number indicates by what factor a relative error in the input is amplified in the solution, so the moment one forms $X^{\mathsf{T}}X$ the amplification factor is squared.

Example 7.1（Centring alone improves the condition number）

For the data of Example 6.1 we had $X^{\mathsf{T}}X = \begin{pmatrix} 5 & 15 \\ 15 & 55\end{pmatrix}$. Being a symmetric $2\times 2$ matrix, its eigenvalues come from the characteristic equation $\lambda^2 - 60\lambda + 50 = 0$ (trace $60$, determinant $50$), giving

$$\lambda = 30 \pm \sqrt{900-50} = 30 \pm \sqrt{850} ,$$

that is, $\lambda_{\max} \approx 59.155$ and $\lambda_{\min} \approx 0.845$. Hence $\kappa_2(X^{\mathsf{T}}X) \approx 70.0$ and $\kappa_2(X) \approx \sqrt{70} \approx 8.4$.

Now centre $x$, replacing the second column by $x_i - \bar{x} = -2,-1,0,1,2$. Its inner product with the first column becomes $\sum_i (x_i-\bar{x}) = 0$, so

$$X_{\mathrm{c}}^{\mathsf{T}}X_{\mathrm{c}} = \begin{pmatrix} 5 & 0 \\ 0 & 10\end{pmatrix}$$

is diagonal and $\kappa_2(X_{\mathrm{c}}^{\mathsf{T}}X_{\mathrm{c}}) = 10/5 = 2$. The data are the same; only the choice of coordinates changed. This is one of the reasons why centring and scaling the features is recommended as preprocessing. (Centring does not change the fitted line itself; see exercise Exercise 8.2.)

Remark 7.2（Solving via the QR decomposition）

When $\operatorname{rank}X = p$, one can factor $X = QR$ with $Q\in\mathbb{R}^{n\times p}$ satisfying $Q^{\mathsf{T}}Q = I_p$ and $R\in\mathbb{R}^{p\times p}$ upper triangular with positive diagonal entries. This is precisely Gram–Schmidt orthogonalisation（Theorem 6.1）［内積空間とグラム・シュミット直交化］. Substituting into the normal equations gives

$$X^{\mathsf{T}}X = R^{\mathsf{T}}Q^{\mathsf{T}}QR = R^{\mathsf{T}}R, \qquad X^{\mathsf{T}}\boldsymbol{y} = R^{\mathsf{T}}Q^{\mathsf{T}}\boldsymbol{y} ,$$

so the normal equations become $R^{\mathsf{T}}R\boldsymbol{w} = R^{\mathsf{T}}Q^{\mathsf{T}}\boldsymbol{y}$. Since $R$ is invertible (positive diagonal), so is $R^{\mathsf{T}}$, and multiplying both sides by $(R^{\mathsf{T}})^{-1}$ turns the system into the upper triangular system

$$R\hat{\boldsymbol{w}} = Q^{\mathsf{T}}\boldsymbol{y} ,$$

which is solved by back substitution. Because $X^{\mathsf{T}}X$ is never formed, the condition number stays at $\kappa_2(X)$. The hat matrix also simplifies to $P = QQ^{\mathsf{T}}$. When rank deficiency is possible, the still more robust approach via the singular value decomposition (the Appendix at the end of this article) is used.

7.3. Beyond this article

Here are the typical situations in which the closed form ceases to be usable, together with the tools that take over.

Situation	What happens	Remedy
$p$ in the tens of thousands or more	factoring a $p\times p$ matrix costs $O(p^3)$	Gradient Descent, conjugate gradients
$\operatorname{rank}X < p$, or nearly dependent columns	indeterminate solution, wild coefficients	ridge regression (exercise Exercise 8.3), pseudoinverse (Appendix)
the response is $0$ or $1$	squared error is unnatural and does not give probabilities	Logistic Regression
one wants to construct the features themselves	choosing $\varphi$ becomes the problem	Principal Component Analysis, Neural Networks and Backpropagation
one wants the uncertainty of the coefficients	a point estimate is not enough	The Role of Probability Theory and Bayesian Statistics

From the statistical side, one can show that under the assumption that the errors $\varepsilon_i$ have mean $0$, variance $\sigma^2$ and are uncorrelated, the least squares estimator has the smallest variance among unbiased linear estimators (the Gauss–Markov theorem). We do not treat it here, since no probabilistic assumptions have been made in this article, but the proof needs nothing beyond the hat matrix and the orthogonality established above.

8. Exercises

Exercise 8.1Easy

Show that if the columns of the design matrix $X$ include $\boldsymbol{1} = (1,\dots,1)^{\mathsf{T}}$, then the residual of a least squares solution satisfies $\sum_i \hat{r}_i = 0$. Then verify with explicit data that this fails for the model without intercept $\hat{y} = w x$ (whose design matrix is the single column $\boldsymbol{x}$).

Solution

First part. By Theorem 3.3, a least squares solution satisfies $X^{\mathsf{T}}(\boldsymbol{y}-X\hat{\boldsymbol{w}}) = X^{\mathsf{T}}\hat{\boldsymbol{r}} = \boldsymbol{0}$. The $j$-th component of $X^{\mathsf{T}}\hat{\boldsymbol{r}}$ is the inner product of the $j$-th column of $X$ with $\hat{\boldsymbol{r}}$. If the $j$-th column is $\boldsymbol{1}$, that component is $\langle\boldsymbol{1},\hat{\boldsymbol{r}}\rangle = \sum_i \hat{r}_i$, which therefore vanishes.

Second part. Take the data $(x_1,y_1) = (1,1)$ and $(x_2,y_2) = (2,0)$. The design matrix is $X = (1, 2)^{\mathsf{T}}$ with $n=2$, $p=1$, and $X^{\mathsf{T}}X = 1^2+2^2 = 5$, $X^{\mathsf{T}}\boldsymbol{y} = 1\cdot 1 + 2\cdot 0 = 1$, so $\hat{w} = 1/5 = 0.2$. The fitted values are $\hat{\boldsymbol{y}} = (0.2, 0.4)^{\mathsf{T}}$ and the residuals $\hat{\boldsymbol{r}} = (0.8, -0.4)^{\mathsf{T}}$, whence $\sum_i \hat{r}_i = 0.4 \ne 0$. Of course $\langle \boldsymbol{x}, \hat{\boldsymbol{r}}\rangle = 1\cdot 0.8 + 2\cdot(-0.4) = 0$ still holds. Orthogonality is to the columns of $X$, and $\boldsymbol{1}$, not being in the model, need not be orthogonal to the residual.

Exercise 8.2Standard

Split $X = (\boldsymbol{1} \ \ \tilde{X})$ with $\tilde{X}\in\mathbb{R}^{n\times d}$, and split the parameters correspondingly as $\boldsymbol{w} = (b, \boldsymbol{u})$ with $\boldsymbol{u}\in\mathbb{R}^d$. Write the column means as $\bar{\boldsymbol{x}} = \frac1n\tilde{X}^{\mathsf{T}}\boldsymbol{1}\in\mathbb{R}^d$ and the centred matrix and vector as $X_{\mathrm{c}} = \tilde{X}-\boldsymbol{1}\bar{\boldsymbol{x}}^{\mathsf{T}}$ and $\boldsymbol{y}_{\mathrm{c}} = \boldsymbol{y}-\bar{y}\boldsymbol{1}$. Show that the normal equations are equivalent to

$$\hat{b} = \bar{y}-\bar{\boldsymbol{x}}^{\mathsf{T}}\hat{\boldsymbol{u}}, \qquad X_{\mathrm{c}}^{\mathsf{T}}X_{\mathrm{c}}\,\hat{\boldsymbol{u}} = X_{\mathrm{c}}^{\mathsf{T}}\boldsymbol{y}_{\mathrm{c}} .$$

Solution

Write the normal equations $X^{\mathsf{T}}X\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$ in block form. Since $\boldsymbol{1}^{\mathsf{T}}\boldsymbol{1} = n$, $\boldsymbol{1}^{\mathsf{T}}\tilde{X} = n\bar{\boldsymbol{x}}^{\mathsf{T}}$ and $\boldsymbol{1}^{\mathsf{T}}\boldsymbol{y} = n\bar{y}$,

$$\begin{pmatrix} n & n\bar{\boldsymbol{x}}^{\mathsf{T}} \\ n\bar{\boldsymbol{x}} & \tilde{X}^{\mathsf{T}}\tilde{X}\end{pmatrix} \begin{pmatrix} b \\ \boldsymbol{u}\end{pmatrix} = \begin{pmatrix} n\bar{y} \\ \tilde{X}^{\mathsf{T}}\boldsymbol{y}\end{pmatrix} .$$

The first block row reads $nb + n\bar{\boldsymbol{x}}^{\mathsf{T}}\boldsymbol{u} = n\bar{y}$, that is, $b = \bar{y}-\bar{\boldsymbol{x}}^{\mathsf{T}}\boldsymbol{u}$. Substituting this into the second block row $n\bar{\boldsymbol{x}}b + \tilde{X}^{\mathsf{T}}\tilde{X}\boldsymbol{u} = \tilde{X}^{\mathsf{T}}\boldsymbol{y}$ gives

$$\bigl(\tilde{X}^{\mathsf{T}}\tilde{X} - n\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}}\bigr)\boldsymbol{u} = \tilde{X}^{\mathsf{T}}\boldsymbol{y} - n\bar{\boldsymbol{x}}\bar{y} .$$

It remains to check that the coefficients on the left and right are the centred quantities. From $\boldsymbol{1}^{\mathsf{T}}\boldsymbol{1} = n$ and $\boldsymbol{1}^{\mathsf{T}}\tilde{X} = n\bar{\boldsymbol{x}}^{\mathsf{T}}$,

$$\begin{aligned} X_{\mathrm{c}}^{\mathsf{T}}X_{\mathrm{c}} &= (\tilde{X}-\boldsymbol{1}\bar{\boldsymbol{x}}^{\mathsf{T}})^{\mathsf{T}}(\tilde{X}-\boldsymbol{1}\bar{\boldsymbol{x}}^{\mathsf{T}}) \\ &= \tilde{X}^{\mathsf{T}}\tilde{X} - \tilde{X}^{\mathsf{T}}\boldsymbol{1}\bar{\boldsymbol{x}}^{\mathsf{T}} - \bar{\boldsymbol{x}}\boldsymbol{1}^{\mathsf{T}}\tilde{X} + \bar{\boldsymbol{x}}(\boldsymbol{1}^{\mathsf{T}}\boldsymbol{1})\bar{\boldsymbol{x}}^{\mathsf{T}} \\ &= \tilde{X}^{\mathsf{T}}\tilde{X} - n\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}} - n\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}} + n\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}} = \tilde{X}^{\mathsf{T}}\tilde{X} - n\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}} \end{aligned}$$

and similarly

$$X_{\mathrm{c}}^{\mathsf{T}}\boldsymbol{y}_{\mathrm{c}} = \tilde{X}^{\mathsf{T}}\boldsymbol{y} - \bar{y}\tilde{X}^{\mathsf{T}}\boldsymbol{1} - \bar{\boldsymbol{x}}\boldsymbol{1}^{\mathsf{T}}\boldsymbol{y} + \bar{y}\bar{\boldsymbol{x}}\boldsymbol{1}^{\mathsf{T}}\boldsymbol{1} = \tilde{X}^{\mathsf{T}}\boldsymbol{y} - n\bar{y}\bar{\boldsymbol{x}} .$$

This gives the two asserted equations, and conversely the original block system can be recovered from them, so they are equivalent. In other words, one may first centre the data, solve a system in $d$ unknowns, and finally recover the intercept as $\hat{b} = \bar{y}-\bar{\boldsymbol{x}}^{\mathsf{T}}\hat{\boldsymbol{u}}$. For $d = 1$ this is exactly the formula of Example 6.1.

Exercise 8.3Standard

Let $\lambda > 0$ and consider the ridge regression objective

$$L_{\lambda}(\boldsymbol{w}) = \lVert \boldsymbol{y}-X\boldsymbol{w}\rVert^2 + \lambda\lVert\boldsymbol{w}\rVert^2 .$$

(1) Show that $X^{\mathsf{T}}X+\lambda I_p$ is invertible regardless of the rank of $X$. (2) Show that $L_{\lambda}$ has exactly one minimiser, given by $\hat{\boldsymbol{w}}_{\lambda} = (X^{\mathsf{T}}X+\lambda I_p)^{-1}X^{\mathsf{T}}\boldsymbol{y}$.

Solution

(1). For $\boldsymbol{v}\ne\boldsymbol{0}$,

$$\boldsymbol{v}^{\mathsf{T}}(X^{\mathsf{T}}X+\lambda I_p)\boldsymbol{v} = \lVert X\boldsymbol{v}\rVert^2 + \lambda\lVert\boldsymbol{v}\rVert^2 \ge \lambda\lVert\boldsymbol{v}\rVert^2 > 0 .$$

If there were a $\boldsymbol{v}\ne\boldsymbol{0}$ with $(X^{\mathsf{T}}X+\lambda I_p)\boldsymbol{v} = \boldsymbol{0}$, the left-hand side would be $0$, a contradiction. Hence the kernel is $\{\boldsymbol{0}\}$ and the matrix is invertible.

(2). Consider the augmented design matrix and observation vector

$$X_{\lambda} = \begin{pmatrix} X \\ \sqrt{\lambda}\,I_p\end{pmatrix}\in\mathbb{R}^{(n+p)\times p}, \qquad \boldsymbol{y}_{\lambda} = \begin{pmatrix} \boldsymbol{y} \\ \boldsymbol{0}\end{pmatrix}\in\mathbb{R}^{n+p} .$$

Computing the norm block by block,

$$\lVert \boldsymbol{y}_{\lambda}-X_{\lambda}\boldsymbol{w}\rVert^2 = \lVert\boldsymbol{y}-X\boldsymbol{w}\rVert^2 + \lVert\sqrt{\lambda}\boldsymbol{w}\rVert^2 = L_{\lambda}(\boldsymbol{w}) ,$$

so ridge regression is an ordinary least squares problem for the augmented data. Since $X_{\lambda}^{\mathsf{T}}X_{\lambda} = X^{\mathsf{T}}X+\lambda I_p$ and $X_{\lambda}^{\mathsf{T}}\boldsymbol{y}_{\lambda} = X^{\mathsf{T}}\boldsymbol{y}$, by Theorem 3.3 its normal equations are $(X^{\mathsf{T}}X+\lambda I_p)\boldsymbol{w} = X^{\mathsf{T}}\boldsymbol{y}$. By (1) and Corollary 5.2 there is exactly one solution, given by the stated formula. As long as $\lambda > 0$, the answer is uniquely determined even for rank-deficient data such as that of Example 5.3; this is the benefit of regularisation.

Exercise 8.4Hard

Let $P\in\mathbb{R}^{n\times n}$ be symmetric and idempotent ($P^{\mathsf{T}} = P$, $P^2 = P$). (1) Show that the eigenvalues of $P$ are $0$ and $1$ only. (2) Show that $\operatorname{tr}P = \operatorname{rank}P$. (3) Verify that this is consistent with $\operatorname{tr}P = p$ for the hat matrix of Proposition 4.3.

Solution

(1). Suppose $P\boldsymbol{v} = \mu\boldsymbol{v}$ with $\boldsymbol{v}\ne\boldsymbol{0}$. Multiplying both sides by $P$ gives $P^2\boldsymbol{v} = \mu P\boldsymbol{v} = \mu^2\boldsymbol{v}$, while $P^2 = P$ makes the left-hand side $P\boldsymbol{v} = \mu\boldsymbol{v}$. Hence $(\mu^2-\mu)\boldsymbol{v} = \boldsymbol{0}$, and since $\boldsymbol{v}\ne\boldsymbol{0}$ we get $\mu^2 = \mu$, that is, $\mu\in\{0,1\}$.

(2). As $P$ is real symmetric, the spectral theorem gives an orthogonal matrix $U$ with $P = U\Lambda U^{\mathsf{T}}$ and $\Lambda = \operatorname{diag}(\mu_1,\dots,\mu_n)$. By (1) each $\mu_i$ is $0$ or $1$. By the cyclic property of the trace, $\operatorname{tr}P = \operatorname{tr}(\Lambda U^{\mathsf{T}}U) = \operatorname{tr}\Lambda = \#\{i : \mu_i = 1\}$. On the other hand $U$ is invertible, so $\operatorname{rank}P = \operatorname{rank}\Lambda$, which is also the number of ones. The two are therefore equal.

(3). The matrix $P = X(X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}$ of Proposition 4.3 is symmetric idempotent by part 1 of that proposition, and its image is $\operatorname{Im}X$ by part 2, so $\operatorname{rank}P = \dim\operatorname{Im}X = \operatorname{rank}X = p$. Combined with (2) this gives $\operatorname{tr}P = p$, agreeing with the result obtained directly from cyclicity in part 4 of the proposition. In Example 4.5 we had $n = 3$, $p = 2$, eigenvalues $1, 1, 0$ and trace $2$.

References

Textbooks (linear algebra and geometry). G. Strang, Introduction to Linear Algebra, 5th ed., Wellesley–Cambridge Press, 2016 — Chapter 4 (orthogonality, projections, least squares). The standard geometric account of projection matrices and the normal equations. Kenichi Kanatani, Kore nara Wakaru Saitekika Sugaku: Kiso Genri kara Keisan Shuho made, Kyoritsu Shuppan, 2005 (in Japanese) — the chapters on least squares and the singular value decomposition; a well-organised introduction for Japanese readers.

Numerical computation. G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013 — Chapter 5 (orthogonalisation and least squares problems). The standard source for QR-based solutions and the discussion of condition numbers. Å. Björck, Numerical Methods for Least Squares Problems, SIAM, 1996 — a monograph covering numerical methods for least squares problems comprehensively.

Connections to statistics and machine learning. T. Hastie, R. Tibshirani, and J. Friedman, The Elements of Statistical Learning, 2nd ed., Springer, 2009 — Chapter 3 (linear methods for regression). Ridge regression, multicollinearity and variable selection are explained as extensions of the framework of this article. A PDF is made available by the authors. C. M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006 — Chapter 3 (linear models for regression); a careful treatment of least squares reread as maximum likelihood.

History. S. M. Stigler, The History of Statistics: The Measurement of Uncertainty before 1900, Harvard University Press, 1986 — Part I. The Legendre–Gauss priority dispute and the process by which least squares came to be accepted are traced carefully from the sources.

Appendix: The pseudoinverse and the minimum-norm solution

Pinning down one answer even when the rank drops. As we saw in Example 5.3, when $\operatorname{rank}X < p$ there are infinitely many least squares solutions. A numerical library must return one of them. The standard choice is the one of smallest norm.

Proposition 8.5（The minimum-norm least squares solution）

Let $X\in\mathbb{R}^{n\times p}$ and $\boldsymbol{y}\in\mathbb{R}^n$, and let $S$ be the set of least squares solutions. Then $S \cap \operatorname{Im}(X^{\mathsf{T}})$ consists of a single point, and that element $\hat{\boldsymbol{w}}^{+}$ is the unique element of $S$ of smallest norm.

Proof（Proposition 8.5）

By Theorem 4.2 we have $S\ne\varnothing$, and by Theorem 3.3 $S = \boldsymbol{w}_0+\operatorname{Ker}X$ for any chosen $\boldsymbol{w}_0\in S$. Applying Lemma 4.1 to $X^{\mathsf{T}}$ gives $(\operatorname{Im}X^{\mathsf{T}})^{\perp} = \operatorname{Ker}X$, so we have the orthogonal decomposition

$$\mathbb{R}^p = \operatorname{Im}(X^{\mathsf{T}}) \oplus \operatorname{Ker}X .$$

Decomposing $\boldsymbol{w}_0 = \boldsymbol{u}+\boldsymbol{z}$ accordingly, with $\boldsymbol{u}\in\operatorname{Im}X^{\mathsf{T}}$ and $\boldsymbol{z}\in\operatorname{Ker}X$, we get $\boldsymbol{u} = \boldsymbol{w}_0-\boldsymbol{z} \in S$ and $\boldsymbol{u}\in\operatorname{Im}X^{\mathsf{T}}$, so $S\cap\operatorname{Im}(X^{\mathsf{T}})\ne\varnothing$. If two elements $\boldsymbol{u},\boldsymbol{u}'$ belong to this intersection, then $\boldsymbol{u}-\boldsymbol{u}'\in\operatorname{Ker}X$ (from the form of $S$) and $\boldsymbol{u}-\boldsymbol{u}'\in\operatorname{Im}X^{\mathsf{T}}$; the two summands of an orthogonal decomposition intersect in $\{\boldsymbol{0}\}$, so $\boldsymbol{u} = \boldsymbol{u}'$ and the intersection is a single point, which we denote $\hat{\boldsymbol{w}}^{+}$.

For minimality, every $\boldsymbol{w}\in S$ can be written $\boldsymbol{w} = \hat{\boldsymbol{w}}^{+}+\boldsymbol{v}$ with $\boldsymbol{v}\in\operatorname{Ker}X$, and $\hat{\boldsymbol{w}}^{+}\perp\boldsymbol{v}$, so the Pythagorean theorem gives

$$\lVert\boldsymbol{w}\rVert^2 = \lVert\hat{\boldsymbol{w}}^{+}\rVert^2 + \lVert\boldsymbol{v}\rVert^2 \ge \lVert\hat{\boldsymbol{w}}^{+}\rVert^2 ,$$

with equality only when $\boldsymbol{v} = \boldsymbol{0}$, that is, $\boldsymbol{w} = \hat{\boldsymbol{w}}^{+}$.

∎

The pseudoinverse. This $\hat{\boldsymbol{w}}^{+}$ can be written using the Moore–Penrose pseudoinverse built from the singular value decomposition $X = U\Sigma V^{\mathsf{T}}$.

Definition 8.6（The Moore–Penrose pseudoinverse）

Let $X = U\Sigma V^{\mathsf{T}}$ be a singular value decomposition of $X\in\mathbb{R}^{n\times p}$, where $U$ and $V$ are orthogonal and the diagonal entries of $\Sigma$ are the singular values $\sigma_1\ge\cdots\ge\sigma_r > 0$ with all other entries $0$. Let $\Sigma^{+}$ be the matrix obtained by replacing the non-zero diagonal entries of $\Sigma$ by their reciprocals and transposing. Then

$$X^{+} = V\Sigma^{+}U^{\mathsf{T}} \in \mathbb{R}^{p\times n}$$

is called the pseudoinverse of $X$. When $\operatorname{rank}X = p$ it agrees with $X^{+} = (X^{\mathsf{T}}X)^{-1}X^{\mathsf{T}}$.

One has $\hat{\boldsymbol{w}}^{+} = X^{+}\boldsymbol{y}$; this is checked by writing out, in terms of $U$, $\Sigma$ and $V$, that $X^{+}\boldsymbol{y}$ satisfies the normal equations and lies in $\operatorname{Im}X^{\mathsf{T}}$ (see Chapter 5 of Golub–Van Loan for details). This is the solution returned by NumPy’s np.linalg.lstsq and np.linalg.pinv. For the solution set $(0.5, 1.5-2t, t)$ of Example 5.3, the squared norm

$$0.25 + (1.5-2t)^2 + t^2 = 5t^2 - 6t + 2.5$$

is minimised at $t = 3/5 = 0.6$, giving $\hat{\boldsymbol{w}}^{+} = (0.5,\ 0.3,\ 0.6)^{\mathsf{T}}$. Indeed this solution is orthogonal to the generator $(0,2,-1)^{\mathsf{T}}$ of the kernel ($2\cdot 0.3 - 0.6 = 0$). Its squared norm is $0.7$, smaller than the $2.5$ of the solution $(0.5, 1.5, 0)$ at $t = 0$.