Neural Networks from Scratch

In 1969 Marvin Minsky and Seymour Papert published a book showing, with mathematical care, that a single-layer perceptron cannot learn something as simple as “exclusive or”: output 1 when exactly one of two switches is on. The fix was already known in principle — add a hidden layer between input and output — but nobody yet had a practical way to train one. By the end of this page you will have built that hidden layer by hand, seen it solve the problem in four lines of arithmetic, and written a general network class whose only missing piece is the learning procedure that took the field another seventeen years to popularise.

From one neuron to a layer

The track overview trained a single neuron: a weighted sum of inputs plus a bias, squashed by a sigmoid. Written as a formula:

As the maths page showed, w · x is a similarity score: the neuron fires strongly when the input points in the direction of its weights. The bias shifts the threshold — how similar the input must be before the neuron responds. The weights and biases together are the neuron’s parameters, the numbers training will adjust.

One detector is not much use. A layer is many neurons looking at the same input, each with its own weight vector. Stack those weight vectors as the columns of a matrix W and every neuron’s weighted sum comes out of a single matrix product. Feed a whole batch of inputs as the rows of X, and one line computes every neuron for every example:

A network is layers in sequence: the output of one becomes the input of the next. The layers between input and output are called hidden layers, because nobody tells them what to compute — they invent their own intermediate features during training. A network with several hidden layers is “deep”, which is all the word means in deep learning.

Why non-linearity matters: XOR

What if we left out φ? Then each layer would be X @ W + b, and a stack of them would be a product of matrices. But a product of matrices is just another matrix. However many linear layers you stack, the result is equivalent to a single linear layer. Let NumPy confirm it:

import numpy as np

rng = np.random.default_rng(42)
W1 = rng.normal(size=(2, 4))     # layer 1: 2 inputs -> 4 units
W2 = rng.normal(size=(4, 1))     # layer 2: 4 units -> 1 output
x = np.array([[0.5, -1.0]])

two_layers = (x @ W1) @ W2       # no activation in between
one_layer = x @ (W1 @ W2)        # a single 2 -> 1 matrix
print("two linear layers:", two_layers.ravel().round(4))
print("one merged layer: ", one_layer.ravel().round(4))
print("merged weights shape:", (W1 @ W2).shape)

two linear layers: [0.1267]
one merged layer:  [0.1267]
merged weights shape: (2, 1)

Matrix multiplication is associative — (xW₁)W₂ = x(W₁W₂) — so the four hidden units bought nothing. (Biases do not help either: they merge into a single bias.) A linear model can only draw straight boundaries: in two dimensions, one line separating “yes” from “no”.

Now plot the four XOR cases on a square: (0, 0) and (1, 1) should give 0; (0, 1) and (1, 0) should give 1. The 1s sit on opposite corners. No single straight line can put both 1-corners on one side and both 0-corners on the other. This is the limitation Minsky and Papert formalised in Perceptrons, and it fed the scepticism described in The First AI Winter.

Add a hidden layer with a non-linear activation and the problem dissolves. The weights below are the textbook solution from chapter 6 of Goodfellow, Bengio and Courville’s Deep Learning, using the ReLU activation, max(0, z):

import numpy as np

X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])   # the four XOR inputs

W1 = np.array([[1, 1],
               [1, 1]])          # both hidden units add the two inputs
b1 = np.array([0, -1])           # ...but the second needs both to be on
w2 = np.array([1, -2])           # output: first unit minus twice the second

H = np.maximum(0, X @ W1 + b1)   # hidden layer with ReLU
y = H @ w2
for inputs, hidden, out in zip(X, H, y):
    print(inputs, "-> hidden", hidden, "-> output", out)

[0 0] -> hidden [0 0] -> output 0
[0 1] -> hidden [1 0] -> output 1
[1 0] -> hidden [1 0] -> output 1
[1 1] -> hidden [2 1] -> output 0

Read the hidden column. The first unit counts how many inputs are on (0, 1, 1, 2). The second unit, thanks to its bias of −1 and the ReLU, fires only when both are on. The output takes “how many” and subtracts twice “both”, which is exactly XOR. In the hidden layer’s coordinates, (0, 1) and (1, 0) have landed on the same point (1, 0), and the problem has become linearly separable.

Activation functions

Any non-linear activation function breaks the collapse into a single linear map, but some train far better than others. Three have mattered most. Here are their values at a few points:

import numpy as np

z = np.array([-3.0, -1.0, 0.0, 1.0, 3.0])
sigmoid = 1 / (1 + np.exp(-z))
tanh = np.tanh(z)
relu = np.maximum(0, z)

print("z:      ", z)
print("sigmoid:", sigmoid.round(3))
print("tanh:   ", tanh.round(3))
print("relu:   ", relu)

z:       [-3. -1.  0.  1.  3.]
sigmoid: [0.047 0.269 0.5   0.731 0.953]
tanh:    [-0.995 -0.762  0.     0.762  0.995]
relu:    [0. 0. 0. 1. 3.]

Compare them tab by tab. The slope of each function matters as much as its shape, because on the next page training will multiply those slopes together.

The forward pass as matrix multiplication

Running inputs through a network to get outputs is called the forward pass. It is a short pipeline of matrix products and element-wise functions, and the fastest way to understand any architecture — including the Transformer — is to follow the shapes through it. Tap each stage of the network we are about to build: 2 inputs, 16 hidden units, 2 output classes, and a batch of 200 examples.

Notice what the weight shapes do not contain: the batch size. The same 82 parameters process one example or a million. That is also why the forward pass is so fast on a GPU — each stage is one large, regular, parallel operation.

The same pattern, scaled up, is a real model. A feed-forward block inside a GPT-2-sized Transformer takes a (batch, tokens, 768) tensor, multiplies by a (768, 3072) matrix, applies an activation, and multiplies by a (3072, 768) matrix to come back down. You already know how to read every step of that sentence.

Softmax outputs and a network in NumPy

For classification, the final layer produces one score per class, and a softmax turns each row of scores into probabilities, exactly as on the maths page. Now we can write the whole network as a small class. It has a constructor that creates the parameters, and a forward method that runs the pipeline above and remembers the intermediate results (we will need them on the next page). The dataset is make_moons from scikit-learn: two interleaving half-circles of points, generated on your machine, that no straight line can separate.

import numpy as np
from sklearn.datasets import make_moons

def softmax(Z):
    Z = Z - Z.max(axis=1, keepdims=True)       # row by row, for stability
    E = np.exp(Z)
    return E / E.sum(axis=1, keepdims=True)

class TwoLayerNet:
    def __init__(self, n_in, n_hidden, n_out, seed=42):
        rng = np.random.default_rng(seed)
        self.W1 = rng.normal(0, 1 / np.sqrt(n_in), size=(n_in, n_hidden))
        self.b1 = np.zeros(n_hidden)
        self.W2 = rng.normal(0, 1 / np.sqrt(n_hidden), size=(n_hidden, n_out))
        self.b2 = np.zeros(n_out)

    def forward(self, X):
        self.X = X
        self.H = np.tanh(X @ self.W1 + self.b1)    # (batch, n_hidden)
        logits = self.H @ self.W2 + self.b2        # (batch, n_out)
        self.P = softmax(logits)                   # each row sums to 1
        return self.P

X, y = make_moons(n_samples=200, noise=0.2, random_state=42)
net = TwoLayerNet(n_in=2, n_hidden=16, n_out=2)
P = net.forward(X)

print("X", X.shape, "-> H", net.H.shape, "-> P", P.shape)
print("first three rows of P:\n", P[:3].round(3))
print("parameters:", sum(a.size for a in (net.W1, net.b1, net.W2, net.b2)))
print("accuracy before training:", (P.argmax(axis=1) == y).mean())
print("cross-entropy before training:", round(-np.log(P[np.arange(len(y)), y]).mean(), 3))

X (200, 2) -> H (200, 16) -> P (200, 2)
first three rows of P:
 [[0.605 0.395]
 [0.352 0.648]
 [0.359 0.641]]
parameters: 82
accuracy before training: 0.55
cross-entropy before training: 0.71

Line by line, the new parts:

• keepdims=True keeps the row maxima and row sums as a (200, 1) column, so broadcasting subtracts and divides per row. Forgetting it is a classic bug: you would normalise over the whole batch instead.
• The initial weights are random, so that the 16 hidden units start different from each other (if they began identical they would stay identical). Scaling the spread by 1/√n_in keeps each weighted sum at a sensible size however many inputs feed it — an idea studied carefully by Xavier Glorot and Yoshua Bengio in 2010 and taken further in this track’s page on training well.
• Biases start at zero, which is the usual choice; the random weights already break the symmetry.
• The last line picks, for each example, the probability assigned to its true class (P[np.arange(len(y)), y]), and averages the cross-entropy over the batch.

The untrained network is guessing. Its accuracy of 0.55 is close to a coin flip, and its loss of 0.71 is close to −log 0.5 ≈ 0.693, the cross-entropy of a model that says 50/50 to everything. That is a useful sanity check whenever you build a classifier: before training, the loss should be roughly log(number of classes). If it is far higher, something is wrong with your initialisation or your loss.

What can a network represent?

The XOR network shows a hidden layer can carve one awkward region. How far does that go? There is a famous answer, often quoted and often over-quoted.

Consensus The universal approximation theorem, proved for sigmoid-type activations by George Cybenko (1989) and by Kurt Hornik, Maxwell Stinchcombe and Halbert White (1989), and extended by Moshe Leshno and colleagues (1993) to any activation that is not a polynomial (ReLU included), says: a network with one hidden layer can approximate any continuous function on a bounded region of input space as closely as you like — provided it has enough hidden units.

You can watch the “enough hidden units” part happen. Below, a hidden layer with random, untrained weights feeds an output layer whose weights are solved exactly by least squares. The target is the wiggly curve sin(2x) between −π and π.

import numpy as np

rng = np.random.default_rng(42)
x = np.linspace(-np.pi, np.pi, 200).reshape(-1, 1)
target = np.sin(2 * x).ravel()                     # a wiggly curve to imitate

for width in [2, 5, 20, 100]:
    W1 = rng.normal(0, 2, size=(1, width))         # random hidden layer, never trained
    b1 = rng.uniform(-np.pi, np.pi, size=width) * np.abs(W1).ravel()
    H = np.tanh(x @ W1 - b1)                       # (200, width)
    H = np.hstack([H, np.ones((200, 1))])          # plus an output bias
    w2, *_ = np.linalg.lstsq(H, target, rcond=None)  # best output weights
    error = np.abs(H @ w2 - target).max()
    print(f"hidden units {width:>3}: worst-case error {error:.4f}")

hidden units   2: worst-case error 1.1182
hidden units   5: worst-case error 0.8329
hidden units  20: worst-case error 0.0020
hidden units 100: worst-case error 0.0000

Each tanh unit is a soft step placed somewhere along the x axis; the output layer adds up scaled steps to build the curve. Two or five steps are too coarse to trace two full waves, while twenty already match it to within 0.002 on these sample points, and a hundred to within the four decimal places printed. That is the theorem’s intuition in miniature.

Debated What the theorem implies about deep learning is much less settled. Read it carefully and notice what it does not say:

• How many units? It promises existence, not economy. For some functions the required width grows exponentially with the number of inputs.
• Can training find them? It says such weights exist, not that gradient descent will reach them. Above, we cheated: the output layer was solved exactly.
• Will it generalise? It concerns fitting a known function, not learning from a finite, noisy sample. Our errors were measured only on the 200 points used for fitting, and only between −π and π; outside that range the approximation is worthless.

The practical lesson: a network with enough capacity can represent what you need far more often than not. Whether it learns it from your data is the question that matters, and that is where we go next. Our TwoLayerNet can compute; it cannot yet improve. On Backpropagation and Gradient Descent we give it a backward method and watch its loss fall from 0.71.

Key terms on this page

Sources & further reading

• Marvin Minsky & Seymour Papert, Perceptrons: An Introduction to Computational Geometry (MIT Press, 1969).
• George Cybenko, “Approximation by superpositions of a sigmoidal function”, Mathematics of Control, Signals and Systems 2 (1989); Kurt Hornik, Maxwell Stinchcombe & Halbert White, “Multilayer feedforward networks are universal approximators”, Neural Networks 2 (1989); Moshe Leshno, Vladimir Lin, Allan Pinkus & Shimon Schocken, “Multilayer feedforward networks with a nonpolynomial activation function can approximate any function”, Neural Networks 6 (1993).
• Ian Goodfellow, Yoshua Bengio & Aaron Courville, Deep Learning (MIT Press, 2016) — chapter 6, “Deep Feedforward Networks”, including the XOR example and universal approximation.
• Xavier Glorot, Antoine Bordes & Yoshua Bengio, “Deep sparse rectifier neural networks”, AISTATS (2011); Vinod Nair & Geoffrey Hinton, “Rectified linear units improve restricted Boltzmann machines”, ICML (2010).
• Matus Telgarsky, “Benefits of depth in neural networks”, COLT (2016); Ronen Eldan & Ohad Shamir, “The power of depth for feedforward neural networks”, COLT (2016).
• Grant Sanderson (3Blue1Brown), Neural Networks (video series, 2017 onwards) — chapter 1, “But what is a neural network?”.