Linear & Logistic Regression

Gradient Descent from Scratch in Python

Suresh Madhra·Jun 24, 2026·12 min read

One idea, three variants

Almost every model you will train is fitted by the same procedure: measure the error, work out which direction increases it, and step the other way. That is gradient descent. The three variants people argue about — batch, stochastic and mini-batch — differ only in how much data they look at before taking each step.

θ_{t + 1} = θ_{t} - η \nabla_{θ} J (θ_{t})

θ are the parameters, J is the loss, ∇J is the gradient, and η is the learning rate — the single most important number in the whole procedure.

1. The loss and its gradient

For linear regression with mean squared error the gradient is easy enough to derive by hand, which makes it the perfect place to start.

J (b) = \frac{1}{n} i = 1 \sum n (b_{0} + b_{1} x_{i} - y_{i})^{2}

\frac{\partial J}{\partial b _{0}} = \frac{2}{n} i \sum e_{i}, \frac{\partial J}{\partial b _{1}} = \frac{2}{n} i \sum e_{i} x_{i}

2. Batch gradient descent

Compute the gradient over the whole dataset, then take one step. Smooth, deterministic, and slow — every step costs a full pass over the data.

batch_gd.py

data = [(1, 3.1), (2, 4.9), (3, 7.2), (4, 8.8), (5, 11.1)]   # y ~ 2x + 1

def loss(b0, b1):
    return sum((b0 + b1 * x - y) ** 2 for x, y in data) / len(data)

def batch_gradient_descent(lr=0.02, epochs=200):
    b0 = b1 = 0.0
    n = len(data)
    for epoch in range(epochs):
        g0 = sum((b0 + b1 * x - y) for x, y in data) * 2 / n
        g1 = sum((b0 + b1 * x - y) * x for x, y in data) * 2 / n
        b0 -= lr * g0
        b1 -= lr * g1
        if epoch % 40 == 0:
            print(">epoch=%3d loss=%.5f b0=%.4f b1=%.4f" % (epoch, loss(b0, b1), b0, b1))
    return b0, b1

b0, b1 = batch_gradient_descent()
print()
print("Fitted: y = %.3f + %.3f x" % (b0, b1))

3. Stochastic gradient descent

Update after every single example. The path to the minimum is noisy, but each step is cheap and the noise itself helps escape shallow local minima in non-convex problems.

sgd.py

from random import seed, shuffle

data = [(1, 3.1), (2, 4.9), (3, 7.2), (4, 8.8), (5, 11.1)]

def loss(b0, b1):
    return sum((b0 + b1 * x - y) ** 2 for x, y in data) / len(data)

def sgd(lr=0.02, epochs=40):
    seed(1)
    b0 = b1 = 0.0
    rows = list(data)
    for epoch in range(epochs):
        shuffle(rows)
        for x, y in rows:
            error = b0 + b1 * x - y
            b0 -= lr * 2 * error
            b1 -= lr * 2 * error * x
        if epoch % 8 == 0:
            print(">epoch=%3d loss=%.5f b0=%.4f b1=%.4f" % (epoch, loss(b0, b1), b0, b1))
    return b0, b1

b0, b1 = sgd()
print()
print("Fitted: y = %.3f + %.3f x  (in a fraction of the passes batch GD needed)" % (b0, b1))

4. Mini-batch, momentum and decay

Mini-batch is the practical compromise: average the gradient over 32–256 examples. Momentum accumulates a velocity so consistent directions accelerate and oscillations cancel. A decaying learning rate lets you move fast early and settle precisely later.

v_{t + 1} = β v_{t} + (1 - β) \nabla J, θ_{t + 1} = θ_{t} - η v_{t + 1}

minibatch_momentum.py

from random import seed, shuffle

data = [(x, 2.0 * x + 1.0 + (0.1 if x % 2 else -0.1)) for x in range(1, 41)]

def loss(b0, b1):
    return sum((b0 + b1 * x - y) ** 2 for x, y in data) / len(data)

def minibatch(lr=0.01, epochs=60, batch_size=8, beta=0.9, decay=0.02):
    seed(1)
    b0 = b1 = 0.0
    v0 = v1 = 0.0
    rows = list(data)
    for epoch in range(epochs):
        shuffle(rows)
        step = lr / (1.0 + decay * epoch)
        for i in range(0, len(rows), batch_size):
            batch = rows[i:i + batch_size]
            m = len(batch)
            g0 = sum((b0 + b1 * x - y) for x, y in batch) * 2 / m
            g1 = sum((b0 + b1 * x - y) * x for x, y in batch) * 2 / m
            v0 = beta * v0 + (1 - beta) * g0
            v1 = beta * v1 + (1 - beta) * g1
            b0 -= step * v0
            b1 -= step * v1
        if epoch % 12 == 0:
            print(">epoch=%3d lr=%.4f loss=%.5f b0=%.4f b1=%.4f" % (epoch, step, loss(b0, b1), b0, b1))
    return b0, b1

b0, b1 = minibatch()
print()
print("Fitted: y = %.3f + %.3f x" % (b0, b1))

from random import seed, shuffle

data = [(x, 2.0 * x + 1.0 + (0.1 if x % 2 else -0.1)) for x in range(1, 41)]

def loss(b0, b1):
    return sum((b0 + b1 * x - y) ** 2 for x, y in data) / len(data)

def minibatch(lr=0.01, epochs=60, batch_size=8, beta=0.9, decay=0.02):
    seed(1)
    b0 = b1 = 0.0
    v0 = v1 = 0.0
    rows = list(data)
    for epoch in range(epochs):
        shuffle(rows)
        step = lr / (1.0 + decay * epoch)
        for i in range(0, len(rows), batch_size):
            batch = rows[i:i + batch_size]
            m = len(batch)
            g0 = sum((b0 + b1 * x - y) for x, y in batch) * 2 / m
            g1 = sum((b0 + b1 * x - y) * x for x, y in batch) * 2 / m
            v0 = beta * v0 + (1 - beta) * g0
            v1 = beta * v1 + (1 - beta) * g1
            b0 -= step * v0
            b1 -= step * v1
        if epoch % 12 == 0:
            print(">epoch=%3d lr=%.4f loss=%.5f b0=%.4f b1=%.4f" % (epoch, step, loss(b0, b1), b0, b1))
    return b0, b1

b0, b1 = minibatch()
print()
print("Fitted: y = %.3f + %.3f x" % (b0, b1))

5. Diagnosing a bad run

Loss increases or becomes nan: learning rate too high, or features unscaled. Divide η by 10 and standardize.
Loss falls then plateaus far above zero: model too simple, or learning rate now too large for fine-tuning — decay it.
Loss falls agonisingly slowly: learning rate too small, or features on wildly different scales.
Loss oscillates every epoch with SGD: normal. Watch a moving average, not the raw value.
Train loss falls, validation loss rises: you are past the sweet spot. Stop early.

Plot the loss curve every single time. It is the cheapest diagnostic in machine learning and it catches most bugs before the metrics do.

← Back to blog