Deriving the Normal Equation (optional)

This page shows where the normal equation in Linear Regression comes from:

\mathbf{w}^* = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}

The idea is simple: MSE is shaped like a bowl, and at the bottom of a bowl the slope is zero. So we write MSE with matrices, take its gradient (the slope in every direction at once), set the gradient to zero, and solve for $\mathbf{w}$ .

Step 1: Write MSE with Matrices

As in the lecture, we absorb the bias into the weights by adding a column of 1s to $\mathbf{X}$ , so $b$ is the first entry of $\mathbf{w}$ . Then the predictions for all $n$ instances are $\mathbf{X}\mathbf{w}$ , and the errors are the vector $\mathbf{X}\mathbf{w} - \mathbf{y}$ .

MSE is the average of the squared errors. The sum of squares of a vector's entries is the vector's dot product with itself, so:

MSE(\mathbf{w}) = \frac{1}{n} \sum_{i=1}^{n} (\mathbf{w}^\top \mathbf{x}^{(i)} - y^{(i)})^2 = \frac{1}{n} (\mathbf{X}\mathbf{w} - \mathbf{y})^\top (\mathbf{X}\mathbf{w} - \mathbf{y})

Step 2: Expand

Multiply out the product, using the transpose rule $(\mathbf{X}\mathbf{w})^\top = \mathbf{w}^\top \mathbf{X}^\top$ :

(\mathbf{X}\mathbf{w} - \mathbf{y})^\top (\mathbf{X}\mathbf{w} - \mathbf{y}) = \mathbf{w}^\top \mathbf{X}^\top \mathbf{X} \mathbf{w} - \mathbf{w}^\top \mathbf{X}^\top \mathbf{y} - \mathbf{y}^\top \mathbf{X} \mathbf{w} + \mathbf{y}^\top \mathbf{y}

The two middle terms are equal. Each is a single number (a $1 \times 1$ matrix), and the second is the transpose of the first; a single number is its own transpose. So they combine:

MSE(\mathbf{w}) = \frac{1}{n} \left( \mathbf{w}^\top \mathbf{X}^\top \mathbf{X} \mathbf{w} - 2 \mathbf{w}^\top \mathbf{X}^\top \mathbf{y} + \mathbf{y}^\top \mathbf{y} \right)

Step 3: Take the Gradient

The gradient $\nabla_{\mathbf{w}}$ is the vector of partial derivatives, one for each entry of $\mathbf{w}$ . We need two rules, the vector versions of $\frac{d}{dw}(aw) = a$ and $\frac{d}{dw}(aw^2) = 2aw$ :

$\nabla_{\mathbf{w}} (\mathbf{w}^\top \mathbf{a}) = \mathbf{a}$ for a constant vector $\mathbf{a}$
$\nabla_{\mathbf{w}} (\mathbf{w}^\top \mathbf{A} \mathbf{w}) = 2\mathbf{A}\mathbf{w}$ for a constant symmetric matrix $\mathbf{A}$ (one that equals its own transpose)

You can check the first rule with two entries: $\mathbf{w}^\top \mathbf{a} = w_1 a_1 + w_2 a_2$ , whose partial derivatives are $a_1$ and $a_2$ , which is the vector $\mathbf{a}$ .

$\mathbf{X}^\top \mathbf{X}$ is symmetric, because $(\mathbf{X}^\top \mathbf{X})^\top = \mathbf{X}^\top \mathbf{X}$ . The last term $\mathbf{y}^\top \mathbf{y}$ does not contain $\mathbf{w}$ , so its derivative is zero. Applying the two rules term by term:

\nabla_{\mathbf{w}} MSE(\mathbf{w}) = \frac{1}{n} \left( 2 \mathbf{X}^\top \mathbf{X} \mathbf{w} - 2 \mathbf{X}^\top \mathbf{y} \right) = \frac{2}{n} \mathbf{X}^\top (\mathbf{X}\mathbf{w} - \mathbf{y})

This is the same gradient that Gradient Descent uses to take each step downhill.

Step 4: Set the Gradient to Zero and Solve

At the bottom of the bowl, the gradient is zero:

\frac{2}{n} \mathbf{X}^\top (\mathbf{X}\mathbf{w} - \mathbf{y}) = \mathbf{0} \quad \Longrightarrow \quad \mathbf{X}^\top \mathbf{X} \mathbf{w} = \mathbf{X}^\top \mathbf{y}

These are the normal equations. Multiplying both sides by the inverse $(\mathbf{X}^\top \mathbf{X})^{-1}$ gives the normal equation:

\mathbf{w}^* = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}

Two notes:

Why "normal"? The equation $\mathbf{X}^\top (\mathbf{y} - \mathbf{X}\mathbf{w}) = \mathbf{0}$ says the errors are perpendicular to every column of $\mathbf{X}$ , and another word for perpendicular is normal.
When it fails. $(\mathbf{X}^\top \mathbf{X})^{-1}$ exists only when no feature is a combination of the others. If a feature is duplicated, or one feature is another times a constant (say, price in dollars and price in cents), $\mathbf{X}^\top \mathbf{X}$ has no inverse and the formula breaks. scikit-learn uses SVD instead, which still finds a solution in that case; see When the Inverse Fails below.

Check It in NumPy

On the simulated data from the lecture, the normal equation finds $b = 5.215$ and $w_1 = 2.770$ , and the gradient of MSE at that point is zero up to rounding error:

import numpy as np

np.random.seed(42)  # the same simulated data as the Linear Regression lecture
X = 2 * np.random.rand(100, 1)
y = 5 + 3 * X + np.random.randn(100, 1)
X_with_bias = np.c_[np.ones((100, 1)), X]  # add a column of 1s for the bias b
n = len(X_with_bias)

# the normal equation
w_best = np.linalg.inv(X_with_bias.T @ X_with_bias) @ X_with_bias.T @ y
print(w_best)
print('*' * 50)

# the gradient of MSE at w_best is (almost exactly) zero: this is the bottom of the bowl
gradient = 2 / n * X_with_bias.T @ (X_with_bias @ w_best - y)
print(gradient)

Output:
[[5.21509616]
 [2.77011339]]
**************************************************
[[ 1.08246745e-15]
 [-9.15933995e-16]]

1.08e-15 means $1.08 \times 10^{-15}$ , which is zero up to rounding error.

When the Inverse Fails

Let's break the normal equation on purpose. We add a second feature that is just 2 times the first, so one feature is a combination of the other and $\mathbf{X}^\top \mathbf{X}$ has no inverse:

import numpy as np

np.random.seed(42)  # the same simulated data as above
X = 2 * np.random.rand(100, 1)
y = 5 + 3 * X + np.random.randn(100, 1)

# add a duplicated feature: the second feature is just 2 times the first
X_dup = np.c_[np.ones((100, 1)), X, 2 * X]

# the normal equation: X^T X has no inverse, so this fails
try:
    w_inv = np.linalg.inv(X_dup.T @ X_dup) @ X_dup.T @ y
    print(w_inv)
except np.linalg.LinAlgError as e:
    print('LinAlgError:', e)
print('*' * 50)

# lstsq uses SVD, so it still works
w_svd = np.linalg.lstsq(X_dup, y, rcond=None)[0]
print(w_svd)
print('*' * 50)

# the predictions match the model without the duplicated feature
X_new = np.array([[0.6], [0.9], [1.3]])
print(np.c_[np.ones((3, 1)), X_new, 2 * X_new] @ w_svd)

Output:
LinAlgError: Singular matrix
**************************************************
[[5.21509616]
 [0.55402268]
 [1.10804535]]
**************************************************
[[6.87716419]
 [7.70819821]
 [8.81624356]]

np.linalg.inv fails. Because of rounding, on some computers it does not raise an error and instead prints huge, meaningless numbers, which is worse because nothing warns you.
np.linalg.lstsq solves the same least-squares problem with SVD, the same approach scikit-learn's LinearRegression uses. It finds $b = 5.215$ and splits the weight between the two copies: $w_1 = 0.554$ and $w_2 = 1.108$ . Since the second feature is $2x$ , the model is $5.215 + 0.554x + 1.108 \cdot 2x = 5.215 + 2.770x$ , the same line as before.
The predictions are identical to the lecture's. The data cannot tell how the weight should be split between two copies of the same information, so many $\mathbf{w}$ give the same predictions; SVD picks the one with the smallest weights.

The idea behind SVD is that it breaks $\mathbf{X}$ into simpler pieces, one per independent direction in the data, and skips any direction that carries no information, here the duplicated feature, so it never has to divide by zero. You do not need the details to use it: LinearRegression().fit() does this for you.