Fitting the Model: Gradient Descent and the Normal Equation

1 Learning objectives

2 Motivating example

Tutorial 3 ended with a formula: ŵ = (XᵀX)⁻¹Xᵀy. Type it into NumPy and you get the answer. No iteration, no learning rate, no choices to make. So why would anyone bother with gradient descent?

There are two answers, and the first one is decisive.

The formula does not exist for most models. Logistic regression — the very next unit — has no closed form. Neither do decision trees, random forests, support vector machines or neural networks. Every single model you will meet after Unit 6 is fitted by some iterative procedure. If you cannot walk downhill, you can fit exactly one model, and you have already fitted it.

The formula also becomes impractical at scale. The cost of forming and inverting XᵀX grows as O(d³) in the number of features. With a thousand features that is a billion operations per inversion, and the matrix alone occupies eight megabytes. With a hundred thousand features — perfectly ordinary for text data — it is a trillion operations and eighty gigabytes. The normal equation becomes unusable long before gradient descent does.

So gradient descent is not a fallback for when the algebra is hard. It is the primary tool of modern machine learning, and the normal equation is a convenient special case that happens to be available here.

3 Concept explanation

3.1 The gradient descent algorithm

From Tutorial 4, the gradient of the MSE cost function is:

∂J/∂w = −(2/n) Σ x_i e_i      ∂J/∂b = −(2/n) Σ e_i

where e_i = y_i − (w x_i + b). The gradient points uphill, so the descent direction is its negative. The gradient descent update rule is therefore:

Expanded into the residuals, the whole step is:

eᵢ  = yᵢ − (w xᵢ + b)              for every i

w  ←  w  +  α · (2/n) Σ xᵢ eᵢ
b  ←  b  +  α · (2/n) Σ eᵢ

Read the update for w as a sentence. Each example contributes to the step in proportion to its feature value x_i and its residual e_i. A large x with a positive residual pushes the slope up strongly; a small x barely contributes regardless of its residual. And if the residuals cancel out — as they do at the optimum — the step becomes zero.

The full algorithm, written out:

1. Initialise Set w and b to some starting values. The conventional choice is zero for both, though random small values are also used.
2. Compute predictions and residuals Using the current parameters, calculate ŷ_i and e_i for every training example.
3. Compute the gradient ∂J/∂w and ∂J/∂b from the residuals.
4. Update both parameters Move each one a small step against its partial derivative, using the same learning rate α.
5. Check for convergence If the gradient is small enough, or the cost has stopped improving, or you have hit the iteration limit — stop. Otherwise return to step 2.

3.2 The learning rate

The learning rate α controls the size of each step. It is the single most important hyperparameter in gradient descent, and the one you will spend the most time tuning in practice.

For a quadratic cost function there is a precise bound. If the largest curvature of the surface — the largest eigenvalue of the Hessian — is λ_max, then gradient descent converges if and only if:

0 < α < 2 / λ_max

Below that bound it converges; above it, it diverges. At exactly the bound it oscillates forever without settling. In practice you rarely compute λ_max explicitly, but the bound explains the empirical rule: if your cost curve is diverging, halve the learning rate and try again.

3.3 Convergence and stopping criteria

Gradient descent approaches the minimum asymptotically. In exact arithmetic it never actually arrives — each step gets you closer, but there is always a step remaining. In floating-point arithmetic the cost eventually stops changing meaningfully and further iterations are wasted. So you need a stopping rule.

Three are common, and they are usually combined:

A practical pattern is to set a generous iteration cap, a small cost-change tolerance, and a patience of ten iterations. Run until any of the three triggers. Then plot the cost curve and confirm visually that the plateau is real and not an artefact of the stopping rule.

3.4 Feature scaling

Tutorial 4 established that the shape of the loss surface depends on the scale of the features. This section makes the consequences concrete.

Suppose one feature is measured in metres (range 0–3) and another in millimetres (range 0–3000). The Hessian entries corresponding to the two coefficients differ by a factor of roughly a million, so the contours become extremely elongated. On such a surface:

• The maximum safe learning rate is set by the steepest direction — the large-scale feature — and must therefore be very small.
• Progress along the flat direction — the small-scale feature — is correspondingly slow. You could need tens of thousands of iterations to resolve a coefficient whose natural scale would have taken ten.

Two standard rescaling methods fix this:

Standardisation is the safer default. It is less distorted by outliers than min–max, and the resulting features have a natural interpretation: a value of +2 means "two standard deviations above average".

3.5 The normal equation

From Tutorial 3, the least-squares solution in matrix form is:

This is not an approximation and it does not iterate. It produces the exact minimiser of the MSE cost in a single computation, and the answer is identical to the limit that gradient descent would converge to given infinite iterations.

Computational cost. The two expensive operations are:

• Forming XᵀX, which is a (d+1) × (d+1) matrix. Each of its (d+1)² entries is a dot product over n examples, giving O(n d²).
• Inverting it, which is O(d³).

The total is O(n d² + d³), and the memory required is O(d²) for the matrix itself. That last term is what kills it: with d = 10,000 features, the matrix has 10⁸ entries — eight hundred megabytes in double precision, before you have done any arithmetic.

In practice, scikit-learn does not form the inverse explicitly. It uses a QR decomposition or a singular value decomposition of X, which is more numerically stable — forming (XᵀX)⁻¹ squares the condition number of the problem and can lose precision badly when features are nearly collinear. The mathematical answer is the same; the numerical route is better. This is why the LinearRegression documentation warns about the case where the design matrix is ill-conditioned.

3.6 Choosing between them

For linear regression specifically, both methods produce the same answer, so the choice is entirely about cost and numerical behaviour. The decision rule is straightforward.

The threshold of around a thousand features is a rough guideline, not a law. The real question is whether the O(d³) term is affordable. On a modern laptop, inverting a 1,000 × 1,000 matrix takes a fraction of a second; a 10,000 × 10,000 matrix takes minutes and uses most of your memory. Compute both sides of the trade-off for your problem rather than assuming.

4 Visual illustration

5 Worked example

We will run gradient descent on the bike-share data, twice — once on the raw temperature values and once on standardised values — and compare. The contrast between the two runs is the whole point of this tutorial.

Σ eᵢ  = 120 + 180 + 240 + 300 + 360 + 430 + 510 + 560  =  2,700

Σ xᵢeᵢ = 8(120) + 12(180) + 15(240) + 18(300)
        + 21(360) + 24(430) + 27(510) + 30(560)
      = 960 + 2,160 + 3,600 + 5,400
        + 7,560 + 10,320 + 13,770 + 16,800
      = 60,570

∂J/∂w = −(2/8)(60,570)  = −15,142.50
∂J/∂b = −(2/8)( 2,700)  =    −675.00

Gradient  =  ( −15,142.50 , −675.00 )

H  =  (2/n) XᵀX

         ⎡ 850.75   38.75 ⎤
   =     ⎢                ⎥
         ⎣  38.75    2.00 ⎦

Eigenvalues:   λ_max ≈ 852.57 ,  λ_min ≈ 0.18

Maximum safe learning rate:  2 / λ_max ≈ 0.00235
Condition number:            852.57 / 0.18 ≈ 4,845

w₁ = 0 − 0.002 × (−15,142.50)  =  30.285
b₁ = 0 − 0.002 × (−675.00)     =   1.350

x̄  = 155 / 8          = 19.375
s   = √(Sxx / n)
    = √(399.875 / 8)
    = √49.984
    ≈ 7.070

H  =  ⎡ 2   0 ⎤        λ_max = λ_min = 2
      ⎣ 0   2 ⎦

Maximum safe learning rate:  2 / 2 = 1.0
Condition number:            2 / 2 = 1

The same computation in code

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler

df = pd.DataFrame({
    "temp_c":  [ 8, 12, 15, 18, 21, 24, 27, 30],
    "rentals": [120, 180, 240, 300, 360, 430, 510, 560],
})
x = df["temp_c"].to_numpy()
y = df["rentals"].to_numpy()
n = len(x)


def mse(w, b, xv=x, yv=y):
    return ((yv - (w * xv + b)) ** 2).mean()


def gradient(w, b, xv=x, yv=y):
    e = yv - (w * xv + b)
    return -2 * (xv * e).mean(), -2 * e.mean()


def descend(xv, w, b, alpha, iters=50000, tol=1e-8):
    history = [(w, b, mse(w, b, xv, y))]
    for _ in range(iters):
        gw, gb = gradient(w, b, xv, y)
        w_new, b_new = w - alpha * gw, b - alpha * gb
        if abs(w_new - w) < tol and abs(b_new - b) < tol:
            w, b = w_new, b_new
            history.append((w, b, mse(w, b, xv, y)))
            break
        w, b = w_new, b_new
        history.append((w, b, mse(w, b, xv, y)))
    return history


# --- Run 1: raw temperature ---
raw = descend(x, 0.0, 0.0, alpha=0.001)
print(f"raw   : {len(raw):>6} iterations, "
      f"w = {raw[-1][0]:.4f}, b = {raw[-1][1]:.4f}, "
      f"MSE = {raw[-1][2]:.4f}")

# --- Run 2: standardised temperature ---
scaler = StandardScaler()
x_std = scaler.fit_transform(x.reshape(-1, 1)).ravel()

std = descend(x_std, 0.0, 0.0, alpha=0.4)
print(f"std   : {len(std):>6} iterations, "
      f"w̃ = {std[-1][0]:.4f}, b̃ = {std[-1][1]:.4f}, "
      f"MSE = {std[-1][2]:.4f}")

# --- Compare against the normal equation ---
X = np.column_stack([np.ones(n), x])
w_ne = np.linalg.solve(X.T @ X, X.T @ y)
print(f"normal: b = {w_ne[0]:.4f}, w = {w_ne[1]:.4f}, "
      f"MSE = {mse(w_ne[1], w_ne[0]):.4f}")

# --- And against scikit-learn ---
sk = LinearRegression().fit(x.reshape(-1, 1), y)
print(f"sklearn: b = {sk.intercept_:.4f}, w = {sk.coef_[0]:.4f}")

raw   :  38515 iterations, w = 20.6473, b = -62.5598, MSE = 103.8703
std   :      6 iterations, w̃ = 146.1260, b̃ = 337.4986, MSE = 103.8700
normal: b = -62.5976, w = 20.6502, MSE = 103.8703
sklearn: b = -62.5976, w = 20.6502

All four routes — gradient descent on raw features, gradient descent on scaled features, the normal equation and scikit-learn — arrive at the same MSE to four decimal places. The scaled run needed six iterations; the raw run needed thirty-eight and a half thousand. Same answer, same algorithm, same data.

6 Interactive activity

7 Self-check questions

Answer each question in your own words before revealing the answer.