Gradient Descent in Depth — and Building It in NumPy

What this lesson delivers

• The gradient of the MSE loss, derived from scratch and then written in one vectorized line.
• A from-scratch batch gradient-descent loop in NumPy that converges to the Lesson 10 answer.
• A working feel for the learning rate: too small crawls, too large explodes, and there's a sweet spot you can find by reading a loss curve.
• Why feature scaling is not optional for gradient descent — the same fit that needs 40,000 steps on raw data lands in one step when standardized.
• Batch vs. mini-batch vs. stochastic gradient descent, each as a small edit to the same loop.
• How to decide when the optimizer has converged and it's time to stop.

1 · Why walk down when you could jump?

Recall the closing tension of Lesson 10. The normal equation θ = (X^TX)⁻¹X^Ty hands you the exact best weights in one shot — but forming and inverting X^TX costs on the order of d³ arithmetic for d features, and needs the whole dataset in memory at once. With a dozen features that's nothing; with fifty thousand features, or ten million rows that don't fit in RAM, it is impossible. And for models whose loss is not a neat convex bowl, there is no formula to invert in the first place.

Gradient descent trades one expensive exact answer for many cheap approximate steps. Each step looks at the slope of the loss and nudges the parameters a little way downhill; enough steps and you arrive at (essentially) the same minimum, having never inverted a matrix. Its per-step cost is only the cost of evaluating the gradient — linear in the data, not cubic in the features — and, as we'll see with the stochastic variant, you don't even need all the data in memory to take a step.

2 · The gradient of the loss, derived by hand

Gradient descent needs one thing the concept lesson took on faith: the actual gradient of this loss. Let's derive it — it is three lines of the same calculus from Lesson 10's closed-form section, stopped one step earlier. We never set the derivative to zero this time; we keep it, because the derivative is the downhill direction.

Write the residual of example i as r⁽ⁱ⁾ = ŷ⁽ⁱ⁾ − y⁽ⁱ⁾ = (w x⁽ⁱ⁾ + b) − y⁽ⁱ⁾. The loss is the mean of the squared residuals:

Differentiate with respect to each parameter. The square's exponent drops to a factor of 2 (chain rule), and the inside's derivative is x⁽ⁱ⁾ for w and 1 for b — exactly the two derivatives you computed in Lesson 10, only now we stop before setting them to zero:

The vectorized gradient — one line for any number of features

Writing two separate sums doesn't scale to fifty features. As in Lesson 10, fold the bias into the weights: add a column of 1s to the data so the parameter vector is θ = (b, w₁, …) and the predictions for all n rows are the single matrix–vector product Xθ. Then the whole residual vector is r = Xθ − y, and both partials collapse into one gradient vector:

That single expression is the entire mathematical content of gradient descent for Linear Regression. Everything else is a loop around it. Notice it reuses the exact quantities from the normal equation — X^TX and X^Ty live inside it — which is why setting ∇J = 0 gives back that formula. The two methods are the same landscape; one jumps to the bottom, the other walks.

3 · The update rule, and the algorithm in full

The update is Lesson 7's rule, now with a gradient we can actually compute. Move the parameters a small step in the direction opposite the gradient (downhill is the negative of "the way the loss grows"), scaled by the learning rate α:

The full algorithm is nothing more than repeating that line:

4 · From scratch in NumPy

We'll fit the same six-house dataset from Lesson 10, so we can check the optimizer against the closed-form answer w ≈ 2.104, b ≈ 47.87. One preparatory step needs a word of justification now and a full one in section 6: we first standardize the feature — subtract its mean and divide by its standard deviation. Skip this and, as we'll see, raw metres force a learning rate near 0.0001 and tens of thousands of steps. Standardized, the same fit is almost instantaneous. Trust the move for now; section 6 shows exactly why.

First the transparent, two-scalar version — a direct transcription of the two partials from section 2, so you can see each gradient by name:

import numpy as np

# Same toy houses as Lesson 10: size (m2) -> price (thousands $)
x = np.array([ 50,  60,  80,  95, 120, 150], dtype=float)
y = np.array([150, 175, 210, 260, 300, 360], dtype=float)

# Standardize the feature: z has mean 0, std 1  (see section 6 for why)
x_mean, x_std = x.mean(), x.std()
z = (x - x_mean) / x_std
n = len(y)

w, b = 0.0, 0.0       # start at the origin
alpha = 0.1          # learning rate

for step in range(500):
    y_hat = w * z + b                 # predictions
    resid = y_hat - y                 # r = y_hat - y
    grad_w = (2 / n) * np.sum(resid * z)   # dJ/dw  (section 2)
    grad_b = (2 / n) * np.sum(resid)       # dJ/db
    w -= alpha * grad_w               # step downhill
    b -= alpha * grad_b

print(w, b)      # ~ 72.3   242.5   (parameters in STANDARDIZED space)

Those aren't the Lesson 10 numbers — they're the parameters of the line in standardized coordinates. That's expected: we changed the feature's units, so the slope and intercept changed with them. Converting back to real metres and dollars is pure algebra (a slope divided by the feature's spread, an intercept shifted to undo the centring):

w_real = w / x_std                       # undo the /std on the slope
b_real = b - w * x_mean / x_std          # undo the centring
print(w_real, b_real)   # ~ 2.104   47.87   <- exactly the Lesson 10 closed form

The vectorized version — the one you'll actually reuse

The two-scalar loop only handles a single feature. Here is the general form, straight from ∇J = (2/n) X^Tr. Swap in a wider X and it fits any number of features without another line changing — and it records the loss each step so we can plot convergence:

# Design matrix: bias column of 1s beside the standardized feature
X = np.column_stack([np.ones_like(z), z])   # shape (6, 2): [1, z] per row

theta = np.zeros(X.shape[1])    # [b, w], generalizes to [b, w1, w2, ...]
alpha = 0.1
history = []

for step in range(500):
    resid = X @ theta - y                   # r = X.theta - y   (all rows at once)
    grad  = (2 / n) * (X.T @ resid)         # the whole gradient in one line
    theta -= alpha * grad
    history.append(np.mean(resid ** 2))    # track the loss to watch it fall

b, w = theta
print(w / x_std, b - w * x_mean / x_std)   # ~ 2.104   47.87  (same answer)

5 · The learning rate: the one knob that decides everything

Lesson 7 warned that the step size α is the single most consequential number in the procedure. Now you can watch it happen. On the standardized house data, run the loop above at several learning rates and look at the loss curve:

6 · Feature scaling: why the bowl's shape is destiny

Now the debt from section 4 comes due. Why standardize at all — the closed form in Lesson 10 didn't need it? Because gradient descent, unlike the formula, is exquisitely sensitive to the shape of the loss bowl, and the shape is set by the scale of the features.

Feed the loop the raw feature (sizes up to 150) instead of the standardized one, and it falls apart. The gradient's weight-component contains x² terms as large as 150² = 22{,}500, while the bias-component is of order 1. The two parameters therefore live on wildly different slopes: the loss surface is not a round bowl but a long, narrow valley — for this data, about 80,000 times longer than it is wide. That is precisely the failure mode Lesson 7 promised would "actually bite you."

# The SAME loop on RAW x instead of standardized z -- watch it struggle
X_raw = np.column_stack([np.ones_like(x), x])   # [1, size], size up to 150

# alpha = 0.1 here -> loss overflows to inf almost immediately.
# To even stay stable you must shrink the step to about:
alpha = 0.0001
# ...and then it needs on the order of 40,000 steps to reach
# the answer that standardized data reached in ~100.

7 · Batch, mini-batch, and stochastic — how much data per step

Every loop so far computed the gradient over all n examples before taking one step. That's batch gradient descent: the most accurate downhill direction, but each step pays for a full pass through the data. When the data is millions of rows, one step becomes painfully expensive — and you'd like to start improving before you've even read everything once. The fix is to estimate the gradient from a subset of the data each step. This is the same update rule; only the rows feeding the gradient change.

Here is stochastic gradient descent as a small edit to the same loop — shuffle the data, then step once per example. An epoch is one full sweep through the dataset:

theta = np.zeros(X.shape[1])
alpha = 0.01          # noisier steps -> usually a smaller rate than batch

for epoch in range(50):                  # 50 sweeps through the data
    order = np.random.permutation(n)      # reshuffle each epoch
    for i in order:
        xi = X[i:i+1]                    # one row, kept 2-D: shape (1, features)
        yi = y[i:i+1]
        resid = xi @ theta - yi           # gradient from a SINGLE example
        grad  = 2 * (xi.T @ resid)          # (2/1) * x_i^T r_i
        theta -= alpha * grad

For mini-batch, change one thing: instead of one index i, slice the shuffled order into chunks and feed each chunk (a small X) through the very same vectorized gradient (2/m) X^Tr, where m is the chunk size. The gradient line you already wrote handles it unchanged — that's the payoff of vectorizing.

Knowing when to stop

The closed form finishes when the algebra ends; an iterative method has to decide it's done. Three common criteria, usually combined:

• The loss stops improving. When the change in loss between steps falls below a small tolerance, the marble has effectively reached the floor — further steps only shuffle floating-point dust.
• The gradient goes (nearly) to zero. A tiny gradient magnitude means the ground is flat — the definition of a minimum, and the same condition the normal equation solved for directly.
• A step budget. A hard cap on iterations (or epochs) so a badly-tuned run can't loop forever — and, later, a way to stop before a model starts overfitting (Lesson 13).

8 · What to keep honest about

How this connects to the next lesson

Exercises & reflection

The first is with pencil; the rest are at the keyboard, extending the code above. If you can do them, gradient descent is genuinely yours.