Chapter 2: Linear Regression

Without a model you can only repeat values you have already observed. A model is a formula with unknown numbers. Here $x$ is the single input feature and $\hat{y}$ — "y-hat" — is the predicted output, hatted to separate it from an observed $y_i$. The parameters $w_0$ and $w_1$ are what we choose: the intercept $w_0$ is the prediction at $x = 0$, the slope $w_1$ is how far $\hat{y}$ moves per one-unit increase in $x$. This is chapter 0's $y = mx + b$ in ML notation. Its assumption is strong: if the truth curves, no $(w_0, w_1)$ repairs it and the leftover errors are systematic, not random.

A model can only be fitted if you can measure how wrong it is at each point; the residual is that per-point measurement. For training example $i$, $y_i$ is the observed target, $\hat{y}_i = w_0 + w_1 x_i$ is what the line predicts at that example's input, and $e_i$ is the signed vertical gap: positive means the model underpredicted, negative means it overpredicted. Residuals are not properties of the data: move the line and every $e_i$ changes, which is what fitting exploits. They measure vertical error only, so $x$ is assumed noise-free; residuals that curve or fan out warn that the linear model above is the wrong shape.

Residuals give $n$ numbers, but choosing parameters needs one number to minimise — and simply adding residuals fails, because positive and negative misses cancel and a terrible line can total zero. Squaring destroys the signs, so the mean of the squared residuals over the $n$ training examples is zero only for a perfect fit. That number is the loss. Dividing by $n$ makes it per-example, so datasets of different sizes stay comparable. Its units are the square of $y$'s units, which is why people report $\sqrt{\text{MSE}}$ (RMSE) instead. The price of squaring is outlier sensitivity: one point ten units off costs as much as a hundred points one unit off.

Squaring residuals is a choice, not a law, with three defences. Probabilistic: assume every observation is the line plus independent noise $\epsilon_i \sim \mathcal{N}(0, \sigma^2)$, where $\sigma^2$ is the noise variance. Each point's likelihood is proportional to $e^{-e_i^2 / 2\sigma^2}$, the dataset's likelihood is their product, and taking the logarithm turns that product into $-\frac{1}{2\sigma^2}\sum e_i^2$ plus constants — so maximising likelihood is literally minimising MSE. Computational: $e^2$ is differentiable everywhere, while $|e|$ has a corner at zero. Behavioural: doubling an error quadruples its cost. That last is also the flaw — MSE estimates the conditional mean, so one outlier drags the line, whereas MAE estimates the median and shrugs.

Because MSE is a convex bowl in $(w_0, w_1)$, its one flat point is the minimum: set both partial derivatives to zero and solve. Differentiating $J = \frac{1}{n}\sum(y_i - w_0 - w_1 x_i)^2$ with respect to $w_0$ gives $-\frac{2}{n}\sum e_i = 0$ — residuals must sum to zero — which rearranges to $w_0 = \bar{y} - w_1 \bar{x}$, where $\bar{x}$ and $\bar{y}$ are sample means. Substituting back and differentiating with respect to $w_1$ produces the formula shown: the covariance over the variance of $x$. It fails only when $\sum(x_i - \bar{x})^2 = 0$ — every $x_i$ identical, so there is no slope to estimate; near-constant $x$ makes it numerically fragile.

MSE answers "how wrong?" in squared units of $y$, so its value alone cannot say whether a fit is good. $R^2$ fixes that by comparing against the laziest model: ignore $x$ and always predict $\bar{y}$. The denominator $\sum(y_i - \bar{y})^2$ is that baseline's squared error — the total variation in $y$; the numerator $\sum(y_i - \hat{y}_i)^2$ is what your line leaves. Their ratio is the fraction of baseline error still unexplained, so $R^2$ is the fraction removed: 0.8 means four fifths of the variation in $y$ is explained by $x$. It goes negative when a model beats nothing, and never falls when you add a feature, even a useless one.

The loss at your current parameters says nothing about which way to move, and a blind grid search over $K$ values per axis costs $K^p$ evaluations. The gradient answers in one calculation. $J$ is the scalar loss, $\theta$ the parameter vector, and $\nabla J(\theta)$ a vector shaped like $\theta$ whose $j$-th entry $\partial J/\partial\theta_j$ is how fast loss rises per unit increase in $\theta_j$, other parameters fixed. Assembled, those slopes point along steepest ascent -- precisely why descent subtracts them. It assumes $J$ is differentiable (mean absolute error is not, at zero), and on a flat plateau it sits near zero, offering no direction. These are chapter 0's partial derivatives doing real work.

A direction is not a move: the gradient is a rate of change, not a distance, and nothing stops you walking past the minimum. This rule supplies the missing scale. $\theta_{old}$ is the current parameter vector, $\nabla J(\theta)$ the gradient at that point (recompute it every iteration), and $\alpha > 0$ the learning rate; the minus sign turns steepest ascent into descent. It is a first-order Taylor step: we replace $J$ by its tangent plane and trust it only briefly, which is why steps stay small. The gradient shrinks near a minimum, so steps shrink automatically. It fails when $\alpha$ outruns the local curvature. This same line, fed by backpropagation, trains every neural network in chapter 9.

Rather than assert this formula, derive it. Start from $J(w) = \frac{1}{n}\sum_i (y_i - \hat{y}_i)^2$ with prediction $\hat{y}_i = wx_i$. The power rule gives the factor 2; the chain rule then multiplies by $\partial\hat{y}_i/\partial w = x_i$ and carries a minus sign because $\hat{y}_i$ enters subtracted. So: 2 from the square, $x_i$ from the chain rule, $\frac{1}{n}$ from the MSE's averaging -- drop it and the gradient grows with dataset size, forcing a new $\alpha$ for every $n$. Here $n$ counts examples, $y_i - \hat{y}_i$ is the residual from the previous topic, $x_i$ the feature. In matrix form, $\frac{2}{n}X^T(Xw - y)$. It is a residual-feature correlation that vanishes when the two are uncorrelated -- the normal equations again.

The gradient fixes direction; $\alpha$ is the free knob you choose, and it decides whether training works. Too small (say $\alpha = 0.001$) still converges, but the curve is a nearly straight decline needing a hundred times more epochs than you will sit through. Too large overshoots the valley and lands higher up the far wall, so error grows by a constant factor per step: loss reads 12, 40, 300, 1e6, then inf, then nan once inf - inf reaches the update -- and nan weights never recover. Just right drops the loss steeply, then flattens. Practically: standardize features, sweep 0.001 to 0.1 logarithmically, keep the largest value whose loss still falls monotonically.

Taken literally, the update rule averages the gradient over all $n$ training examples, giving the exact gradient of the loss rather than an estimate. One pass over the data -- one epoch -- yields exactly one update, costing $O(nd)$ for $n$ examples, $d$ features. Nothing is sampled, so the trajectory is deterministic and with a sane $\alpha$ the loss falls smoothly and monotonically -- the easiest variant to debug. It breaks on scale: with $n$ in the millions, a whole pass buys one step, and $X$ must fit in memory. A few hundred updates is a short journey, so wall-clock convergence is far slower than its clean curve suggests.

Batch GD wastes work: after a thousand similar examples the direction is clear, yet it reads the remaining 999,000 before moving. SGD takes the opposite extreme: draw one random example $i$, compute $\nabla J_i$ from it, update, repeat -- so one epoch delivers $n$ updates at $O(d)$ each. The estimate is unbiased, $E[\nabla J_i] = \nabla J$, so on average it moves the right way, but variance is large and individual steps point wrong constantly. The loss curve is jagged, and at fixed $\alpha$ parameters never settle but rattle inside a noise ball whose radius scales with $\alpha$. Decaying $\alpha$ shrinks it. Shuffle every epoch, or ordering becomes systematic bias.

Neither extreme is what anyone runs. Mini-batch averages the gradient over $B$ randomly drawn examples: $B = n$ recovers batch GD, $B = 1$ recovers SGD, one dial spanning both. Averaging $B$ independent estimates divides gradient variance by $B$, so noise falls only as $\sqrt{B}$ -- diminishing returns past a hundred or so. The usual 32 to 256 dominates for hardware reasons, not statistical ones: a batch is one $B \times d$ matrix multiply, and vectorized CPU or GPU units chew through 32 rows in nearly the time of one, so that variance reduction is near-free. Push $B$ into the thousands and gradients get so quiet that $\alpha$ must be raised and warmed up.

"Stop when the gradient is zero" is useless in practice: floating-point arithmetic rarely lands on zero, and under SGD the gradient is noisy and never settles there. Real criteria are three, usually combined. A gradient norm test, $\|\nabla J\| < \epsilon$, scale-dependent and meaningless unless features are standardized. A relative loss test, $|J_t - J_{t-1}| / |J_{t-1}| < \text{tol}$ over several consecutive epochs -- scikit-learn's tol with n_iter_no_change is exactly this. And a hard iteration cap, so a diverging run still terminates. Better still, early-stop on a validation loss with patience. Always plot the curve: flat from step one means $\alpha$ too small or features unscaled; rising or nan means too large.

With one predictor you can only ask how price moves with area, never how it moves at a fixed bedroom count. Multiple regression fixes that by giving every feature its own weight. Here $x_j$ is the $j$-th of $d$ features, $w_j$ is its weight, and $w_0$ is the intercept, the prediction when all features are zero. Defining $\mathbf{x} = [1, x_1, \ldots, x_d]^T$ folds the intercept in, so the model is one dot product $\mathbf{w}^T\mathbf{x}$. Each $w_j$ is an all-else-equal slope, meaningful only if the others can be held fixed — which fails when predictors move together. The objective remains the previous topic's mean squared error.

Computing $\mathbf{w}^T\mathbf{x}_i$ in a Python loop over $n$ examples is slow and hides the structure; stacking them turns the dataset into one multiplication. Row $i$ of the design matrix $\mathbf{X}$ is example $i$, so $\mathbf{X}$ is $n \times d$: $n$ rows, one per example, $d$ columns, one per feature. To absorb the intercept you prepend a column of ones, making $\mathbf{X}$ $n \times (d+1)$ and $\mathbf{w}$ $(d+1) \times 1$ with $w_0$ first, so each row contributes $w_0 \cdot 1$. Then $\mathbf{X}\mathbf{w}$ is $n \times 1$, matching $\mathbf{y}$. Mismatched shapes are the usual bug; worse, omitting the ones column raises no error and silently forces the fit through the origin.

Gradient descent walks downhill to the minimum; for squared error you can jump straight to the bottom. Setting the gradient of $J = \|\mathbf{y} - \mathbf{X}\mathbf{w}\|^2$ to zero gives $\mathbf{X}^T\mathbf{X}\mathbf{w} = \mathbf{X}^T\mathbf{y}$, hence the closed form. $\mathbf{X}^T\mathbf{X}$ is the $(d+1) \times (d+1)$ Gram matrix and $\mathbf{X}^T\mathbf{y}$ is $(d+1) \times 1$, so $\mathbf{w}$ is too: no learning rate, no iterations, no convergence check. The inverse needs full column rank, so it fails when a feature is an exact combination of others, or when $d+1 > n$. Inversion costs $O(d^3)$, which is why past a few thousand features gradient descent's $O(nd)$ per step wins.

Gradient descent zigzagged in the previous topic because one feature ranged over thousands and another over single digits, stretching the loss contours into a narrow valley. Standardization removes the units: subtract the feature's mean $\mu$ and divide by its standard deviation $\sigma$, column by column, leaving every feature with mean 0 and variance 1. A scaled value reads as a number of standard deviations from the mean, which makes weights comparable in the next concept. Prefer it for roughly bell-shaped or unbounded features. The non-negotiable rule: fit $\mu$ and $\sigma$ on the training set only and reuse those numbers on test data; recomputing them there leaks information.

Standardization leaves values unbounded, which is wrong when a model expects a fixed input range, or a feature has hard limits like pixel intensity 0 to 255. Min-max normalization maps the training range onto $[0, 1]$: subtract the column minimum $x_{min}$, divide by the range $x_{max} - x_{min}$, so the smallest training value becomes 0 and the largest 1. It also leaves exact zeros at zero when the minimum is 0, which matters for sparse data. Its weakness: both endpoints are single observations, so one outlier stretches the denominator and squashes every real value into a thin band — standardize instead. Again both statistics come from training data only.

Once fitted, you want to know which features drive the prediction. The sign of $w_j$ gives direction, and on standardized features $|w_j|$ gives the effect of a one-standard-deviation change, so magnitudes can be ranked. On raw features they cannot: the same measurement in millimetres rather than metres carries a weight 1000 times smaller with identical predictive power, so a large coefficient may only mean a small-numbered unit. Two caveats. A weight carries estimation uncertainty, so weigh it against its standard error (the $t$-statistic) first. And every weight is conditional on the other features present, which is what makes it unreliable under the collinearity described next.

Two features carrying nearly the same information, like living area and total area, leave no way to divide credit between them. Their columns are almost linearly dependent, so $\mathbf{X}^T\mathbf{X}$ has a near-zero eigenvalue, $(\mathbf{X}^T\mathbf{X})^{-1}$ blows up, and $\text{Var}(\mathbf{w}) = \sigma^2 (\mathbf{X}^T\mathbf{X})^{-1}$ inflates. The symptom is distinctive: huge coefficients with huge standard errors that flip sign when you refit on a slightly different sample, while $R^2$ stays good. It damages interpretation, not accuracy, so it hides unless you check a variance inflation factor (above 5 to 10 is trouble). Fix it by dropping or merging columns, or by adding a ridge penalty making $\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I}$ invertible — topic 4.

A straight line has the same slope everywhere, so it cannot bend toward the data in one region without pulling away in another. Polynomial regression buys curvature by adding powers of the input as columns. Here $\hat{y}$ is the prediction, $x$ the input feature, $d$ the degree, and $w_0, \ldots, w_d$ the $d+1$ parameters — one per power, $w_0$ being the intercept. Each $x^k$ is treated as its own feature, so the design matrix $\mathbf{X}$ from Matrix Form gains columns $[1, x, \ldots, x^d]$ and the Normal Equation recovers $\mathbf{w}$ unchanged. The catch: past the training range the highest power dominates and predictions diverge, so extrapolation is unusable.

Powers of $x$ are only one option; the general move is to keep the model linear and put the flexibility into a fixed map $\phi$ applied before fitting. $\phi$ takes a $p$-dimensional input row and returns a $D$-dimensional one with $D > p$; the estimator learns one weight per transformed column, so a curve in $x$-space is a flat hyperplane in $\phi$-space. Nothing forces $\phi$ to be polynomial — $\log x$ and $\sqrt{x}$ are valid too. The failure mode is scale: $x \in [0,10]$ gives $x^3 \in [0,1000]$, which stalls gradient descent and ill-conditions $\mathbf{X}^T\mathbf{X}$, so apply Feature Scaling: Standardization after expanding, and apply the identical $\phi$ at test time.

With two inputs and no product term the model is additive: $\hat{y} = w_0 + w_1x_1 + w_2x_2$ is a flat plane, and the return on $x_1$ is $w_1$ whatever $x_2$ does. Adding the column $x_1x_2$ with weight $w_3$ tilts that plane, because the marginal effect becomes $\partial\hat{y}/\partial x_1 = w_1 + w_3x_2$ — a slope that is itself a function of the other feature. Positive $w_3$ means the features reinforce each other, negative $w_3$ means one saturates the other. The danger is Multicollinearity: the product column correlates strongly with $x_1$ and $x_2$ unless both are centered, inflating the weight variances. Centering before forming the product largely removes it.

Nothing picks $d$ for you: it is the largest lever on how the model behaves. Because a degree-$d$ model contains every degree-$(d-1)$ model — set the top weight to zero — training Mean Squared Error only falls as $d$ rises, so it is worthless here. Validation error, measured on rows the fit never saw, traces a U instead: high at $d = 1$ where the model is too rigid to follow the curve (underfitting), lowest near the data's true complexity, then climbing as spare powers fit noise (overfitting). The widening gap between the two curves is the signal. Degrees 2–3 usually suffice. Chapter 4 develops this machinery properly.

With one input, degree 5 costs six columns and nobody notices. With many inputs the same degree is ruinous. A degree-$d$ expansion of $p$ features produces $\binom{p+d}{d}$ columns — one per multiset of at most $d$ factors drawn from the $p$ features. For fixed $d$ that grows like $p^d / d!$, so each extra degree multiplies the column count by roughly $p/d$. Once $D$ approaches the sample count $n$, $\mathbf{X}^T\mathbf{X}$ becomes singular, the Normal Equation has no unique solution, and weight variances explode well before that. The rule of thumb is 10–20 rows per parameter. The fixes are interaction_only=True, feature selection, or the ridge penalty from the next topic.

Writing the expansion by hand is tedious and easy to get wrong between training and inference: a different column order at predict time pairs each weight with the wrong feature. PolynomialFeatures is a transformer — fit records the column layout for the given degree and transform replays it on later rows. Two arguments matter. include_bias=False drops the leading column of ones, which you want when the estimator fits its own intercept, since two constant columns make $\mathbf{X}^T\mathbf{X}$ singular. interaction_only=True keeps cross-terms like $x_1x_2$ but drops pure powers, blunting the Curse of Dimensionality. Chain it inside a Pipeline with Feature Scaling: Standardization so the scaler sees the expanded columns.

Ordinary least squares cares about error, not about the size of a coefficient, so with two nearly collinear features it happily pairs a huge $+w_1$ against a huge $-w_2$ that almost cancel — the multicollinearity failure from Multiple Regression. Ridge attaches a price tag. $J(\mathbf{w})$ is the objective minimised, $\mathbf{w} = (w_1, \ldots, w_p)$ holds the $p$ coefficients, MSE is the mean squared error, and $\lambda \ge 0$ is the penalty strength. Because the cost is quadratic, a weight of 10 costs a hundred times a weight of 1, so the optimiser prefers many small weights to one large one. Ridge shrinks every coefficient but zeroes none.

Ridge shrinks but never eliminates: with a thousand candidate features you still carry a thousand nonzero coefficients and no idea which matter. Lasso penalises $\sum_{j=1}^{p} |w_j|$, the L1 norm, where $w_j$ is feature $j$'s coefficient, $p$ the feature count, and $\lambda$ the strength. Absolute value has a kink at zero and no derivative there, so there is no Normal-Equation-style closed form; solvers use coordinate descent, soft-thresholding one coefficient at a time. The payoff is exact zeros — feature selection done by the fit itself, making Feature Importance from Multiple Regression automatic. Caveats: among correlated features Lasso keeps one arbitrarily, and when $p > n$ it selects at most $n$.

Both penalties shrink, so why does only one land on exact zeros? Reframe the fit as constrained optimisation: minimise MSE subject to $\sum_j |w_j| \le t$ (L1) or $\sum_j w_j^2 \le t$ (L2), the budget $t$ tightening as $\lambda$ grows. In two dimensions the L1 region is a diamond with corners on the axes; the L2 region is a circle with none. The MSE contours are ellipses centred on the OLS solution from Simple Linear Regression, expanding until they touch the region. A sharp corner catches an expanding ellipse from many directions, and a corner is exactly where a coefficient is zero. A circle touches an axis only by chance.

Lasso's arbitrary choice among correlated features is a real problem: with two near-duplicate sensors, resampling noise picks the survivor, and a report built on it is not reproducible. Elastic Net charges both. $\lambda$ sets the overall strength and the mixing parameter $\alpha \in [0,1]$ divides it — $\alpha = 1$ is pure Lasso, $\alpha = 0$ is pure Ridge, in between both. The L1 part zeroes weak features; the L2 part supplies a grouping effect that makes correlated features share weight instead of compete, so a group enters or leaves together. The costs: a second hyperparameter, and a naming trap — scikit-learn calls the strength alpha and this ratio l1_ratio.

Every formula above hinges on one number you supply. $\lambda \ge 0$ multiplies the penalty — the exchange rate between fitting the data and keeping coefficients small. At $\lambda = 0$ the penalty vanishes and you recover ordinary least squares — minimum bias, maximum variance, free to overfit. As $\lambda \to \infty$ every coefficient is crushed to zero and the model predicts the training mean $\bar{y}$ — zero variance, maximum bias. Useful values sit between, and you do not pick one by eye: sweep a log-spaced grid and keep the best cross-validated error, as RidgeCV and LassoCV do. Judging $\lambda$ on training error always picks zero, as in Degree Selection.

The penalty is computed on raw coefficient values, which makes it scale-dependent, unlike plain OLS. Halve a feature's units and its coefficient must double for the same prediction; the fit is unchanged, but the L2 charge on it quadruples. Record a length in kilometres, not millimetres, and its coefficient grows a millionfold, so its L2 penalty grows by $10^{12}$ — and the model effectively stops using it. The fix is the standardization from Multiple Regression: subtract each column's mean and divide by its standard deviation, so one $\lambda$ means the same thing everywhere. One exception: the intercept $w_0$ goes unpenalised, since shrinking it only drags predictions off the data's centre.