A Network from Scratch: no framework, in this tab

Every code block on this page runs where you are standing. No install, no GPU, no account. By the end of it you will have trained a two-layer neural network to 97% test accuracy on 1,797 hand-written digits, in about a second, using nothing but NumPy and about thirty lines of maths you wrote yourself.

That is not a toy claim to make the lesson sound exciting. Jump to Train It, press Run, and time it. The first Run on the page downloads NumPy and scikit-learn (the dataset lives inside it), which takes a few seconds once; after that, the training itself takes about a second.

There is a reason to do this the hard way once. model.fit() is four characters of typing and teaches you nothing about what went wrong when the loss goes to nan at 2am. Write the backward pass by hand a single time and every framework afterwards stops being magic — you will know exactly which line Keras is calling on your behalf.

1. Real Data, No Download

scikit-learn ships a digits dataset inside the package. 1,797 images, 8×8 pixels, greyscale 0–16, already labelled. No network request, which is why it works in a browser tab.

It is not MNIST. MNIST is 28×28 and 70,000 images; this is its smaller, older cousin from the UCI archive. For learning the mechanics that is an advantage — 64 inputs instead of 784 means the weight matrices are small enough to print and stare at.

from sklearn.datasets import load_digits
import numpy as np

X, y = load_digits(return_X_y=True)
print(X.shape, y.shape)          # (1797, 64) (1797,)
print("pixel range:", X.min(), "to", X.max())
print("first label:", y[0])

# Each row is a flattened 8x8 image. Fold one back up to look at it.
print(X[0].reshape(8, 8).astype(int))

Run that. The 8×8 grid of numbers is a zero — squint and you can see the hole in the middle where the values drop back to 0.

Two things to do before any network sees this data. Scale the pixels into 0–1, because a weight update sized for inputs of 16 will be sixteen times too large for inputs of 1. And turn each label into a one-hot row — the digit 3 becomes [0,0,0,1,0,0,0,0,0,0] — because the network is going to output ten numbers, not one.

from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split
import numpy as np

X, y = load_digits(return_X_y=True)
X = X / 16.0                       # 0..16 -> 0..1
Y = np.eye(10)[y]                  # one-hot: shape (1797, 10)

Xtr, Xte, Ytr, Yte, ytr, yte = train_test_split(
    X, Y, y, test_size=0.25, random_state=0)

print("train", Xtr.shape, "test", Xte.shape)
print("label 3 one-hot:", np.eye(10)[3])

np.eye(10)[y] is the whole one-hot encoder. np.eye(10) is a 10×10 identity matrix; indexing it with an array of labels picks one row per label. No loop, no OneHotEncoder import.

2. The Forward Pass

Two layers. 64 inputs → 64 hidden units → 10 outputs. In matrix form the entire forward pass is four lines:

z1 = X @ W1 + b1        # linear
a1 = relu(z1)           # non-linear
z2 = a1 @ W2 + b2       # linear
p  = softmax(z2)        # probabilities

@ is matrix multiply. X is (batch, 64), W1 is (64, 64), so z1 is (batch, 64) — one row of hidden activations per image. The shapes are the thing to hold onto; almost every bug you will hit here is a shape bug, and NumPy will tell you about it loudly.

Two details that matter more than they look:

ReLU is np.maximum(0, z). That is it. The most important activation function of the last fifteen years is a single call, and its derivative is z > 0 — a boolean array, which NumPy happily multiplies as 1s and 0s.

Softmax must subtract its own max first. np.exp(800) is inf, and inf/inf is nan. Subtracting the row maximum before exponentiating changes nothing mathematically — the terms cancel — and stops the overflow dead. Skip this line and your loss becomes nan somewhere around epoch 12, which is a miserable afternoon.

import numpy as np

def relu(z):
    return np.maximum(0, z)

def softmax(z):
    z = z - z.max(axis=1, keepdims=True)      # the line people forget
    e = np.exp(z)
    return e / e.sum(axis=1, keepdims=True)

# Prove the overflow claim to yourself.
big = np.array([[800.0, 799.0, 1.0]])
naive = np.exp(big) / np.exp(big).sum()
print("naive: ", naive)                        # [[nan nan nan]]
print("stable:", softmax(big))                 # [[0.73 0.27 0.  ]]
print("sums to:", softmax(big).sum())

How the weights start matters. Initialise everything to zero and every hidden unit computes the same thing forever — they receive identical gradients and stay identical, so a 64-unit layer does the work of one. Initialise too large and the activations saturate. He initialisation — normal noise scaled by sqrt(2 / fan_in) — is the right default for ReLU and is what Keras uses when you write kernel_initializer="he_normal".

import numpy as np
from sklearn.datasets import load_digits

X, _ = load_digits(return_X_y=True)
X = X / 16.0
rng = np.random.default_rng(0)
H = 64

W1 = rng.normal(0, np.sqrt(2 / 64), (64, H)); b1 = np.zeros(H)
W2 = rng.normal(0, np.sqrt(2 / H), (H, 10));  b2 = np.zeros(10)

def relu(z): return np.maximum(0, z)
def softmax(z):
    z = z - z.max(axis=1, keepdims=True)
    e = np.exp(z); return e / e.sum(axis=1, keepdims=True)

p = softmax(relu(X[:4] @ W1 + b1) @ W2 + b2)
print("output shape:", p.shape)
print("row sums:", p.sum(axis=1))              # all 1.0 — they are probabilities
print("untrained confidence in each class:")
print(p[0].round(3))

Every row sums to 1 and the untrained network is spreading its bet across all ten digits at roughly 10% each. That is exactly right — it knows nothing yet.

3. The Loss, and Why Its Gradient Is So Clean

Cross-entropy loss for one image is -log(p[correct_class]). If the network gives the right digit a probability of 0.9, the loss is 0.105. If it gives it 0.01, the loss is 4.6. Confidence in the wrong answer is punished hard, and that asymmetry is the point.

Now the part that looks like a coincidence and is not. The gradient of cross-entropy loss with respect to z2 — the raw scores before softmax — is:

dz2 = p - Y

Predicted probabilities minus the one-hot truth. That is the entire derivative. The softmax derivative and the log derivative cancel each other almost completely, which is why every framework in existence fuses these two operations into one softmax_cross_entropy op rather than computing them separately. When you see that fused name in a stack trace, this cancellation is why.

import numpy as np

def softmax(z):
    z = z - z.max(axis=1, keepdims=True)
    e = np.exp(z); return e / e.sum(axis=1, keepdims=True)

Y = np.eye(10)[[3, 7]]                          # truth: digits 3 and 7
z2 = np.zeros((2, 10)); z2[0, 3] = 4.0; z2[1, 7] = 0.2   # confident, then not
p = softmax(z2)

loss = -np.log((p * Y).sum(axis=1))
print("per-image loss:", loss.round(3))         # confident guess costs far less
print("gradient dz2 for image 0:", (p - Y)[0].round(3))

Look at the gradient row: the correct class carries a negative number (push that score up) and every wrong class a small positive one (push those down). The whole of learning is that sign, repeated.

4. Backpropagation

Backprop is the chain rule with bookkeeping. You have dz2; you want the gradient for every weight. Walk backwards through the four forward lines, reversing each one:

The pattern repeats: to get a weight's gradient, matrix-multiply the input to that layer (transposed) by the gradient arriving at its output. To pass the gradient further back, multiply by the weights (transposed). Biases just sum, because a bias is added to every row identically.

Do not trust that table — or me. Check it against a numerical gradient. Nudge one weight by a tiny amount, see how much the loss moves, and compare that to what backprop claimed. If they agree to six decimal places, the maths is right.

import numpy as np
from sklearn.datasets import load_digits

X, y = load_digits(return_X_y=True)
X = X[:32] / 16.0
Y = np.eye(10)[y[:32]]

rng = np.random.default_rng(1)
W1 = rng.normal(0, 0.18, (64, 16)); b1 = np.zeros(16)
W2 = rng.normal(0, 0.35, (16, 10)); b2 = np.zeros(10)

def softmax(z):
    z = z - z.max(axis=1, keepdims=True)
    e = np.exp(z); return e / e.sum(axis=1, keepdims=True)

def loss_of(W1, W2):
    a1 = np.maximum(0, X @ W1 + b1)
    p = softmax(a1 @ W2 + b2)
    return -np.log((p * Y).sum(axis=1) + 1e-12).mean()

# Analytic gradient, straight from the table above.
z1 = X @ W1 + b1
a1 = np.maximum(0, z1)
p = softmax(a1 @ W2 + b2)
dz2 = (p - Y) / len(X)
dW2 = a1.T @ dz2
dz1 = (dz2 @ W2.T) * (z1 > 0)
dW1 = X.T @ dz1

# Check the weight with the largest gradient. Picking one at random is a trap:
# plenty of these weights have a gradient of EXACTLY zero, and 0.0 == 0.0 proves
# nothing. See below for why.
i, j = np.unravel_index(np.abs(dW1).argmax(), dW1.shape)
eps = 1e-6
Wp = W1.copy(); Wp[i, j] += eps
Wm = W1.copy(); Wm[i, j] -= eps
numeric = (loss_of(Wp, W2) - loss_of(Wm, W2)) / (2 * eps)

print("checking W1[%d, %d]" % (i, j))
print("backprop says: %+.6e" % dW1[i, j])
print("nudging says:  %+.6e" % numeric)
print("ratio:          %.6f" % (dW1[i, j] / numeric))
print("agree to 1e-9:", abs(dW1[i, j] - numeric) < 1e-9)

# Why some gradients are exactly zero: a pixel that is 0 in every image of the
# batch contributes nothing, so every weight leaving it stays untouched forever.
dead = np.where(X.max(axis=0) == 0)[0]
print()
print("pixels that are 0 across this batch:", len(dead), "of 64")
print("their weights get exactly zero gradient:", np.abs(dW1[dead]).max())

That check is worth keeping in your toolbox. When a custom layer misbehaves in a real framework years from now, gradient checking is still how you find out whether the bug is in your maths or your code.

5. Train It

Everything assembled. Mini-batches of 32, plain stochastic gradient descent, learning rate 0.5, 60 passes over the training set.

import numpy as np, time
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split

X, y = load_digits(return_X_y=True)
X = X / 16.0
Y = np.eye(10)[y]
Xtr, Xte, Ytr, Yte, ytr, yte = train_test_split(X, Y, y, test_size=0.25, random_state=0)

rng = np.random.default_rng(0)
H = 64
W1 = rng.normal(0, np.sqrt(2 / 64), (64, H)); b1 = np.zeros(H)
W2 = rng.normal(0, np.sqrt(2 / H), (H, 10));  b2 = np.zeros(10)

def softmax(z):
    z = z - z.max(axis=1, keepdims=True)
    e = np.exp(z); return e / e.sum(axis=1, keepdims=True)

def predict(A):
    return softmax(np.maximum(0, A @ W1 + b1) @ W2 + b2).argmax(axis=1)

lr, epochs, batch = 0.5, 60, 32
n = len(Xtr)
start = time.time()

for epoch in range(epochs):
    order = rng.permutation(n)
    for i in range(0, n, batch):
        idx = order[i:i + batch]
        xb, yb = Xtr[idx], Ytr[idx]

        z1 = xb @ W1 + b1                     # forward
        a1 = np.maximum(0, z1)
        p = softmax(a1 @ W2 + b2)

        dz2 = (p - yb) / len(idx)             # backward
        dW2 = a1.T @ dz2; db2 = dz2.sum(0)
        dz1 = (dz2 @ W2.T) * (z1 > 0)
        dW1 = xb.T @ dz1; db1 = dz1.sum(0)

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

    if epoch % 15 == 0 or epoch == epochs - 1:
        print("epoch %2d  train %.3f  test %.3f"
              % (epoch, (predict(Xtr) == ytr).mean(), (predict(Xte) == yte).mean()))

print("\ntrained in %.2fs" % (time.time() - start))
print("final test accuracy: %.3f" % (predict(Xte) == yte).mean())

pred = predict(Xte)
wrong = np.where(pred != yte)[0]
print("\n%d of %d wrong:" % (len(wrong), len(yte)))
for i in wrong:
    print("  true %d, guessed %d" % (yte[i], pred[i]))

It prints its own timing. On a recent laptop that is about a second, and the test score lands at 97.1%. Ninety-seven percent on hand-written digits, from a network whose every line you can point at.

Watch the two numbers in the epoch printout. Train accuracy climbing while test accuracy stalls is overfitting, live, in a model small enough to understand. That gap is the single most useful thing on the page.

And read the last lines of that output, the images it got wrong. The mistakes are never random. Thirteen out of 450, and they cluster: an 8 read as a 1 three times, a 7 and a 4 swapped once each way, and 1 is the wrong guess in five of the thirteen. At 8×8 an 8 loses the pinch in its waist and a 4 loses the gap in its top, so those pairs genuinely share most of their lit pixels. The model is not being stupid; it is being shown too few pixels. Across the whole dataset only 3 of the 64 pixels are zero in every single image — this data is denser than it looks, which is why 64 inputs get you to 97%.

6. When to Stop Doing This

You have now written a neural network by hand. Do it once. Then stop, and here is the honest reason.

Everything above is 64 inputs and 4,810 parameters. Scale the same code to a real problem and it falls over in four specific ways:

• No autodiff. Add a layer and you derive its backward pass by hand, again, correctly. PyTorch's loss.backward() does that for an arbitrary graph you invented five minutes ago.
• No GPU. NumPy is CPU-only. A model that takes 40 minutes on CPU takes 30 seconds on a GPU, and that difference decides whether you iterate ten times a day or once.
• Plain SGD only. Adam, learning-rate schedules, gradient clipping, weight decay — every one is another chunk of maths to get right, and the defaults in a framework are defaults for good reasons.
• No convolutions, no attention. Those are not small additions.

So: write it from scratch to understand it, use a framework to ship it. Anyone who tells you to build production models on raw NumPy is selling something.

Why PyTorch and TensorFlow are not on this page

They have no WebAssembly build, so they genuinely cannot run in a browser tab — not a policy we have chosen, a fact about the packages. Everything on this page runs here precisely because NumPy and scikit-learn do ship one.

For the framework lessons you will want a real machine:

pip install torch            # or: pip install tensorflow

That is a fair trade. The mechanics are what transfer between frameworks anyway, and the mechanics are what you just built. When nn.Linear(64, 64) shows up in the next lesson, you already know it holds a (64, 64) weight matrix and a length-64 bias, and you know precisely what its backward pass computes — because you wrote it.

Recap

• load_digits gives you 1,797 real labelled images with no download — the reason this page runs at all.
• The forward pass is four lines. Nearly every bug in it is a shape bug.
• Subtract the row max inside softmax, or meet nan at epoch 12.
• Never initialise weights to zero; He initialisation (sqrt(2/fan_in)) is the ReLU default.
• Cross-entropy after softmax has gradient p - Y. Frameworks fuse the two ops because of that cancellation.
• Backprop: dW = input.T @ grad_out, pass back with grad_out @ W.T, biases sum.
• Check any gradient you derive against a numerical nudge before trusting it.
• 97% test accuracy in about a second — then reach for a framework, because autodiff and a GPU are not optional at scale.