Real problems have many predictors at once: price depends on size and location and age. Multiple linear regression extends the line to a plane, then a hyperplane, and — crucially — introduces the matrix algebra and geometric picture that unify all of linear modeling. This chapter is the structural heart of the book: the normal equations and the projection view here reappear in regularization, GLMs, and the output layer of every neural network.

The model in matrix form

With \(p\) predictors, multiple linear regression models \[y_i=\beta_0+\beta_1 x_{i1}+\dots+\beta_p x_{ip}+\varepsilon_i .\] Stacking observations into a design matrix \(\mX\in\R^{n\times(p+1)}\) (a leading column of ones absorbs the intercept), a response vector \(\vy\in\R^n\), and coefficients \(\vbeta\in\R^{p+1}\), the whole dataset becomes \[\vy=\mX\vbeta+\veps .\]

Each row of \(\mX\) is one observation; each column is one feature across all observations. The product \(\mX\vbeta\) forms, for every row, the weighted sum of features — the prediction. Compressing \(n\) scalar equations into one matrix equation is not just shorthand: it lets us solve for all coefficients at once and exposes the geometry that makes least squares intuitive.

The normal equations

Minimizing \(\RSS(\vbeta)=\norm{\vy-\mX\vbeta}^2\) leads, by setting the gradient to zero, to the normal equations \[\mX^\top\mX\,\bhat=\mX^\top\vy\quad\Longrightarrow\quad \boxed{\;\bhat=(\mX^\top\mX)^{-1}\mX^\top\vy\;}\] (whenever \(\mX^\top\mX\) is invertible). This single formula is the workhorse of classical regression — it generalizes the slope/intercept formulas of Chapter [ch:simple] to any number of predictors.

The matrix \(\mX^{+}=(\mX^\top\mX)^{-1}\mX^\top\) is the Moore–Penrose pseudoinverse. If \(\mX\) were square and invertible we would just solve \(\vy=\mX\vbeta\) exactly; but \(\mX\) is tall (\(n>p\)), the system is overdetermined, and no exact solution exists. The pseudoinverse delivers the least-squares solution: the \(\bhat\) whose predictions come closest to \(\vy\).

The geometry: regression is projection

Here is the picture that makes least squares click.

The fitted values \(\vyhat=\mX\bhat\) are the orthogonal projection of \(\vy\) onto the column space of \(\mX\) — the subspace of all vectors reachable as \(\mX\vbeta\). Least squares finds the point in that subspace closest to \(\vy\); the residual \(\ve=\vy-\vyhat\) is the leftover, and it is perpendicular to every predictor: \(\mX^\top\ve=\mathbf 0\). That orthogonality is the normal equations.

Drop a perpendicular from the tip of \(\vy\) to the plane spanned by the predictors; the foot of that perpendicular is \(\vyhat\). “Least squares” literally means “shortest residual,” and the shortest connection from a point to a plane is the perpendicular one. The matrix \(\mH=\mX(\mX^\top\mX)^{-1}\mX^\top\) that performs this projection (\(\vyhat=\mH\vy\)) is called the hat matrix because it “puts the hat on \(\vy\).”

Interpreting coefficients: “holding others fixed”

In multiple regression, \(\beta_j\) is the expected change in \(y\) for a one-unit increase in \(x_j\) while all other predictors are held constant. This ceteris paribus interpretation is what makes multiple regression a tool for adjustment and control.

The phrase “holding others fixed” is treacherous when predictors are correlated (multicollinearity). If \(x_1\) and \(x_2\) move together in the data, the model cannot cleanly separate their effects: coefficients become unstable, have huge standard errors, and can even flip sign. A coefficient is not the raw association between \(x_j\) and \(y\) — it is the association after removing what the other predictors explain. This is also why a predictor can look important alone yet insignificant in the full model, and vice versa. We diagnose this with the variance inflation factor in Chapter [ch:diagnostics] and fix it with regularization in Chapter [ch:regularization].

When \(\mX^\top\mX\) is not invertible

If predictors are linearly dependent (e.g. redundant columns, or \(p>n\) as in genomics and modern ML), \(\mX^\top\mX\) is singular and the inverse does not exist — there are infinitely many least-squares solutions. Two remedies dominate: use the pseudoinverse (which picks the minimum-norm solution) or add a penalty that restores invertibility — exactly ridge regression, \(\bhat=(\mX^\top\mX+\lambda\mathbf I)^{-1}\mX^\top\vy\) (Chapter [ch:regularization]).

This chapter is the linear layer of a neural network. A Dense/Linear layer computes \(\mX\mathbf W\) — the identical matrix multiply, with the weight matrix playing the role of \(\vbeta\) and many outputs at once. The final layer of a regression network is multiple linear regression on learned features; the projection geometry explains why the network’s last hidden layer is a learned feature space in which the target is (approximately) a linear function. And the singular-\(\mX^\top\mX\) problem is ubiquitous in deep learning, where \(p\gg n\) — which is precisely why weight decay (ridge) is a default. Modern libraries avoid forming \((\mX^\top\mX)^{-1}\) explicitly, solving via QR or SVD for numerical stability, but the math is this.

import numpy as np
rng = np.random.default_rng(0)
n = 200
X = np.column_stack([np.ones(n), rng.normal(size=n), rng.normal(size=n)])
beta_true = np.array([1.0, 2.0, -0.5])
y = X @ beta_true + rng.normal(0, 0.5, size=n)
# Normal equations (use lstsq in practice for stability, not the explicit inverse)
beta_hat = np.linalg.solve(X.T @ X, X.T @ y)
print(beta_hat.round(3))                       # ~[1.0, 2.0, -0.5]
resid = y - X @ beta_hat
print(np.allclose(X.T @ resid, 0, atol=1e-9))  # residual orthogonal to columns -> True

Chapter summary

  • Multiple regression writes the whole dataset as \(\vy=\mX\vbeta+\veps\) with a design matrix \(\mX\).

  • The normal equations give \(\bhat=(\mX^\top\mX)^{-1}\mX^\top\vy\) — the multivariable least-squares solution.

  • Geometrically, least squares projects \(\vy\) onto the column space of \(\mX\); the residual is orthogonal to all predictors (\(\mX^\top\ve=\mathbf 0\)).

  • Each \(\beta_j\) is the effect of \(x_j\) holding others fixed; correlated predictors (multicollinearity) make this unstable.

  • Singular \(\mX^\top\mX\) (e.g. \(p>n\)) needs the pseudoinverse or a ridge penalty — the same linear-layer math used throughout deep learning.