The Normal Equation: Solving Regression in One Shot

The previous post left a hook: you can solve linear regression in a single computation, no gradient descent needed. That's exactly what lectures 2b and 2c show, and along the way the professor compares the result against five other very different kinds of regressor.

Zeroing the derivative instead of walking toward it

Gradient descent works because, at every step, it walks a little in the direction that reduces the error. But for linear regression, the error function (the same sum of squares I'd already been using) has a special property: it's a parabola with respect to the weights, a smooth bowl with no false valleys. And a smooth bowl has exactly one point where the derivative is zero, its bottom. Instead of walking there step by step, you can compute that point directly.

Bishop writes the linear model in a more general form than mine, with a matrix Φ (the design matrix, one row per patient, one column per input variable, plus a column of 1's for the bias term). The error function is

$E_D(\mathbf{w}) = \frac{1}{2}\sum_{n=1}^{N}\{t_n - \mathbf{w}^{\mathsf{T}}\boldsymbol{\phi}(\mathbf{x}_n)\}^2$

Setting this function's gradient to zero and solving for w leaves

$\mathbf{w}_{\text{ML}} = (\boldsymbol{\Phi}^{\mathsf{T}}\boldsymbol{\Phi})^{-1}\boldsymbol{\Phi}^{\mathsf{T}}\mathbf{t}$

This is the normal equation. One matrix, two multiplications, and an inversion, and the optimal fit falls right out, no learning rate, no choosing a number of iterations, none of the care gradient descent demands.

The trick of folding bias in

Notice the formula above has no separate + b, the bias lives inside w itself (Bishop calls it $w_0$). That only works because he defines a fake "basis function," $\phi_0(\mathbf{x}) = 1$, an entire column of 1's, so multiplying it by $w_0$ gives exactly $w_0$ for every patient, the same effect as adding a fixed bias.

The code does this by hand, gluing a column of 1's onto the front of X:

def include_bias(X):
    return np.hstack((np.ones((X.shape[0], 1)), X))

That changes the class's shape: instead of storing coefs_ and intercept_ as two separate things (like in the previous post), there's now a single vector w_, where w_[0] is the bias and the rest are each variable's coefficient.

Three versions, one real comparison

The notebook runs three variations of the same LinearRegressor class, all on the same diabetes dataset (now using scikit-learn's real train_test_split, so the numbers don't match the previous post's exactly, which used my own hand-rolled split):

The normal equation beats both, without me having to pick a learning_rate or max_iter. That's not a coincidence: gradient descent is a way to approximate this exact same answer by iterating, and with enough iterations and a good learning rate it converges to the same place. The normal equation just skips straight to the end.

The pseudo-inverse is pinv

$(\Phi^{\mathsf{T}}\Phi)^{-1}\Phi^{\mathsf{T}}$ has its own name: the Moore-Penrose pseudo-inverse, denoted $\Phi^\dagger$. It's a generalization of "matrix inverse" to matrices that aren't square (which is always the case here: Φ has one row per patient and one column per variable, almost never equal). np.linalg.pinv(X) computes exactly this, so the whole class shrinks to:

class LinearRegressor(BaseEstimator, RegressorMixin):
    def fit(self, X, y):
        X = include_bias(X)
        self.w_ = np.linalg.pinv(X) @ y
        return self

    def predict(self, X):
        X = include_bias(X)
        return (X @ self.w_).reshape(X.shape[0],)

An entire linear regression class, solved in one line. And Bishop already flags the risk: if two input columns are too similar to each other (collinearity), $\Phi^{\mathsf{T}}\Phi$ gets close to singular and the computation becomes numerically unstable. That'll come back once the course reaches feature selection.

Matches scikit-learn to the decimal

The proof the math is right: I ran the same pseudo-inverse approach and scikit-learn's ready-made LinearRegression() side by side.

regressor = LinearRegressor()  # the pinv-based class above
regressor.fit(X_train, y_train)
print(mean_squared_error(y_train, regressor.predict(X_train)))

from sklearn.linear_model import LinearRegression
sk_regressor = LinearRegression().fit(X_train, y_train)
print(mean_squared_error(y_train, sk_regressor.predict(X_train)))

Output: my version: train MSE 2868.5497028355776. scikit-learn's LinearRegression: train MSE 2868.549702835577.

The difference only shows up in the last decimal place, floating-point noise, not a difference in method. LinearRegression() computes exactly this under the hood.

The regressor zoo

Lecture 2c takes this same normal-equation regressor and also measures MSE on the test set (the 20% held out from training), then compares it against five very different kinds of model, all with scikit-learn's default parameters:

I'm not going to explain how each of these models works internally yet (KNN, decision trees, and random forests each get their own lecture later in the course, and that's where I'll come back to them properly). But three lessons come out of the table alone:

1. Low training MSE means nothing on its own. The decision tree zeroed out its training error (it literally memorized every patient) and was the worst of all of them on the test set. That's overfitting in its purest form, the same phenomenon Bishop showed back in chapter 1 with the degree-9 polynomial.
2. The test-set winner wasn't the model that fit training best. The random forest fits training almost as well as the full tree (2530.82, quite close to the tree's absurd zero), but without going as far, and that's exactly why it generalizes better: 2785.98 on test, the lowest test MSE in the whole table. A forest is many trees trained on different slices of the data, with the final prediction being the average across all of them, and that average cancels out a lot of the excess each individual tree commits.
3. SGDRegressor beat the exact normal equation on the test set (2863.35 versus 2900.19), even with a slightly worse training MSE. That's not a coincidence: as I already saw in the specialization playlist, scikit-learn's SGDRegressor ships with L2 regularization on by default. Here that regularization, without me asking for it, ended up helping it generalize a bit better.

LinearSVR came out visibly worse than everything else, on both train and test, but that's more about its default hyperparameters not suiting this dataset than about the method itself, a story for another day.

Wrapping up

Practical application

I use the normal-equation model (the same pinv-based LinearRegressor from above) and look at the test set, the 89 patients the model never saw during fitting, to visualize what that 2900.19 MSE actually means case by case.

regressor = LinearRegressor()
regressor.fit(X_train, y_train)
y_pred_test = regressor.predict(X_test)

Notice the cloud looks a lot like the one from the previous post (which was on training data), without getting visibly worse on test. That confirms numerically what the table already showed: train MSE 2868.55 versus test MSE 2900.19, a small gap. The model didn't memorize training, it genuinely generalized, it's just that being a plain linear model, it still misses quite a bit case by case, the same ceiling I already saw in the previous post.