Linear regression with many features

Two measurements, one number to predict. Below, a flat plane is fitted to forty points in two ways at once. Before it plays: will the step-by-step fit and the one-step formula end on the same plane?

The model and the data matrix

Each example now has $d$ features, $x_{i1},\dots,x_{id}$, collected in the vector $\mathbf{x}_i$, and one target $y_i$, a number. There are $m$ examples; statistics writes $n$ for a sample size, and on these pages $m$ counts examples. The model is a weighted sum of the features plus a constant:

The weights $\mathbf{w}=(w_1,\dots,w_d)$ and the bias term $b$ are the parameters, $\boldsymbol\theta=(\mathbf{w},b)$. With one feature this is the line of Linear regression; with two it is a plane over the floor of the two features, the plane in the picture; with more, the same sum, though no longer something to draw.

The bias term is a weight on a feature that is always 1. So put a 1 into every $\mathbf{x}_i$ as an extra feature and $b$ into $\mathbf{w}$ as its weight; from here on $\mathbf{w}$ holds $b$ too. Stack the examples as the rows of the data matrix $\mathbf{X}$, $m$ rows and $d+1$ columns, one of them all ones, and the targets as $\mathbf{y}=(y_1,\dots,y_m)$. All $m$ predictions are then one product:

The loss and its gradient

The loss is the average squared residual of Models and loss, now over a plane:

Its gradient comes entry by entry, as in Partial derivatives and the gradient. The weight $w_q$ appears in each prediction $\hat y_i$ once, multiplied by $x_{iq}$, so the slope of $(\hat y_i-y_i)^2$ along $w_q$ is $2(\hat y_i-y_i)x_{iq}$. Averaging over the examples and stacking the $d+1$ entries,

The loss is a bowl in the weights, and the gradient is zero at its bottom, where

Two routes to the bottom

In one step. Those are the normal equations, $d+1$ equations in $d+1$ unknowns. Solve them and the answer is exact, at once.

Step by step. Or walk down the bowl by Gradient descent: $\mathbf{w}_{t+1}=\mathbf{w}_t-\eta\nabla J(\mathbf{w}_t)$, with $t$ counting steps. Every step costs one pass through the data.

Now the answer to the question. Both routes end on the same plane: the bowl has one bottom, and the normal equations name it while gradient descent walks to it. In the clip, after 30 steps the two readouts agree to three decimals.

In the clip the two features were centred first: their averages were subtracted, so each feature column adds up to zero. That makes the column of ones square to the two feature columns, so the bias term separates from the weights. One step with $\eta=0.5$ puts $b$ at the average of the targets, 3.109, and from then on the bowl in $(w_1,w_2)$ on the right is the whole story. The readings were made from $y=3+2x_1-x_2$ plus noise, so the fitted weights land near 2 and −1, at 1.975 and −0.886, not on them: forty noisy readings do not pin them exactly.

Which route to take depends on the size of the problem. With a few features, solving the normal equations is quick and exact. With very many features, or data too large to hold at once or arriving all the time, steps are cheaper, and Stochastic gradient descent takes them on a few examples at a time. And many models have no normal equations at all: from Logistic regression on, stepping is the only route.

Reading the weights

Each weight says how much the prediction changes per unit of its feature, with the others held still: $w_1$ is in units of the target per unit of $x_1$. So the size of a weight depends on the units. Measure a distance in kilometres instead of metres and its weight becomes 1000 times larger, for exactly the same model. A bigger weight does not mean a more important feature.

Features that move together are a second trap. If two columns of $\mathbf{X}$ are equal, the model can raise one weight and lower the other by the same amount without changing a single prediction: there are infinitely many best answers, and $\mathbf{X}^\top\mathbf{X}$ cannot be inverted. If two columns are only nearly equal, the bowl becomes a long, flat-bottomed valley:

The two features in the valley differ only in their middle two readings, yet the weights they lead to, −5.75 and 8.75, are far from the round bowl’s 1.25 and 1.75, and of opposite signs. Along the valley floor the loss hardly changes, so the data cannot tell the two weights apart: a little new noise slides the bottom far along the valley, and the weights swing while the predictions barely move. Conditioning measures how long the valley is. On such a bowl, gradient descent zig-zags as it did on the stretched bowl of Gradient descent.

Does another feature help?

Add a feature to a least-squares fit and the training error can only fall or stay: the fit could always set the new weight to 0 and do as well as before. That is exactly why the training error is a poor judge of a feature. Here is a feature of pure noise, a random number with nothing to do with the target, added to the two real ones.

The noise feature helped on the training readings in every one of the hundred repeats, because the fit found some chance pattern in it to use. On new readings that pattern is not there, and the error rose, in 67 of the 100 repeats and on average. Bias and variance of a model explains the size of both changes: each fitted coefficient takes about $\sigma^2/m$ off the training error and adds about as much to the error on new data, where $\sigma^2$ is the variance of the noise, here $0.3^2=0.09$, and $0.09/40=0.00225$.

Noise and least squares

Why squares? Suppose each target is the model’s value plus noise, $y_i=\mathbf{w}^\top\mathbf{x}_i+\varepsilon_i$, where the $\varepsilon_i$ are independent, bell-shaped, with mean 0 and the same spread. Then the weights under which the readings are most likely are exactly the least-squares weights.

Worked example

Four examples $(x_1,x_2,y)$: $(1,1,6)$, $(1,-1,2)$, $(-1,1,3)$ and $(-1,-1,0)$. With columns ones, $x_1$, $x_2$, and $\mathbf{w}=(b,w_1,w_2)$:

The normal equations. The three columns are at right angles to each other, and each has squared length 4, so $\mathbf{X}^\top\mathbf{X}=4\mathbf{I}$. And $\mathbf{X}^\top\mathbf{y}=(6+2+3+0,\ 6+2-3-0,\ 6-2+3-0)=(11,5,7)$. So $4\mathbf{w}=(11,5,7)$:

Check. The predictions are $2.75+1.25+1.75=5.75$, then 2.25, 3.25 and −0.25. The residuals are $0.25$, $-0.25$, $-0.25$ and $0.25$, so the loss is $0.0625$, down from $(36+4+9+0)/4=12.25$ with every weight zero.

Gradient descent. At zero, $\nabla J=\frac24\mathbf{X}^\top(\mathbf{0}-\mathbf{y})=-\tfrac12(11,5,7)=(-5.5,-2.5,-3.5)$. One step with $\eta=0.25$ goes to $(1.375,0.625,0.875)$, exactly halfway. One step with $\eta=0.5$ lands on $(2.75,1.25,1.75)$, the answer. That works because this bowl is round: its curvature is $\frac2m\mathbf{X}^\top\mathbf{X}=2\mathbf{I}$, 2 in every direction, so every step multiplies the distance to the bottom by $1-2\eta$, which is 0 at $\eta=0.5$. Any $\eta<1$ converges; at $\eta=1$ it bounces for ever.

A prediction. At $(x_1,x_2)=(0,0)$ the model predicts the bias term, 2.75.

Units. If $x_1$ were in kilometres rather than metres, $w_1$ would be $1.25\times1000=1250$ per kilometre for the same predictions.

Where you’ll meet this

Energy use. A building’s daily energy use predicted from the outdoor temperature and the number of people inside, one weight each.

Sensor correction. A sensor whose reading drifts with temperature and supply voltage is corrected by a weighted sum of both at once.

Filters. The Wiener filter and linear prediction choose filter taps by the same normal equations, with past samples as the features.

Reference card