Given a set of training data \(\mathcal{D}=\{(x_i,y_i)\}_{i=1}^{n}\) with \(x_i,y_i\in\mathbb{R}\), regression asks for a function \(f\) that best explains how \(y\) depends on \(x\). This is one of the oldest and most useful problems in applied mathematics: predicting a continuous outcome from an observed input, fitting a trend line through noisy measurements, or quantifying how strongly two quantities are related.
In the simplest case, \(f\) is assumed to be a straight line, \(f(x)=mx+b\). What looks like a two-line-of-algebra problem at first turns out to be a small, self-contained piece of applied mathematics: it needs a model, a notion of error, an optimization criterion, and a closed-form solution, all before a single line of code is written. What follows builds all of that from the ground up, then generalizes it to many features, curved relationships, and ill-conditioned data.
What Is Regression?
A regression model is a function \(y=f(x;\,a,b,\dots)\) with parameters \(a,b,\dots\) that is fit to \(\mathcal{D}\) so that its predictions match the observed \(y_i\) as closely as possible. The parameters are the unknowns; \(\mathcal{D}\) is given. Once fitted, \(f\) can be evaluated at new inputs to predict an outcome, or simply inspected to quantify how \(x\) and \(y\) relate.
The simplest and most common regression model is linear regression, where \(f\) is an affine function of its parameters. Everything below is developed for the one-feature case first, since it is fully solvable by hand and builds the intuition needed for the general, multi-feature case afterward.
The Simple Linear Model
A straight line is determined by two points. With three or more data points, a line will in general miss every single one of them, since real measurements are noisy and rarely fall exactly on a common line. To express this honestly, an error term \(\epsilon\) is attached to the model for every data point, describing the gap between the line and the actual measurement:
\[ f(x) = mx+b+\epsilon. \]
For a given \(x_i\), the model should predict \(y_i\). The residual, i.e. the error actually made at that point, is
\[ \epsilon_i = f(x_i) - y_i = mx_i+b-y_i. \]
A useful mental picture is a rigid rod resting against a set of springs, one anchored at each data point. The rod settles into the position that balances all the springs at once—this is exactly the line produced by the fitting procedure derived next, and it is why the method is comparatively insensitive to any single point but very sensitive to a stretched-out spring far away from the rest.
Deriving the Least-Squares Estimators
To turn "the errors should be small" into a solvable problem, the individual residuals are combined into a single number to minimize. Squaring each residual before summing serves two purposes: it removes the sign, so positive and negative errors cannot cancel out, and it penalizes large deviations more than small ones, which favors a solution with many small errors over one with a single very large error. The result is the loss function
\[ \mathcal{L}(m,b) = \sum_{i=1}^n\epsilon_i^2 = \sum_{i=1}^n\bigl(f(x_i)-y_i\bigr)^2 = \sum_{i=1}^n(mx_i+b-y_i)^2. \]
As a sum of squares in \(m\) and \(b\), \(\mathcal{L}\) is a convex quadratic function, i.e. it is shaped like a bowl that opens upward. Consequently, it has exactly one stationary point, and that stationary point is its global minimum. Finding the best-fitting line therefore reduces to setting both partial derivatives of \(\mathcal{L}\) to zero:
\[ \frac{\partial}{\partial m}\mathcal{L}(m,b) = -2\sum_{i=1}^n x_i(y_i-mx_i-b) \overset{!}{=} 0, \qquad \frac{\partial}{\partial b}\mathcal{L}(m,b) = -2\sum_{i=1}^n (y_i-mx_i-b) \overset{!}{=} 0. \]
It is easiest to start with the second equation and solve for \(b\), since it does not involve an extra factor of \(x_i\):
\[ \sum_{i=1}^n y_i - m\sum_{i=1}^n x_i - nb = 0 \quad\Longleftrightarrow\quad b = \frac{1}{n}\sum_{i=1}^n y_i - m\cdot\frac{1}{n}\sum_{i=1}^n x_i = \overline{y}-m\overline{x}, \]
where \(\overline{x}\) and \(\overline{y}\) denote the arithmetic means of the \(x_i\) and \(y_i\). This already makes intuitive sense: the fitted line always passes through the point \((\overline{x},\overline{y})\), the "center of mass" of the data, regardless of the slope.
Substituting this expression for \(b\) into the first equation eliminates \(b\) and leaves a single equation in \(m\):
\[ \sum_{i=1}^n x_i\bigl(y_i - mx_i - (\overline{y}-m\overline{x})\bigr) = 0 \;\Longleftrightarrow\; \sum_{i=1}^n x_i(y_i-\overline{y}) = m\sum_{i=1}^n x_i(x_i-\overline{x}). \]
Because \(\sum_{i=1}^n(\overline{x}y_i-\overline{x}\,\overline{y})=0\) and \(\sum_{i=1}^n(x_i\overline{x}-\overline{x}^2)=0\) (both sums telescope to zero, since they only involve constants times a deviation from the mean, which itself sums to zero), the identity above can be shifted to a fully centered form without changing its value:
\[ \sum_{i=1}^n (x_i-\overline{x})(y_i-\overline{y}) = m\sum_{i=1}^n (x_i-\overline{x})^2. \]
Solving for \(m\) gives the slope of the least-squares line:
\[ \boxed{\,m = \frac{\displaystyle\sum_{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}{\displaystyle\sum_{i=1}^n(x_i-\overline{x})^2} = \frac{\operatorname{Cov}(x,y)}{\operatorname{Var}(x)}\,}, \qquad b = \overline{y}-m\overline{x}. \]
Here \(\operatorname{Cov}(x,y)\) is the covariance, measuring how \(x\) and \(y\) vary together, and \(\operatorname{Var}(x)\) is the variance of \(x\) alone. Both are ordinarily divided by \(n\) (or \(n-1\) for the unbiased sample estimate), but since that factor appears in both the numerator and the denominator of \(m\), it cancels and can be dropped entirely; the two raw sums above are all that is needed.
Worked Example
Consider five measurements:
| \(i\) | 1 | 2 | 3 | 4 | 5 |
|---|---|---|---|---|---|
| \(x_i\) | 1 | 2 | 3 | 4 | 5 |
| \(y_i\) | 2 | 4 | 5 | 4 | 5.5 |
The means are \(\overline{x}=3\) and \(\overline{y}=4.1\). Subtracting the means and combining terms:
| \(x_i-\overline{x}\) | \(-2\) | \(-1\) | \(0\) | \(1\) | \(2\) |
|---|---|---|---|---|---|
| \(y_i-\overline{y}\) | \(-2.1\) | \(-0.1\) | \(0.9\) | \(-0.1\) | \(1.4\) |
| Product | \(4.2\) | \(0.1\) | \(0\) | \(-0.1\) | \(2.8\) |
| \((x_i-\overline{x})^2\) | \(4\) | \(1\) | \(0\) | \(1\) | \(4\) |
Summing the last two rows gives \(\sum(x_i-\overline{x})(y_i-\overline{y}) = 7\) and \(\sum(x_i-\overline{x})^2 = 10\), so
\[ m = \frac{7}{10} = 0.7, \qquad b = 4.1 - 0.7\cdot 3 = 2.0. \]
The fitted line is \(y=0.7x+2.0\), drawn above together with the residual "springs" connecting each point to the line it was fitted against.
Goodness of Fit: Correlation and \(R^2\)
The slope alone does not say how well the line actually fits the data. The correlation coefficient
\[ r = \frac{\operatorname{Cov}(x,y)}{\sqrt{\operatorname{Var}(x)\operatorname{Var}(y)}} \]
rescales the covariance to always lie in \([-1,1]\), by the Cauchy–Schwarz inequality applied to the centered vectors \((x_i-\overline{x})\) and \((y_i-\overline{y})\). A value near \(\pm1\) indicates that the points lie close to a line, and the sign of \(r\) matches the sign of the slope \(m\). A value near \(0\) means there is little to no linear relationship, though a strong nonlinear one may still exist.
Squaring \(r\) gives the coefficient of determination \(R^2=r^2\), interpreted as the fraction of the variance in \(y\) that is explained by the linear relationship with \(x\); the remaining \(1-R^2\) is unexplained, left in the residuals. For the worked example, \(r\approx 0.82\) and \(R^2\approx 0.68\): about 68% of the spread in \(y\) is accounted for by the fitted line.
This interpretation is not just a rule of thumb, it follows directly from the least-squares construction itself. Write each observation as \(y=\hat y+\epsilon\), the sum of its prediction \(\hat y=mx+b\) and its residual. The covariance between predictions and residuals vanishes:
\[ \operatorname{Cov}(\hat y,\epsilon) = \operatorname{Cov}(mx+b,\,y-mx-b) = m\operatorname{Cov}(x,y)-m^2\operatorname{Var}(x) = m\bigl(\operatorname{Cov}(x,y)-m\operatorname{Var}(x)\bigr) = 0, \]
because \(m=\operatorname{Cov}(x,y)/\operatorname{Var}(x)\) makes the term in parentheses zero, the very same condition used to derive \(m\) in the first place. Since predictions and residuals are uncorrelated, the variance of \(y=\hat y+\epsilon\) splits additively,
\[ \operatorname{Var}(y) = \operatorname{Var}(\hat y) + \operatorname{Var}(\epsilon), \]
and \(R^2=\operatorname{Var}(\hat y)/\operatorname{Var}(y)\) is exactly the share of that total variance contributed by the fitted line, with the rest, \(\operatorname{Var}(\epsilon)/\operatorname{Var}(y)=1-R^2\), left unexplained in the residuals.
Regression Toward the Mean
Standardizing both variables to \(z\)-scores, mean \(0\) and variance \(1\), makes the slope formula collapse to something remarkably simple. With \(\operatorname{Var}(x)=1\), the least-squares slope \(m=\operatorname{Cov}(x,y)/\operatorname{Var}(x)\) is just the correlation coefficient \(r\) itself, so the fitted line for standardized data is
\[ \hat z_y = r\,z_x. \]
Because \(|r|\le 1\), a predictor that is \(k\) standard deviations away from its own mean predicts an outcome only \(rk\) standard deviations away from its mean, strictly closer to the mean whenever \(|r|<1\). This is regression toward the mean, and it is where the whole technique gets its name: an unusually extreme input predicts a less extreme output.
Francis Galton observed exactly this in 1886, studying the heights of parents and their adult children: tall parents tend to have children who are taller than average, but shorter than their parents; short parents tend to have children shorter than average, but taller than themselves. Nothing shrinks from one generation to the next, a family's height is not "regressing" anywhere. What regresses is the prediction, precisely because \(|r|<1\) whenever height is not perfectly heritable, and part of a child's height comes from sources uncorrelated with the parent's.
Partial Correlation
A strong correlation between \(x\) and \(y\) does not mean one causes the other; both may instead be driven by some third variable \(z\). Partial correlation makes this precise: it measures what remains of the association between \(x\) and \(y\) once the part explained by \(z\) is regressed out of both.
Standardize all three variables, regress \(x\) on \(z\) and \(y\) on \(z\) separately using the slope formula from above, and keep only the residuals \(e_x=x-r_{xz}z\) and \(e_y=y-r_{yz}z\), the parts of \(x\) and \(y\) that \(z\) does not explain. Since standardized variables have \(\operatorname{Var}(z)=1\) and \(\operatorname{Cov}(x,z)=r_{xz}\),
\[ \operatorname{Cov}(e_x,e_y) = \operatorname{Cov}(x,y) - r_{xz}\operatorname{Cov}(z,y) - r_{yz}\operatorname{Cov}(x,z) + r_{xz}r_{yz}\operatorname{Var}(z) = r_{xy}-r_{xz}r_{yz}, \]
while \(\operatorname{Var}(e_x)=\operatorname{Var}(x)-2r_{xz}\operatorname{Cov}(x,z)+r_{xz}^2 \operatorname{Var}(z) = 1-r_{xz}^2\), and likewise \(\operatorname{Var}(e_y)=1-r_{yz}^2\). The partial correlation is the ordinary correlation between these two residuals,
\[ \boxed{r_{xy\cdot z} = \frac{\operatorname{Cov}(e_x,e_y)}{\sqrt{\operatorname{Var}(e_x)\operatorname{Var}(e_y)}} = \frac{r_{xy}-r_{xz}r_{yz}}{\sqrt{(1-r_{xz}^2)(1-r_{yz}^2)}}}. \]
A classic illustration: the number of fire trucks sent to a house fire correlates strongly with the property damage the fire causes, say \(r_{xy}=0.95\). This is not because trucks cause damage; both are driven by the size of the fire \(z\), which correlates with the number of trucks dispatched (\(r_{xz}=0.97\), bigger fires get a bigger response) and with the damage (\(r_{yz}=0.99\), bigger fires do more damage). Controlling for fire size:
\[ r_{xy\cdot z} = \frac{0.95-0.97\cdot 0.99}{\sqrt{(1-0.97^2)(1-0.99^2)}} \approx -0.30. \]
Once the size of the fire is held fixed, sending more trucks is if anything mildly associated with less damage, the opposite sign of the raw correlation. The original \(r_{xy}=0.95\) said nothing about trucks causing damage; it only reflected the shared dependence of both variables on a fire's size.
Anscombe's Quartet
Every quantity introduced so far, the slope, the intercept, the correlation, and \(R^2\), is a summary statistic: a single number computed from the full scatter of points. Anscombe (1973) constructed four small data sets that make a memorable point about summary statistics: all four share the same mean, variance, correlation, regression line, and \(R^2\), to two decimal places, yet look completely different once plotted.
| \(i\) | \(x\) | \(y_1\) | \(y_2\) | \(y_3\) | \(x_4\) | \(y_4\) |
|---|---|---|---|---|---|---|
| 1 | 10 | 8.04 | 9.14 | 7.46 | 8 | 6.58 |
| 2 | 8 | 6.95 | 8.14 | 6.77 | 8 | 5.76 |
| 3 | 13 | 7.58 | 8.74 | 12.74 | 8 | 7.71 |
| 4 | 9 | 8.81 | 8.77 | 7.11 | 8 | 8.84 |
| 5 | 11 | 8.33 | 9.26 | 7.81 | 8 | 8.47 |
| 6 | 14 | 9.96 | 8.10 | 8.84 | 8 | 7.04 |
| 7 | 6 | 7.24 | 6.13 | 6.08 | 8 | 5.25 |
| 8 | 4 | 4.26 | 3.10 | 5.39 | 19 | 12.50 |
| 9 | 12 | 10.84 | 9.13 | 8.15 | 8 | 5.56 |
| 10 | 7 | 4.82 | 7.26 | 6.42 | 8 | 7.91 |
| 11 | 5 | 5.68 | 4.74 | 5.73 | 8 | 6.89 |
All four data sets share \(\overline{x}\approx 9.00\), \(\overline{y}\approx 7.50\), \(\operatorname{Var}(x)\approx 11.0\), \(\operatorname{Var}(y)\approx 4.13\), \(r\approx 0.816\), the fitted line \(\hat y=3.00+0.500x\), and \(R^2\approx 0.67\). Plotted, they tell four different stories:
- Set I is what the least-squares assumptions were built for: a genuinely linear relationship with roughly evenly scattered noise.
- Set II is a smooth curve, not a line. The straight fit systematically overshoots in the middle and undershoots at both ends, a pattern that residual plots and the feature maps introduced next are meant to catch.
- Set III is nearly a perfect line for ten of its eleven points; a single outlier pulls the fitted line away from what would otherwise be an almost exact fit.
- Set IV has no variability in \(x\) at all except for one point. That single observation, far out on the horizontal axis, single-handedly determines the entire slope, an extreme case of a high-leverage point.
The lesson is not that summary statistics are useless, but that they are not a substitute for looking at the data: always plot the scatterplot and the residuals before trusting a fitted line.
Multiple Linear Regression
Real problems rarely depend on a single input. With \(D\) features, each observation becomes a vector \(\mathbf{x}=(x_1,\dots,x_D)\), the slope becomes a weight vector \(\mathbf{w}=(w_1,\dots,w_D)\), and the model generalizes to
\[ y = \sum_{i=1}^D w_i x_i + b + \epsilon = \langle\mathbf{w},\mathbf{x}\rangle + b + \epsilon, \]
which fits a hyperplane through the \((D+1)\)-dimensional point cloud instead of a line through a plane.
With several candidate predictors, a single scatterplot no longer shows the whole picture. A natural first step is a scatterplot matrix: every pair of variables plotted against each other in a grid, making it easy to spot which predictors correlate with the target, and with one another, before fitting anything.
To avoid carrying the bias term \(b\) separately through every derivation, it is folded into the weight vector: define \(\hat{\mathbf{w}}=(\mathbf{w}^\top,b)^\top\) and augment every feature vector with a constant \(1\), \(\hat{\mathbf{x}}=(\mathbf{x}^\top,1)^\top\), so that \(\langle\hat{\mathbf{w}},\hat{\mathbf{x}}\rangle = \langle\mathbf{w},\mathbf{x}\rangle+b\). Stacking all \(n\) augmented feature vectors as rows of a matrix \(\mathbf{X}\) turns the whole data set into a single linear map:
\[ \mathbf{X} = \left(\begin{array}{c}\hat{\mathbf{x}}^{(1)\top}\\\hat{\mathbf{x}}^{(2)\top}\\\vdots\\ \hat{\mathbf{x}}^{(n)\top}\end{array}\right), \qquad \mathbf{y} = \left(\begin{array}{c}y^{(1)}\\y^{(2)}\\\vdots\\y^{(n)}\end{array}\right). \]
The vector of all residuals is \(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}\), and the sum-of-squares loss from before becomes a single quadratic form,
\[ \mathcal{L}(\hat{\mathbf{w}}) = (\mathbf{y}-\mathbf{X}\hat{\mathbf{w}})^\top(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}). \]
Differentiating with respect to \(\hat{\mathbf{w}}\) and setting the result to zero,
\[ \nabla_{\hat{\mathbf{w}}}\mathcal{L} = -2\mathbf{X}^\top\mathbf{y} + 2\mathbf{X}^\top\mathbf{X}\hat{\mathbf{w}} \overset{!}{=} 0, \]
gives the normal equation
\[ \mathbf{X}^\top\mathbf{X}\hat{\mathbf{w}} = \mathbf{X}^\top\mathbf{y}. \]
If the Gramian matrix \(\mathbf{X}^\top\mathbf{X}\) is invertible, this yields the least-squares solution directly:
\[ \hat{\mathbf{w}} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}, \]
with \(\mathbf{X}^\top\mathbf{y}\) called the moment matrix. Setting \(D=1\) and expanding this expression by hand reproduces exactly the \(m\) and \(b\) derived above, so the simple linear case is the \(D=1\) special case of this more general result, not a separate theory.
The normal equation \(\mathbf{X}^\top\mathbf{X}\hat{\mathbf{w}}=\mathbf{X}^\top\mathbf{y}\) can be rewritten as \(\mathbf{X}^\top(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}})=\mathbf{0}\), which says exactly that the residual vector is orthogonal to every column of \(\mathbf{X}\). This is the multivariate version of the orthogonality used above: predictions and residuals remain uncorrelated, and the same variance decomposition \(\operatorname{Var}(y)=\operatorname{Var}(\hat y)+\operatorname{Var}(\epsilon)\) carries over unchanged.
Forming \(\mathbf{X}^\top\mathbf{X}\) explicitly squares the condition number of \(\mathbf{X}\), which can make the normal equation numerically fragile when \(\mathbf{X}\) is ill-conditioned. The singular value decomposition of \(\mathbf{X}\) gives the same \(\hat{\mathbf{w}}\) through the Moore–Penrose pseudo-inverse, \(\hat{\mathbf{w}}=\mathbf{X}^{+}\mathbf{y}\), without ever forming \(\mathbf{X}^\top\mathbf{X}\), and is the numerically preferred route whenever the Gramian is close to singular.
Batch Gradient Descent
The normal equation gives an exact solution, but it requires inverting the \((D+1)\times(D+1)\) matrix \(\mathbf{X}^\top\mathbf{X}\), which costs roughly \(O(D^3)\) arithmetic operations. For a handful of features this is negligible, but once \(D\) reaches into the thousands, or the Gramian is nearly singular, an iterative method that only ever multiplies by \(\mathbf{X}\) is often preferable, or even necessary. Since \(\mathcal{L}\) is a convex quadratic bowl in \(\hat{\mathbf{w}}\), as established above, any method that always steps downhill is guaranteed to converge to the same unique minimum the normal equation finds directly, just approximately and one step at a time.
For this iterative route it is convenient to rescale the loss by \(\frac{1}{2n}\), where \(n\) is again the number of training points. This does not move the minimum, since scaling a function by a positive constant does not change where it is smallest, but it cancels the leading factor of \(2\) out of the gradient derived earlier:
\[ \mathcal{L}(\hat{\mathbf{w}}) = \frac{1}{2n}\bigl\|\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}\bigr\|^2, \qquad \nabla_{\hat{\mathbf{w}}}\mathcal{L} = \frac1n\mathbf{X}^\top\bigl(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y}\bigr). \]
Starting from any initial guess \(\hat{\mathbf{w}}^{(0)}\), batch gradient descent repeatedly takes a small step opposite the gradient, using the entire training set at every single step—hence "batch":
\[ \boxed{ \hat{\mathbf{w}} := \hat{\mathbf{w}} - \frac{\alpha}{n}\,\mathbf{X}^\top\bigl(\mathbf{X}\hat{\mathbf{w}}-\mathbf{y}\bigr) } \]
or, written out for a single weight \(w_j\) using one augmented feature \(\hat{x}_j^{(i)}\) at a time,
\[ w_j := w_j - \alpha\cdot\frac1n\sum_{i=1}^n\Bigl(\langle\hat{\mathbf{w}},\hat{\mathbf{x}}^{(i)}\rangle - y^{(i)}\Bigr)\hat{x}_j^{(i)}. \]
The step size \(\alpha\), the learning rate, controls how far each update moves: too large and the iterates overshoot the minimum and can diverge; too small and convergence is correct but slow. Updates are repeated for a fixed number of iterations, or until \(\mathcal{L}\) stops decreasing by more than some small tolerance between consecutive steps.
In practice, features living on very different scales, say one column in the thousands and another between \(0\) and \(1\), distort the bowl into a long, narrow valley: a small \(\alpha\) is then needed to avoid overshooting along the steep direction, which makes progress along the flat direction painfully slow. Rescaling every feature to zero mean and unit variance before running gradient descent restores a more evenly shaped bowl and usually speeds up convergence considerably. Variants that use only a small random subset of the data per step, stochastic or mini-batch gradient descent, trade a noisier descent direction for a much cheaper update, which matters once \(n\) itself becomes very large.
Gradient descent toward the least-squares solution
The blue reference is the exact closed-form optimum.
Data and current line
Loss by iteration
Loss contours in parameter space
Polynomial Regression and Feature Maps
When the relationship between \(x\) and \(y\) is visibly curved, a straight line underfits no matter how the data is measured. A common fix is a low-degree polynomial,
\[ y = w_3x^3+w_2x^2+w_1x+w_0+\epsilon. \]
This still counts as linear regression, because the model is linear in the parameters \(w_0,\dots,w_3\), even though it is cubic in \(x\). Concretely, define the feature map
\[ \phi(x) = \left(\begin{array}{c}1\\x\\x^2\\x^3\end{array}\right), \]
so that \(y=\langle\mathbf{w},\phi(x)\rangle+\epsilon\) is again exactly the multiple linear regression problem solved above, just with \(\phi(x)\) taking the place of the raw feature vector. Any fixed, nonlinear transformation of \(x\) can be used this way—logarithms, trigonometric terms, interaction terms between several original features—without changing a single line of the normal-equation solution.
The same idea also works in the other direction: instead of, or in addition to, transforming the predictor, the target itself can be transformed to straighten out a curved relationship. Forbes (1857) found that atmospheric pressure and the boiling point of water are related non-linearly, but that \(\log(\text{pressure})\) is very nearly linear in boiling point; fitting the transformed target with the same ordinary least-squares machinery, then undoing the transform on the prediction, often works better than adding polynomial terms to the predictor alone.
When no single global polynomial degree or feature map fits the data well everywhere, an alternative is to give up on a single global model entirely and instead fit many small local ones, one per neighborhood of the data, an idea developed in the chapter on LOESS and LOWESS.
Choosing the polynomial degree, or more generally which features to include, is itself a modeling decision. A degree that is too low underfits, missing real curvature in the data; a degree that is too high overfits, chasing noise instead of signal and generalizing poorly to new points. The degree becomes a hyperparameter, tuned by splitting the data into training, validation, and test sets (or by \(k\)-fold cross-validation): fit on the training set for each candidate degree, compare performance on the validation set, and report the final number only on the untouched test set.
Regularization: Ridge Regression
With many features, few data points, or features that are themselves strongly correlated, \(\mathbf{X}^\top \mathbf{X}\) can become singular or nearly singular, and the normal-equation solution becomes unstable or undefined. Ridge regression addresses this by adding an \(L^2\) penalty on the weights to the loss function,
\[ \mathcal{L}_\lambda(\hat{\mathbf{w}}) = \|\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}\|^2 + \lambda\|\mathbf{w}\|^2, \qquad \lambda>0, \]
which shifts the normal equation to
\[ \hat{\mathbf{w}} = (\mathbf{X}^\top\mathbf{X}+\lambda\mathbf{I})^{-1}\mathbf{X}^\top\mathbf{y}. \]
Adding \(\lambda\mathbf{I}\) makes the matrix being inverted strictly positive definite for any \(\lambda>0\), so it is always invertible, at the cost of shrinking the fitted weights toward zero. This trades a small amount of bias for a often much larger reduction in variance, which tends to generalize better whenever the un-regularized solution would otherwise be dominated by noise.
Common Pitfalls
- Extrapolating the fitted line far outside the range of the training data; the model only describes the relationship where it has seen evidence for it.
- Treating a high \(R^2\) or strong correlation as proof of causation, rather than of a linear association. Specially designed studies, such as a regression discontinuity design that assigns treatment purely by whether a running variable crosses a fixed cutoff, are a notable exception where a causal claim can be justified even from observational data.
- Assuming the regression line is symmetric in \(x\) and \(y\). Minimizing vertical residuals gives the slope \(m_{Y|X}=\operatorname{Cov}(x,y)/\operatorname{Var}(x)\) derived above, but swapping the roles of \(x\) and \(y\) and minimizing horizontal residuals instead gives a different line, with slope \(m_{X|Y}=\operatorname{Cov}(x,y)/\operatorname{Var}(y)\). The two slopes multiply to \(m_{Y|X}\cdot m_{X|Y}=r^2\), so the two lines coincide only in the degenerate case \(r=\pm1\), where every point already lies exactly on one line.
- Forgetting that least squares is sensitive to outliers: a single far-away point can dominate the sum of squared errors and pull the whole line toward it.
- Fitting a straight line to visibly curved data instead of checking the residuals first and adding a feature map or higher polynomial degree where needed.
- Adding more features without regularization or validation, which shrinks the training error but can severely hurt performance on new data through overfitting or a near-singular Gramian matrix.