Machine Learning from scratch: Bare bones implementations in Python

Machine Learning From Scratch: Bare Bones Implementations in Python

Machine learning from scratch means writing the training loop yourself — the parameter update, the loss function, and the backward pass — instead of calling model.fit(). The point is not to replace scikit-learn. The point is that once you have written the update rule by hand, every framework stops being magic and starts being convenience.

Here is the problem this tutorial solves: you understand what a gradient is, you can read the math, but you have never watched a loss number actually go down because of a line of code you wrote. So below are three complete, runnable implementations — a raw gradient descent loop, a single sigmoid neuron with cross-entropy, and a two-layer network with a hand-written backward pass. Every one is under 40 lines. Copy them, run them, break them.

Setup: NumPy only, and why

You need NumPy and nothing else. No PyTorch, no TensorFlow, no scikit-learn. That constraint is the whole exercise. If you import a framework you will be tempted to use its autograd, and autograd is exactly the thing you are trying to understand.

import numpy as np
np.random.seed(0)

Seeding matters more than people expect. Without it, two runs of the same script give different losses and you cannot tell whether your change to the learning rate helped or whether you just got a luckier initialization.

Implementation 1: Gradient descent on a line

Start with the smallest possible learning problem: fit y = w*x + b to data you generate yourself. You know the true answer, so you can verify the code found it. That is the entire value of this first step — a known ground truth.

# True relationship: y = 3x + 2, plus noise
x = np.linspace(-1, 1, 100)
y = 3.0 * x + 2.0 + np.random.normal(0, 0.1, size=x.shape)

w, b = 0.0, 0.0
lr = 0.1

for step in range(200):
    y_hat = w * x + b
    error = y_hat - y
    loss = np.mean(error ** 2)

    # d(loss)/dw and d(loss)/db, derived by hand
    grad_w = np.mean(2 * error * x)
    grad_b = np.mean(2 * error)

    w -= lr * grad_w
    b -= lr * grad_b

    if step % 50 == 0:
        print(f"step {step:3d}  loss {loss:.4f}  w {w:.3f}  b {b:.3f}")

Run it and the printed w and b walk toward 3.0 and 2.0. The loss falls. That is gradient descent, and there is nothing else in it.

Two things break this if you are careless. First, the learning rate: at lr = 0.1 this converges; push it well past 1.0 with unnormalized inputs and the loss will oscillate and then explode into nan. Second, the 2 * in the gradient is not optional — it is the derivative of the squared error. Drop it and you are silently training at half your intended learning rate. That kind of bug is invisible in the output, which is why the next section exists.

The habit that catches silent gradient bugs: numerical checking

A gradient check compares your hand-derived gradient against a finite-difference estimate. Nudge one parameter by a tiny epsilon, measure how much the loss moves, divide. If your analytic gradient is right, the two numbers agree to several decimal places.

def numerical_grad(f, params, idx, eps=1e-5):
    p = params.copy()
    p[idx] += eps
    loss_plus = f(p)
    p[idx] -= 2 * eps
    loss_minus = f(p)
    return (loss_plus - loss_minus) / (2 * eps)

def loss_fn(params):
    w_, b_ = params
    return np.mean((w_ * x + b_ - y) ** 2)

params = np.array([w, b])
analytic = np.array([grad_w, grad_b])
numeric = np.array([numerical_grad(loss_fn, params, i) for i in range(2)])

print("analytic:", analytic)
print("numeric: ", numeric)
print("relative diff:", np.abs(analytic - numeric) / (np.abs(analytic) + np.abs(numeric) + 1e-8))

The central difference — nudging both up and down and dividing by 2 * eps — is worth the extra forward pass. A one-sided difference is only accurate to about eps; the central version is accurate to about eps squared. With eps = 1e-5 that is the difference between three and five trustworthy digits. A relative difference under roughly 1e-6 means your gradient is correct. Anything above 1e-3 means you have a real bug, not floating-point noise.

Do this once per new loss function, not on every training run. It is a debugging instrument, and it costs a full forward pass per parameter, which is why nobody runs it inside the loop.

Implementation 2: One neuron, sigmoid, cross-entropy

Now a classifier. One neuron, a sigmoid squashing the output into a probability, cross-entropy as the loss. The reason this pairing is the standard teaching example — and the reason it is worth writing by hand once — is what happens when you differentiate it.

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

# Two clusters: class 0 near (-1,-1), class 1 near (1,1)
X = np.vstack([np.random.normal(-1, 0.5, (50, 2)),
               np.random.normal( 1, 0.5, (50, 2))])
Y = np.concatenate([np.zeros(50), np.ones(50)])

w = np.zeros(2)
b = 0.0
lr = 0.5

for step in range(500):
    z = X @ w + b
    p = sigmoid(z)
    loss = -np.mean(Y * np.log(p + 1e-12) + (1 - Y) * np.log(1 - p + 1e-12))

    # The cancellation: dL/dz collapses to (p - Y)
    dz = (p - Y) / len(Y)
    grad_w = X.T @ dz
    grad_b = np.sum(dz)

    w -= lr * grad_w
    b -= lr * grad_b

    if step % 100 == 0:
        acc = np.mean((p > 0.5) == Y)
        print(f"step {step:3d}  loss {loss:.4f}  acc {acc:.2f}")

Look at the line dz = (p - Y) / len(Y). The sigmoid derivative p*(1-p) and the cross-entropy derivative -Y/p + (1-Y)/(1-p) cancel algebraically, leaving just the difference between prediction and label. That is not a coincidence and it is not a shortcut — it is the reason sigmoid and cross-entropy are paired everywhere. Pair sigmoid with mean squared error instead and your gradient carries a p*(1-p) factor that shrinks toward zero whenever the neuron is confidently wrong, which is precisely when you most want a large gradient. The network stalls. If you have ever watched a classifier's loss flatten while accuracy stays near chance, that is usually the cause.

The 1e-12 inside the log is not decoration. log(0) is negative infinity, and a single saturated prediction will turn your loss into nan and poison every weight on the next step.

Implementation 3: Two layers, backward pass by hand

This is the one that makes backprop click. Two layers, a ReLU hidden layer, a sigmoid output, and manual chain rule. The XOR problem is the right test because a single neuron provably cannot solve it — if your two-layer net learns it, the hidden layer is genuinely doing work.

X = np.array([[0,0],[0,1],[1,0],[1,1]], dtype=float)
Y = np.array([0, 1, 1, 0], dtype=float).reshape(-1, 1)

W1 = np.random.randn(2, 4) * 0.5
b1 = np.zeros((1, 4))
W2 = np.random.randn(4, 1) * 0.5
b2 = np.zeros((1, 1))
lr = 0.5

for step in range(5000):
    # Forward
    Z1 = X @ W1 + b1
    A1 = np.maximum(0, Z1)          # ReLU
    Z2 = A1 @ W2 + b2
    A2 = sigmoid(Z2)

    loss = -np.mean(Y * np.log(A2 + 1e-12) + (1 - Y) * np.log(1 - A2 + 1e-12))

    # Backward
    dZ2 = (A2 - Y) / len(X)         # same cancellation as before
    dW2 = A1.T @ dZ2
    db2 = np.sum(dZ2, axis=0, keepdims=True)

    dA1 = dZ2 @ W2.T
    dZ1 = dA1 * (Z1 > 0)            # ReLU derivative: 1 where z>0, else 0
    dW1 = X.T @ dZ1
    db1 = np.sum(dZ1, axis=0, keepdims=True)

    W2 -= lr * dW2; b2 -= lr * db2
    W1 -= lr * dW1; b1 -= lr * db1

    if step % 1000 == 0:
        print(f"step {step:5d}  loss {loss:.4f}")

print("predictions:", np.round(A2.ravel(), 3))

The (Z1 > 0) term is the ReLU derivative, and it is the one people get wrong. ReLU is a kink, not a curve — its slope is exactly 1 for positive inputs and exactly 0 for negative ones. Multiply the incoming gradient by that mask and you are done. Forget it and the hidden layer's gradients become meaningless.

There is a real failure mode hiding here. If a hidden unit's weights drift so that its output is negative for every training example, its gradient is zero for every example, forever. The unit is dead and never recovers. With four hidden units on a four-row dataset you will see this occasionally depending on the seed. Leaky ReLU — np.maximum(0.01 * Z1, Z1) — leaves a small slope on the negative side and prevents it.

When to stop hand-rolling

Writing gradients by hand is a learning tool with a short useful life. The decision rule is concrete: hand-roll when the total number of parameters is small enough that you can write every gradient expression on one screen, and when the architecture is a straight feed-forward stack. That covers roughly the first three things you will ever want to try.

Switch to a framework the moment any of these is true: you have a branch in the computation (an if inside the forward pass), you need convolutions or attention, or you have more than about three layers so the chain rule becomes a transcription exercise with no insight left in it. At that point the hand-written backward pass stops teaching you anything and starts costing you hours per bug.

The reason to care about the mechanism rather than just the API is that the framework will not tell you why the loss is not moving. Knowing that dz should collapse to p - Y, that a dead ReLU has a permanently zero gradient, and that a saturated sigmoid kills the gradient through the p*(1-p) factor is what turns "my model doesn't train" from a mystery into a checklist. That knowledge is portable across every framework you will ever use, which is more than can be said for any single API.

One honest limit: these implementations are teaching artifacts. They process the whole dataset in one batch, use a fixed learning rate, and have no regularization, no validation split, and no numerical safeguards beyond the epsilon inside the log. On real data with real scale they will be slow and fragile. That is fine — they were never meant to ship.

Sources

• AI Tool Database (internally verified snapshot), 2026. Internal database of 360 AI tools with pricing and capability snapshots, most recently verified 2026-09-18.