A Network from Scratch: no framework, in this tab

Every code block on this page runs here, in your browser tab. 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.

You can check this yourself. Go 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 one short line, and it teaches you nothing about what went wrong when the loss suddenly becomes nan. If you write the backward pass by hand once, frameworks stop being a mystery: you will know what Keras is computing for you.

1. Real Data, No Download

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

It is not MNIST. MNIST is 28×28 and 70,000 images, but this is a smaller, older dataset from the UCI archive. For learning how a network works, that is an advantage: 64 inputs instead of 784 means the weight matrices are small enough to print and study.

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. Look closely 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 will 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                       # scales 0..16 down to 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, and 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. Keep track of the shapes. Almost every bug you will hit here is a shape bug, and NumPy will raise a clear error when it happens.

Two details are more important than they look:

ReLU is np.maximum(0, z). The most widely used activation function today is a single call. Its derivative is z > 0, a boolean array, which NumPy 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 does not change the result (the extra factor cancels out), and it prevents the overflow. Skip this line and, as soon as any score gets large enough, the softmax returns nan and the loss becomes nan with it. The small network in this lesson stays below that point, but larger inputs or a higher learning rate can reach it, so always subtract the max.

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  0.]]
print("stable:", softmax(big))                 # [[0.73 0.27 0.  ]]
print("sums to:", softmax(big).sum())

How the weights start matters. If you initialise everything to zero, 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. If the starting weights are too large, 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, because 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 correct, because it has learned 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. A confident wrong answer gets a much larger loss, and that is what we want.

The next result looks like a coincidence, but it 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 whole derivative. The softmax derivative and the log derivative cancel each other almost completely. This is why the major frameworks combine 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   # a confident guess, then an unsure one
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). Training is this sign, applied again and again.

4. Backpropagation

Backprop is the chain rule, applied carefully, step by step. You have dz2, and 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.

Don't just trust that table (or this lesson). 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 remembering. If a custom layer in a real framework gives strange results, gradient checking is still how you find out whether the bug is in your maths or your code.

5. Train It

Here is everything together. 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%. That is 97% on hand-written digits, from a network where you wrote every line.

Watch the two numbers in the epoch printout. When train accuracy keeps rising but test accuracy stops improving, that is overfitting, and here you can watch it happen in a model small enough to understand. That gap is the 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 really do share most of their lit pixels. The model's mistakes come from having too few pixels, not from a bug. Across the whole dataset only 3 of the 64 pixels are zero in every single image. This data has more information in it than it seems, 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. Here is why.

Everything above is 64 inputs and 4,810 parameters. If you use the same approach on a real problem, it fails in four specific ways:

• No autodiff. Every time you add a layer, you must derive its backward pass by hand again, without mistakes. PyTorch's loss.backward() does that automatically for any model you design.
• No GPU. NumPy is CPU-only. Large models train many times faster on a GPU, and that decides whether you can try ten ideas a day or one.
• Plain SGD only. Adam, learning-rate schedules, gradient clipping and weight decay each add more maths you must get right. Frameworks already include them, with well-tested defaults.
• No convolutions, no attention. Those are not small additions.

So: write it from scratch to understand it, use a framework to ship it. Don't build production models on raw NumPy.

Why PyTorch and TensorFlow are not on this page

They have no WebAssembly build, so they cannot run in a browser tab. This is a limit of the packages, not a choice we made. 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 the same in every framework, and the mechanics are what you just built. When Dense(64) (Keras) or nn.Linear(64, 64) (PyTorch) shows up later, you already know it holds a (64, 64) weight matrix and a length-64 bias, and you know exactly what its backward pass computes, because you wrote it.

What You Learned

• load_digits gives you 1,797 real labelled images with no download, which is why this page can run in the browser.
• The forward pass is four lines. Nearly every bug in it is a shape bug.
• Subtract the row max inside softmax, or large scores turn the loss into nan.
• 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. After that, use a framework, because autodiff and a GPU are not optional at scale.

Next: Building Neural Networks with Keras.