Gradient descent with many features

You already know how to train a one-feature model: measure the cost, find the slope, step downhill, repeat. With $n$ features there are $n + 1$ parameters ($w_1$ to $w_n$ and $b$), so the cost landscape has $n + 1$ directions instead of two. That sounds harder to picture, and it is, but the algorithm hardly changes. Each parameter simply gets its own slope.

The cost is the same idea

The cost still measures the average squared error over the training set. The only change is that the prediction now uses the dot product:

With vectorisation that’s three short lines: predict everything, take the errors, average their squares.

import numpy as np

X = np.array([
    [1850, 3, 2, 12],
    [1200, 2, 1, 40],
    [2400, 4, 2, 8],
    [950, 2, 1, 55],
    [1600, 3, 1, 25],
    [3100, 5, 2, 3],
], dtype=float)
y = np.array([412, 236, 538, 178, 330, 690], dtype=float)
m, n = X.shape

def compute_cost(X, y, w, b):
    err = X @ w + b - y              # one error per house
    return (err ** 2).sum() / (2 * len(y))

print("cost with all zeros:", compute_cost(X, y, np.zeros(n), 0.0))
print("cost with a guess:  ", compute_cost(X, y, np.array([0.2, 10, 5, -1]), 30.0))

One slope per parameter

In Math for ML, section 3, you met partial derivatives: the slope of a function along one input while the others stay put, and the gradient, the vector of all of them. For our cost, the partial derivative with respect to weight $w_j$ is:

Compare it with Section 1: it’s the same formula, with $x^{(i)}$ replaced by feature $j$ of example $i$. Each weight’s slope is “the errors, weighted by how much of that feature each house had”. If big houses are under-priced, the error is negative where size is large, the slope for the size weight is negative, and gradient descent pushes that weight up.

Then every parameter steps downhill at the same time, a simultaneous update exactly as before:

The whole gradient in one line

Writing a loop over $j$ would work, but look at what the sum does: for each feature, multiply the error vector by that feature’s column and add up. That’s a dot product of the error vector with each column of $X$, and doing it for all columns at once is a matrix product with the transpose $X^T$:

Shapes check: $X^T$ is $(n, m)$, the error vector is $(m,)$, so the result is $(n,)$, one slope per weight. In NumPy:

def compute_gradient(X, y, w, b):
    err = X @ w + b - y              # shape (m,)
    dj_dw = X.T @ err / len(y)       # shape (n,): one slope per weight
    dj_db = err.mean()               # a single number
    return dj_dw, dj_db

dj_dw, dj_db = compute_gradient(X, y, np.zeros(n), 0.0)
print("slopes for w:", dj_dw)
print("slope for b: ", dj_db)

All the slopes are negative at $\vec{w} = 0$: every prediction is zero, far below every price, so increasing any weight would help. Notice how different their sizes are, though. The size slope is in the hundreds of thousands, the floors slope is in the hundreds. Keep that in mind.

Train it

Put the pieces in a loop. We record the cost after each step so we can check that it falls.

def gradient_descent(X, y, w, b, alpha, iters):
    history = []
    for _ in range(iters):
        dj_dw, dj_db = compute_gradient(X, y, w, b)
        w = w - alpha * dj_dw        # all weights at once
        b = b - alpha * dj_db
        history.append(compute_cost(X, y, w, b))
    return w, b, history

w, b, history = gradient_descent(X, y, np.zeros(n), 0.0, alpha=1e-7, iters=1000)
print("cost: start", round(history[0], 1), "-> end", round(history[-1], 1))
print("weights:", np.round(w, 4), " b:", round(b, 4))
print("predictions:", np.round(X @ w + b), "\nactual:     ", y)

It learned something: the cost fell and the predictions are in the right ballpark. But look at the weights. The size weight moved to about 0.22, while the bedrooms and floors weights are still essentially zero after a thousand steps, and $b$ has barely moved. Gradient descent has learnt “price is about 0.22 times size” and almost nothing else.

Why such a tiny learning rate? Try it: change alpha to 1e-6 and run again. The cost explodes. Size is measured in thousands, so its slope is huge, and any bigger step overshoots wildly along the size direction. But that same tiny step is far too small for features measured in single digits. One learning rate can’t suit features on such different scales. The next lesson fixes this, and the fix is pleasingly simple.

Your turn

import numpy as np

X = np.array([
    [1850, 3, 2, 12],
    [1200, 2, 1, 40],
    [2400, 4, 2, 8],
    [950, 2, 1, 55],
    [1600, 3, 1, 25],
    [3100, 5, 2, 3],
], dtype=float)
y = np.array([412, 236, 538, 178, 330, 690], dtype=float)

def compute_cost(X, y, w, b):
    m = len(y)
    # your code here
    return None

def compute_gradient(X, y, w, b):
    m = len(y)
    # your code here
    dj_dw = None
    dj_db = None
    return dj_dw, dj_db

print(compute_cost(X, y, np.zeros(4), 0.0))
print(compute_gradient(X, y, np.zeros(4), 0.0))

Check yourself