1.1 Least squares and residuals

Regression begins with a modelling choice: which combinations of the measured variables should be allowed to explain the response? Once those combinations are fixed, least squares selects their coefficients by measuring the remaining error. The resulting calculation connects a data table, a quadratic objective, and an orthogonal projection. Understanding that connection makes it possible to check a fitted model rather than treating its coefficients as unexplained software output.

Motivation

Suppose an experiment records a predictor $x_i$ and a response $y_i$ on each of $n$ occasions. A line $a+bx$ offers a compact description of how the response changes with the predictor. Real observations usually do not lie on one line, so the equations $a+bx_i=y_i$ may be inconsistent. We can still ask for the line whose predictions come closest to the observations under a specified measure of error.

For this model, the residual at observation $i$ is the vertical difference $r_i=y_i-(a+bx_i)$. Its sign records whether the observed response lies above or below the prediction. Adding signed residuals alone would allow positive and negative errors to cancel. Squaring before adding produces a nonnegative objective and gives larger discrepancies greater influence. This choice also makes the objective a quadratic function of the unknown coefficients.

The data determine where the model is evaluated; the chosen features determine which prediction vectors are possible. A straight line, a polynomial, and a model with several measured predictors therefore share the same linear algebra once their features have been collected into a matrix.

Definitions

Definition

Design matrix, fitted values, and residuals

Let $A\in\mathbb R^{n\times p}$ contain one row per observation and one column per chosen feature. Let $y\in\mathbb R^n$ be the response vector and $\beta\in\mathbb R^p$ the coefficient vector. Define

$$\widehat y=A\beta,\qquad r=y-A\beta, \qquad J(\beta)=\|y-A\beta\|_2^2 =\sum_{i=1}^n r_i^2.$$

A least-squares coefficient vector minimizes $J$ over $\mathbb R^p$.

For a line, the feature functions are $1$ and $x$, so

$$A=\begin{bmatrix}1&x_1\\ \vdots&\vdots\\1&x_n\end{bmatrix}, \qquad \beta=\begin{bmatrix}a\\b\end{bmatrix}.$$

The column of ones supplies the intercept: its coefficient shifts every fitted response by the same amount. Omitting that column specifies a different model, one constrained to pass through the origin. The model should make that choice explicit before any calculation begins.

With $d$ measured predictors, row $i$ is $(1,x_{i,1},\ldots,x_{i,d})$ and $p=d+1$. The same construction allows fixed feature functions $\phi_1,\ldots,\phi_p$, with $A_{ij}=\phi_j(x_i)$. A polynomial model uses powers of the input as its features. Although its graph may curve, it remains linear in the unknown coefficients. Throughout the least-squares optimization, $A$ and $y$ are fixed; only $\beta$ varies.

The vector $\widehat y$ belongs to $\operatorname{col}(A)\subseteq\mathbb R^n$. This is a space of possible responses at all observed inputs, rather than a space of individual predictor values. Its ambient dimension is the number of observations. This distinction explains where the projection geometry lives.

Theorem/proposition

Proposition

Normal equations and orthogonal residuals

For any real matrix $A$, a vector $\widehat\beta$ minimizes $J$ exactly when

$$A^T(A\widehat\beta-y)=0, \qquad\text{equivalently}\qquad A^TA\widehat\beta=A^Ty.$$

Its residual $\widehat r=y-A\widehat\beta$ is perpendicular to every column of $A$. Consequently, $A\widehat\beta$ is the orthogonal projection of $y$ onto $\operatorname{col}(A)$.

Proposition

When coefficients are unique

The matrix $A^TA$ is invertible exactly when $A$ has linearly independent columns. In that case the unique least-squares coefficient vector is

$$\widehat\beta=(A^TA)^{-1}A^Ty.$$

For an arbitrary $A$, the fitted vector is still unique. If $\widehat\beta$ is one minimizer, all minimizers are $\widehat\beta+z$ with $z\in\ker(A)$.

These statements separate two questions that arise in modelling. The fitted responses may be determined even when the data cannot determine how their contributions should be divided among redundant features. Having more observations than coefficients does not by itself remove that redundancy: independence of the columns is the decisive condition.

Proof sketch or proof idea

For a candidate $\widehat\beta$, write $\widehat r=y-A\widehat\beta$. A change $h$ in the coefficients gives the exact identity

$$J(\widehat\beta+h) =\|\widehat r-Ah\|_2^2 =\|\widehat r\|_2^2-2h^TA^T\widehat r+\|Ah\|_2^2.$$

If $A^T\widehat r=0$, the cross term disappears. Every coefficient change therefore adds the nonnegative amount $\|Ah\|_2^2$ to the objective, proving global minimality. Conversely, at a minimizer the function of one real variable $t\mapsto J(\widehat\beta+th)$ has derivative zero at $t=0$ for every $h$. Its derivative there is $-2h^TA^T\widehat r$, forcing $A^T\widehat r=0$. This proves both directions, including cases with dependent columns.

Proof

Rank, existence, and the unique fitted vector

For every $z\in\mathbb R^p$,

$$z^TA^TAz=\|Az\|_2^2.$$

Thus $A^TA$ has zero kernel precisely when $A$ does. For independent columns this quadratic form is positive for every nonzero $z$, and the inverse formula follows from the normal equations.

For existence, project $y$ onto the finite-dimensional subspace $\operatorname{col}(A)$. The projection has the form $A\widehat\beta$, and its orthogonal residual satisfies the normal equations. If two coefficient vectors satisfy those equations, their difference $z$ satisfies $A^TAz=0$. Multiplying by $z^T$ gives $\|Az\|_2^2=0$, so both vectors produce exactly the same fitted responses. Adding a vector in $\ker(A)$ leaves those responses unchanged, which describes the entire family of minimizers.

An intercept has an immediate diagnostic consequence. Its column is the vector $\mathbf1=(1,\ldots,1)^T$, so the normal equations include $\mathbf1^T\widehat r=0$. Hence the residuals sum to zero and the average fitted response equals the average observed response. Other feature columns impose their own orthogonality equations. These equations test optimality jointly; checking the residual sum alone tests only the intercept equation.

Worked examples

Worked example

Example 1: Fit a line and verify its residuals

Use the three observations $(0,1)$, $(1,2)$, and $(2,2)$. Then

$$A=\begin{bmatrix}1&0\\1&1\\1&2\end{bmatrix}, \qquad y=\begin{bmatrix}1\\2\\2\end{bmatrix}, \qquad A^TA=\begin{bmatrix}3&3\\3&5\end{bmatrix}, \qquad A^Ty=\begin{bmatrix}5\\6\end{bmatrix}.$$

The normal equations are $3a+3b=5$ and $3a+5b=6$. Subtracting yields $2b=1$, followed by $a=7/6$. Therefore

$$\widehat y=\begin{bmatrix}7/6\\5/3\\13/6\end{bmatrix}, \qquad \widehat r=\begin{bmatrix}-1/6\\1/3\\-1/6\end{bmatrix}, \qquad J(\widehat\beta)=\frac16.$$

The residual sum is zero. Its dot product with the predictor column is also zero: $0(-1/6)+1(1/3)+2(-1/6)=0$. These two checks verify the full normal equations. The determinant of $A^TA$ is $6$, so the coefficients are unique. The positive residual sum of squares quantifies the fact that the three observations do not lie on one straight line.

The line also passes through the sample-mean point $(1,5/3)$, as follows from the intercept equation. This is a consequence of the fit, rather than an extra constraint added to the optimization. It supplies a useful arithmetic check when solving a small regression problem by hand.

Worked example

Example 2: Unique fitted values with ambiguous coefficients

Suppose all three inputs equal $2$, while the responses are $1$, $2$, and $3$. The model still has the form $a+bx$, but

$$A=\begin{bmatrix}1&2\\1&2\\1&2\end{bmatrix}, \qquad A^TA=\begin{bmatrix}3&6\\6&12\end{bmatrix}.$$

Every prediction equals the single value $t=a+2b$. The objective becomes

$$(1-t)^2+(2-t)^2+(3-t)^2=2+3(t-2)^2.$$

It is minimized at $t=2$, so all coefficient vectors $(a,b)=(2-2s,s)$, with $s\in\mathbb R$, are minimizers. Each gives $\widehat y=(2,2,2)^T$, residual $(-1,0,1)^T$, and objective value $2$. The normal equations consistently reduce to $a+2b=2$.

The nullspace direction is $(-2,1)^T$: changing the coefficients along it has no effect at the observed input. At the new input $x=3$, however, these models predict $2+s$. Thus agreement of fitted responses does not identify a unique prediction rule away from the observed inputs. The data contain no variation in $x$ from which to determine the slope.

Worked example

Example 3: Weights on residuals and weights in the objective

Fit a constant $c$ to responses $0$, $2$, and $8$, with objective weights $q=(1,1,1/4)^T$. The objective is

$$J_q(c)=c^2+(2-c)^2+\frac14(8-c)^2.$$

Its derivative vanishes when

$$\left(1+1+\frac14\right)c =0+2+\frac14\,8=4, \qquad \widehat c=\frac{16}{9}.$$

The residuals are $(-16/9,2/9,56/9)^T$. Their weighted sum is zero: $-16/9+2/9+(1/4)(56/9)=0$. The weighted objective value is $116/9$. For comparison, equal objective weights give the ordinary mean $10/3$. Reducing the influence of the largest response moves the fitted constant toward the other two observations.

To express this problem as ordinary least squares, multiply the residuals by $W=\operatorname{diag}(1,1,1/2)$. Squaring their lengths then produces the objective weights $1$, $1$, and $1/4$. The distinction matters: multiplying the last residual by $1/4$ would instead give it objective weight $1/16$ and solve a different problem.

Weighted models and computation

In general, positive objective weights $q_i$ define

$$J_q(\beta)=\sum_i q_i(y_i-(A\beta)_i)^2 =\|W(y-A\beta)\|_2^2, \qquad W=\operatorname{diag}(\sqrt{q_1},\ldots,\sqrt{q_n}).$$

Apply the normal equations to the transformed matrix $WA$ and response $Wy$:

$$A^TW^2A\widehat\beta=A^TW^2y.$$

Equivalently, if $Q=\operatorname{diag}(q_i)$, the equation is $A^TQA\widehat\beta=A^TQy$. This convention keeps the objective weights and the factors multiplying individual residuals distinct. Positive weights make $W$ invertible, so $WA$ has the same column rank as $A$. The uniqueness condition is consequently unchanged. Orthogonality now means $A^TQ\widehat r=0$, with respect to the weighted inner product.

For numerical computation, solve the least-squares problem directly. The inverse formula expresses the mathematical solution when full column rank holds; a direct least-squares routine can use orthogonal factorizations without explicitly forming that inverse. For the first example:

import numpy as np

x = np.array([0.0, 1.0, 2.0])
y = np.array([1.0, 2.0, 2.0])
A = np.column_stack((np.ones_like(x), x))
beta, _, rank, singular_values = np.linalg.lstsq(A, y, rcond=None)
residual = y - A @ beta
rss = residual @ residual
# beta is approximately [1.16666667, 0.5]; rss is approximately 1/6

The returned rank describes how many independent feature directions the routine detects. Computing the residual vector explicitly makes its sign convention and the objective visible. For positive objective weights stored in q, the corresponding call uses np.sqrt(q)[:, None] * A and np.sqrt(q) * y. The singular values also connect this computation to dimensionality reduction.

QR and triangular least-squares solving

When $A$ has full column rank and $n\ge p$, orthonormalizing its columns gives a thin QR factorization

$$A=QR,\qquad Q\in\mathbb R^{n\times p},\quad Q^TQ=I_p,\qquad R\in\mathbb R^{p\times p}.$$

Here $R$ is invertible and upper triangular. Split the response into its projection $QQ^Ty$ and its orthogonal remainder. Because $Q^T(I-QQ^T)=0$ and $Q$ preserves lengths of coefficient vectors,

$$\|A\beta-y\|_2^2 =\|R\beta-Q^Ty\|_2^2+\|(I-QQ^T)y\|_2^2.$$

The second term does not depend on $\beta$. The first can be made zero, so the unique minimizer solves $R\widehat\beta=Q^Ty$. Back substitution starts with the last row of this triangular system, then uses the known coefficient in the preceding row. This route gives the same normal-equation solution while avoiding explicit construction of $A^TA$ or its inverse. The full column-rank condition makes $R$ invertible; a rank-deficient design instead needs the least-squares or pseudoinverse treatment above.

Worked example

A triangular regression calculation

Use rows $(1,0),(1,1),(1,2)$ and response $y=(1,1,4)^T$. Choose columns of $Q$ as $(1,1,1)^T/\sqrt3$ and $(-1,0,1)^T/\sqrt2$. Then

$$R=\begin{pmatrix}\sqrt3&\sqrt3\\0&\sqrt2\end{pmatrix},\qquad Q^Ty=\begin{pmatrix}2\sqrt3\\3/\sqrt2\end{pmatrix},\qquad \widehat\beta=\begin{pmatrix}1/2\\3/2\end{pmatrix}.$$

The last triangular equation gives slope $3/2$; the first gives intercept $1/2$. The fitted response is $(1/2,2,7/2)^T$, and the residual is $(1/2,-1,1/2)^T$. Its squared norm is $3/2$, exactly the orthogonal remainder term in the decomposition.

Feature choice and the meaning of the objective

A quadratic response model uses the three features $1$, $x$, and $x^2$. For the inputs in Example 1 its design matrix is

$$A_2=\begin{bmatrix}1&0&0\\1&1&1\\1&2&4\end{bmatrix}.$$

The coefficients of the polynomial $1+(3/2)x-(1/2)x^2$ reproduce all three responses, giving zero training residual. Each former straight-line model remains available by setting the quadratic coefficient to zero, so this smaller minimum also follows from inclusion of the feature spaces. The extra feature changes the set over which we minimize; it does not improve the solution of the original two-feature optimization problem.

At a new input $x=3$, this quadratic predicts $1$, whereas the fitted line predicts $8/3$. Both calculations are determined by their chosen models, but the original three observations provide no response at that new input with which to compare them. Selecting a feature space and estimating coefficients within it are therefore different tasks. Least squares solves the latter once the former is fixed.

Transforming a response also changes what an error means. For observations satisfying $x_i>0$ and $y_i>0$, fit the power model $y=ax^b$ with $a>0$. Logarithms give

$$\log y=\log a+b\log x.$$

Least squares with response $\log y_i$ and features $1$, $\log x_i$ minimizes squared differences on the logarithmic response scale. Exponentiating the fitted intercept recovers the coefficient $a$, but does not turn that objective into squared differences on the original response scale. On the logarithmic scale, equal multiplicative ratios give equal residual differences; on the original scale, equal additive differences do. This is a mathematical change to the fitting criterion, not simply a convenient rewrite of the same loss. The positivity assumptions specify the domain in which the transformation is defined.

Other linearized models

The same design-matrix method handles several other models. The response column determines the scale on which squared errors are minimized:

For the exponential model, solve $\log y=c+dx$, then recover $a=e^c$ and $b=e^d$. Its inputs can be any real numbers; only its observed responses must be positive for this logarithmic fit. Logarithmic regression instead transforms only the inputs, so its responses may be negative or zero. The input-reciprocal model also preserves the original response SSE. Transforming the response by a logarithm or reciprocal selects a different loss.

For each two-column design above, unique coefficients require that its transformed input values are not all equal. Linearization makes coefficients linear in the transformed equation; it does not remove a rank deficiency.

Worked example

Example 4: Exact exponential and input-reciprocal fits

For $(0,2)$, $(1,6)$, and $(2,18)$, the log responses lie on

$$\log y=\log2+x\log3.$$

Thus the transformed least-squares fit has $c=\log2$, $d=\log3$, and the recovered model is $y=2\cdot3^x$. The input zero is allowed here because the input itself is not logged.

For $(1,5)$, $(2,3)$, and $(4,2)$, use the transformed inputs $1$, $1/2$, and $1/4$. They give the exact line $y=1+4/x$. This fit minimizes original response SSE, which is zero. It is defined for every $x\ne0$, including negative inputs; there is no response-positivity restriction.

Worked example

Example 5: Taking response reciprocals changes the objective

Take $x=(0,1,2)^T$ and $y=(1,1/2,1/4)^T$. Fitting the reciprocal responses $(1,2,4)^T$ with features $1,x$ gives

$$\widehat a=\frac56,\qquad \widehat b=\frac32, \qquad \sum_i\left(\frac1{y_i}-(\widehat a+\widehat b x_i)\right)^2=\frac16.$$

Restoring the response gives predictions $(6/5,3/7,6/23)^T$. Their original response SSE is $1/25+1/196+1/8464$, rather than $1/6$. In general,

$$\left(\frac1{y_i}-\frac1{\widehat y_i}\right)^2 =\frac{(y_i-\widehat y_i)^2}{y_i^2\widehat y_i^2}, \qquad y_i\ne0,\quad\widehat y_i\ne0.$$

The denominator depends on the prediction, so this is not the original SSE or a predetermined observation-weighted SSE. The fitted denominator is nonzero at all three observations, but vanishes at $x=-5/9$, where the model has no defined response. Evaluate a model's domain as well as its chosen loss before using a transformed fit for prediction.

Finally, interpret residual checks relative to the actual feature matrix. If Example 1 is fitted through the origin using only the feature $x$, its one coefficient is $6/5$. The residual becomes $(1,4/5,-2/5)^T$. Its dot product with the predictor column is zero, as optimality requires, while its entries sum to $7/5$. There is no intercept column imposing a zero-sum condition. Changing the model changes the diagnostic equations along with the fitted values. Recording the matrix next to the coefficients keeps those equations unambiguous.

Common mistakes

Common mistake

Confusing three different claims about a fit

The normal equations characterize the smallest training residual sum of squares within the chosen feature space. A zero residual additionally means that the observations lie in that space. Unique coefficients additionally require independent feature columns. These are three separate mathematical claims, each with its own test.

Adding features enlarges the possible fitted-response space and can only decrease its minimum training error, because the old model remains available by setting the new coefficients to zero. Prediction on unseen observations is a further modelling question. A smaller training error alone does not compare those prediction errors.

A second frequent error is to assume that every residual is perpendicular to the fitted line drawn in the predictor-response plane. Ordinary least squares uses vertical response differences. The perpendicularity theorem concerns the whole residual vector in $\mathbb R^n$ and the columns of the design matrix. For weighted problems, the corresponding theorem uses the weighted inner product. Keeping the ambient space explicit resolves this apparent geometric conflict.

Summary

Choose features first, assemble one design-matrix row per observation, and minimize the squared residual norm. Check a candidate with the normal equations. Column rank determines coefficient uniqueness; orthogonal projection determines the unique fitted vector. With weighted observations, take square roots of objective weights when scaling rows for an ordinary least-squares solver.

Exercises

Checkpoint

Exercise 1: Fit a constant to responses $1$, $4$, and $7$. Find the minimizer, residual vector, and residual sum of squares.

The design matrix has a single column of ones. Complete a square or use its one normal equation.

Checkpoint

Exercise 2: Let $A=\begin{bmatrix}1&0\\1&1\\1&2\end{bmatrix}$ and propose residual $r=(1,-1,0)^T$. Can it be the residual of an ordinary least-squares fit with this design matrix?

Check its inner product with every column, including the predictor column.

Checkpoint

Exercise 3: A constant model has responses $0$ and $6$, with objective weights $1$ and $2$. Find the fitted constant and the diagonal matrix that scales residuals for ordinary least squares.

Write the weighted objective before differentiating. Distinguish an objective weight from its square root.

Checkpoint

Exercise 4: Fit $y=a+b\log x$ to $(1,-1)$, $(e,1)$, and $(e^2,3)$. Find the exact model. Does the negative response prevent this fit?

Identify which variable is transformed before applying any sign condition.

Solutions

Solution · Exercise 1 solution

The normal equation is $3c=1+4+7=12$, so $\widehat c=4$. The residual vector is $(-3,0,3)^T$ and the residual sum of squares is $9+0+9=18$. Its entries sum to zero, as required by the intercept column. Completing a square gives $J(c)=18+3(c-4)^2$, which also proves uniqueness.

Solution · Exercise 2 solution

The residual sum is zero, but

$$A^Tr=\begin{bmatrix}1-1+0\\0(1)+1(-1)+2(0)\end{bmatrix} =\begin{bmatrix}0\\-1\end{bmatrix}.$$

The second normal equation fails, so this cannot be a least-squares residual for the specified matrix. Orthogonality to the intercept column alone is insufficient.

Solution · Exercise 3 solution

The objective is $c^2+2(6-c)^2$. Its derivative is $2c+4(c-6)$, giving $\widehat c=4$. The residual vector $(-4,2)^T$ has weighted sum $1(-4)+2(2)=0$. The scaling matrix is $W=\operatorname{diag}(1,\sqrt2)$, because $W^TW=\operatorname{diag}(1,2)$.

Solution · Exercise 4 solution

The transformed inputs are $0,1,2$, with responses $-1,1,3$. Thus $a=-1$, $b=2$, and the exact model is $y=-1+2\log x$. Every input is positive. The response is not logged, so its negative value is allowed; the minimized loss is original response SSE.

Related notes

SVD and principal components uses singular vectors to describe important directions in a data matrix. Gradient descent on quadratics studies iterative optimization of objectives such as the one derived here. The course reading list collects the books for this course.