46 - Layered neural networks: forward pass, shapes and universal approximation

Imagine an assembly line. The raw material enters on one side; each station transforms it a little, passing it on to the next; at the end the finished product comes out. No station, on its own, does anything impressive, but the sequence does. A layered neural network works like this: the input (the features) goes through a series of layers, and each produces a new representation of it, better suited to the final task of the last one. In the previous lessons you studied the single neuron and saw that two layers solve XOR.

Now you really stack many neurons, and you learn the most important thing in everyday work with networks: following the shapes of the tensors precisely from one end to the other. The path leads to a complete hand-written forward pass, to counting the parameters, to why weight initialization is not a detail, and to a network that approximates a function with an error you can measure and drive down at will, to understand what the universal approximation theorem states (and what it does not).

The vocabulary of a network

• Layer: a group of neurons that all receive the same input and each produce an output. The output of a layer is a vector with as many numbers as it has neurons.
• Input layer: it computes nothing, it is simply the vector of features.
• Hidden layers: the intermediate layers. They are called this because their value is not observed directly, neither in the data nor in the output.
• Output layer: the last one, whose activation depends on the problem (linear for regression, sigmoid for binary classification, softmax for multiclass).
• Dense network (fully connected, feed-forward): every neuron of a layer is connected to all those of the next layer, and information flows only forward.
• Batch: a group of m examples processed together, in a single matrix multiplication. It is the reason tensors have one more dimension than "one example".

The six entries can be read in any order, but together they draw a single structure: the input comes in, the hidden layers transform it one after the other, the output layer converts it into what the problem needs. The batch is not a level of the network but a dimension of the data: when you later see the shape (m, n), m counts the examples and n the neurons of the layer, and the two dimensions must never be confused.

Notation and conventions on shapes

It is best to choose a convention and respect it everywhere. We use the usual one with NumPy, in which the first dimension is always the batch: a batch of m examples with n features is a matrix of shape (m, n). Layer l has n_l neurons. We denote the input by A⁰ = X, the weight matrix of layer l by W^l and its bias vector by b^l. The rule to remember is only one: W^l has shape (n_{l−1}, n_l), that is (neurons of the previous layer, neurons of this layer). With this convention the product A^{l−1} · W^l has shape (m, n_{l−1}) × (n_{l−1}, n_l) = (m, n_l): the inner dimensions cancel and the outer ones remain.

The drawing shows a 4 → 5 → 3 → 1 network: under each layer is the shape of the activation for a batch of m examples, and between two layers the shape of the weight matrix that connects them. If you can check by eye that every pair of shapes fits together, you are already able to spot the most frequent error in working with networks.

A note on libraries, verified with PyTorch 2.12.0 (see the section on shape and type traps): nn.Linear(n_in, n_out) stores the weights with shape (n_out, n_in), that is transposed with respect to our convention, and computes x · Wᵀ + b. The mathematical result is the same; if you extract the weights of a library to inspect them, check which of the two conventions it uses.

The forward pass

For each layer l = 1, ..., L two things are computed, exactly as for the single neuron but for all the neurons and all the examples together:

Z^l = A^{l−1} · W^l + b^l         pre-activation, shape (m, n_l)       [b^l is added to every row: broadcasting]
A^l = f_l(Z^l)                     activation,     shape (m, n_l)       [f_l applied element by element]

The output of the network is A^L. The name "forward pass" indicates that information flows from the input to the output; in the next lesson, to train, it will take the reverse path (backpropagation). Broadcasting deserves attention: b^l has shape (n_l,) and A^{l−1}W^l has shape (m, n_l); NumPy "stretches" the vector by adding it to each of the m rows. It is convenient, but it is also the mechanism through which a shape error can slip through silently.

A complete example with numbers, by hand

Take the smallest network with a hidden layer: 2 inputs, 2 hidden neurons with ReLU, 1 linear output. A single example, x = (1, 2). Weights and biases chosen for ease of calculation:

W¹ = [[ 0.5, −1.0],      b¹ = [0.0, 1.0]       shape (2, 2): row i = input i, column j = neuron j
      [ 0.25, 0.5]]

Z¹ = x · W¹ + b¹ = [1·0.5 + 2·0.25 + 0.0 ,  1·(−1.0) + 2·0.5 + 1.0] = [1.0, 1.0]
A¹ = ReLU(Z¹) = [1.0, 1.0]

W² = [[2.0],             b² = [0.5]             shape (2, 1)
      [−1.0]]

Z² = A¹ · W² + b² = 1.0·2.0 + 1.0·(−1.0) + 0.5 = 1.5       linear output: ŷ = 1.5

Work through the products calmly: column j of W¹ contains the weights arriving at neuron j, so Z¹_j is the dot product between x and column j, plus the bias. Let us check with NumPy, also showing the shapes of each step:

import numpy as np

x = np.array([[1.0, 2.0]])                      # a batch of m = 1 example: shape (1, 2)
W1 = np.array([[0.5, -1.0], [0.25, 0.5]])       # (2, 2)
b1 = np.array([0.0, 1.0])                       # (2,)
W2 = np.array([[2.0], [-1.0]])                  # (2, 1)
b2 = np.array([0.5])                            # (1,)

Z1 = x @ W1 + b1                                # (1, 2) @ (2, 2) -> (1, 2), then + b1 with broadcasting
A1 = np.maximum(0, Z1)                          # ReLU
Z2 = A1 @ W2 + b2                               # (1, 2) @ (2, 1) -> (1, 1)

print("Z1 =", Z1, " shape", Z1.shape)
print("A1 =", A1, " shape", A1.shape)
print("output =", Z2, " shape", Z2.shape)

Z1 = [[1. 1.]]  shape (1, 2)
A1 = [[1. 1.]]  shape (1, 2)
output = [[1.5]]  shape (1, 1)

What the code does. The @ operator is the matrix product. The Z1 line verifies the hand calculation: 1.0 and 1.0. np.maximum(0, Z1) applies ReLU to every element. The output 1.5 matches the calculation. Get used to printing shapes: in a larger program, the line print(Z1.shape) is the cheapest debugging tool there is.

Counting the parameters

The parameters are the numbers that training has to learn: all the weights and all the biases. Layer l has a matrix W^l with n_{l−1}·n_l elements and a bias with n_l elements, so:

parameters of layer l = n_{l−1} · n_l + n_l = n_l · (n_{l−1} + 1)
total parameters = sum over all layers

For the network in the drawing (4 → 5 → 3 → 1): 5·(4+1) = 25, then 3·(5+1) = 18, then 1·(3+1) = 4, for a total of 47 parameters. The number of parameters is the measure of a network's "capacity", and it determines memory and computation time. A 784 → 128 → 64 → 10 network (a classifier for 28×28 images) has 128·785 + 64·129 + 10·65 = 100 480 + 8 256 + 650 = 109 386 parameters: more than the examples of many datasets, and a first reason why overfitting is a concrete risk.

A complete network by hand: 4 → 8 → 3 on the Iris dataset

Let us put everything in a small class, with a reasonable initialization (we will see shortly why) and a forward pass for a 3-class problem, with softmax at the output. We do not train it yet: we want to check the shapes and the correctness of an untrained model, a check that is done before any serious training.

from sklearn.datasets import load_iris
from sklearn.preprocessing import StandardScaler

Xi, yi = load_iris(return_X_y=True)
Xi = StandardScaler().fit_transform(Xi)                  # for a shape check; in a real project, fit on the training set only

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

class Network:
    def __init__(self, sizes, seed=0):
        rng = np.random.default_rng(seed)
        self.W, self.b = [], []
        for n_in, n_out in zip(sizes[:-1], sizes[1:]):
            self.W.append(rng.normal(size=(n_in, n_out)) * np.sqrt(2.0 / n_in))   # "He" scale: see below
            self.b.append(np.zeros(n_out))

    def forward(self, X, verbose=False):
        A = X
        for l, (W, b) in enumerate(zip(self.W, self.b), start=1):
            Z = A @ W + b
            last = (l == len(self.W))
            A = softmax(Z) if last else np.maximum(0, Z)                 # ReLU hidden, softmax at the output
            if verbose:
                print("layer %d: W %s  Z %s  -> A %s" % (l, W.shape, Z.shape, A.shape))
        return A

    def n_parameters(self):
        return sum(W.size + b.size for W, b in zip(self.W, self.b))

network = Network([4, 8, 3])
P = network.forward(Xi, verbose=True)
loss = -np.log(P[np.arange(len(yi)), yi]).mean()
print("parameters:", network.n_parameters(), "  (8*(4+1) + 3*(8+1) = %d)" % (8 * 5 + 3 * 9))
print("the rows of P sum to 1:", np.allclose(P.sum(axis=1), 1.0), "  shape of P:", P.shape)
print("untrained cross-entropy: %.4f   (ln 3 = %.4f)" % (loss, np.log(3)))

layer 1: W (4, 8)  Z (150, 8)  -> A (150, 8)
layer 2: W (8, 3)  Z (150, 3)  -> A (150, 3)
parameters: 67   (8*(4+1) + 3*(8+1) = 67)
the rows of P sum to 1: True   shape of P: (150, 3)
untrained cross-entropy: 1.5728   (ln 3 = 1.0986)

What the code does. The class stores a list of W matrices and b vectors, one per layer, built from the list of sizes [4, 8, 3]: zip(sizes[:-1], sizes[1:]) produces the pairs (4, 8) and (8, 3), that is (input neurons, output neurons) of each layer. forward goes through the layers: A @ W + b, then ReLU for the hidden ones and softmax for the last. With verbose=True it prints the shapes of each step: 150 Iris examples, a (150, 4) matrix that becomes (150, 8) and then (150, 3). The cross-entropy of an untrained network has as its reference ln 3 ≈ 1.0986, the value you would get if the model assigned each of the 3 classes an equal probability of 1/3.

Here we get 1.57: a bit higher, because with He initialization the output scores are not exactly zero and the model is "confident by chance" about some wrong class. That is fine: the order of magnitude is right. If the initial loss were very different (for example 15, or nan) you would have a bug or a wrong initialization, and it is better to notice before spending hours of training. It is one of the most useful sanity checks in the whole practice.

Weight initialization is not a detail

Why does the class use np.sqrt(2.0 / n_in) as the weight scale? Because at each layer the signal is multiplied by a weight matrix. If the weights are too small, the signal fades layer after layer until it vanishes; if they are too large, it explodes. In both cases there is no useful information left to learn from, nor a gradient signal coming back. Let us measure it: we pass 1000 vectors of 100 random numbers through 10 layers of 100 neurons, with different weight scales, and look at the standard deviation of the activations at the first, fifth and tenth layer.

rng = np.random.default_rng(0)
X100 = rng.normal(size=(1000, 100))

scale = {"std 0.01": lambda n: 0.01,
         "std 1": lambda n: 1.0,
         "Xavier  sqrt(1/n)": lambda n: np.sqrt(1.0 / n),
         "He  sqrt(2/n)": lambda n: np.sqrt(2.0 / n)}

print("weight scale          activation   std layer 1    std layer 5    std layer 10")
for name, weight_scale in scale.items():
    for act in ("tanh", "relu"):
        A, std = X100, []
        for _ in range(10):
            W = rng.normal(size=(100, 100)) * weight_scale(100)
            Z = A @ W
            A = np.tanh(Z) if act == "tanh" else np.maximum(0, Z)
            std.append(A.std())
        print("%-20s  %-10s   %12.5f   %12.5f   %13.5f" % (name, act, std[0], std[4], std[9]))

weight scale          activation   std layer 1    std layer 5    std layer 10
std 0.01              tanh              0.09867        0.00001         0.00000
std 0.01              relu              0.05850        0.00000         0.00000
std 1                 tanh              0.95916        0.95650         0.95732
std 1                 relu              5.82850    15393.49385   323956922.93796
Xavier  sqrt(1/n)     tanh              0.62508        0.31668         0.22479
Xavier  sqrt(1/n)     relu              0.58026        0.15518         0.02418
He  sqrt(2/n)         tanh              0.71849        0.56174         0.55186
He  sqrt(2/n)         relu              0.82523        1.04275         1.41455

What the code does. It generates the input data, then for each combination of weight scale and activation it stacks 10 random layers and records how "large" the activations are (their standard deviation). With weights of standard deviation 0.01 the signal vanishes already at the fifth layer (practically zero): training cannot start. With standard deviation 1 and ReLU the activations grow by enormous factors (on the order of 10⁸ at the tenth layer): the network explodes.

With tanh and standard deviation 1 the outputs stick to the values ±1 (saturation, with an almost zero derivative). The Xavier initializations (or Glorot; here in the simplified variant with variance 1/n_in, while the original by Glorot and Bengio uses 2/(n_in + n_out), designed for tanh and sigmoid) and He (with variance 2/n_in, designed for ReLU) keep the activations in a reasonable range: with He and ReLU the standard deviation stays between 0.8 and 1.4 for all 10 layers.

The factor 2 of He compensates for the fact that ReLU zeroes about half of the activations; the table confirms it, because Xavier applied to a ReLU lets the standard deviation fall from 0.58 to 0.02 in ten layers. Each activation has its own suitable initialization.

The universal approximation theorem

How complex can the function that a network is able to represent be? A fundamental answer is the universal approximation theorem (Cybenko, 1989, for sigmoid activations; Hornik, 1991, for a broader class of activations). Stated in words: a network with a single hidden layer and a finite, but sufficiently large, number of neurons with a non-polynomial activation can approximate any continuous function on a bounded set with arbitrary precision.

The statement is powerful, but it must be read carefully. What it states: that for every continuous function there exists a network (with one hidden layer) that approximates it with the error you want. What it does not state: (1) how many neurons are needed, which can be an enormous number; (2) how to find the weights, nor that gradient descent finds them; (3) that the network generalizes: approximating the training data well is not approximating new data well; (4) that a one-layer network is the efficient way to represent the function. It is an existence theorem, not a recipe.

An explicit construction with ReLU

Instead of trusting it, build by hand a network with a single hidden layer of ReLUs that approximates a function, and measure the error. The idea: a continuous function is approximated by a piecewise linear curve with H segments, and a piecewise linear curve can be written as a sum of ReLUs. Choose H+1 "knots" t₀ < t₁ < ... < t_H and let y_k = f(t_k). The slope of segment k (between t_k and t_{k+1}) is s_k = (y_{k+1} − y_k)/(t_{k+1} − t_k). The piecewise linear curve that interpolates the points is:

g(x) = y₀ + s₀ · (x − t₀) + Σ_{k=1..H−1} (s_k − s_{k−1}) · ReLU(x − t_k)

Why is it right? For x between t₀ and t₁ all the terms ReLU(x − t_k) with k ≥ 1 are zero, and what remains is the line y₀ + s₀(x − t₀), that is the first segment. When x passes t₁ the term ReLU(x − t₁) "switches on", adding slope (s₁ − s₀): the total slope becomes s₀ + (s₁ − s₀) = s₁, exactly that of the next segment; and so on at every knot. Each ReLU is a hidden neuron with weight 1, bias −t_k and an output weight (s_k − s_{k−1}). With H neurons you have a piecewise linear curve of H segments.

def piecewise_network(f, a, b, H):
    t = np.linspace(a, b, H + 1)                      # knots
    y = f(t)
    s = np.diff(y) / np.diff(t)                       # slopes of the H segments
    output_weights = np.diff(s)                       # (s_k - s_{k-1}) for k = 1..H-1
    def network(x):
        output = y[0] + s[0] * (x - t[0])             # initial linear term
        for k in range(1, H):
            output = output + output_weights[k - 1] * np.maximum(0, x - t[k])   # one ReLU hidden neuron per knot
        return output
    return network

xs = np.linspace(0, 2 * np.pi, 2001)
print("H (neurons)   maximum error on [0, 2*pi]   ratio to the previous error")
previous = None
for H in (2, 4, 8, 16, 32, 64):
    error = np.abs(piecewise_network(np.sin, 0, 2 * np.pi, H)(xs) - np.sin(xs)).max()
    print("%11d   %27.5f   %s" % (H, error, "-" if previous is None else "%.2f" % (previous / error)))
    previous = error

H (neurons)   maximum error on [0, 2*pi]   ratio to the previous error
          2                       1.00000   -
          4                       0.21051   4.75
          8                       0.07038   2.99
         16                       0.01885   3.73
         32                       0.00479   3.93
         64                       0.00120   3.99

What the code does. piecewise_network receives the function to approximate (here the sine on [0, 2π]), computes the knots, the slopes and the output weights with the formula above, and returns a function that evaluates the "network" with H ReLU hidden neurons. It is evaluated on 2001 points and compared with the true sine. The maximum error decreases as H grows: from about 0.21 with 4 neurons to 0.0012 with 64. The ratio to the previous error tends to 4 when H doubles: the error decreases as 1/H², as expected from a piecewise linear interpolation of a smooth function, whose error is proportional to the square of the width of the segments.

You therefore have a constructive, numerical demonstration that a single hidden layer with enough neurons gets where it wants. You also have the measure of the price: to obtain an error of one thousandth you need tens of neurons even for a simple function like the sine. In high dimensions, a grid construction like this one would require a number of neurons that grows exponentially with the dimension (reasoning, not measured here).

Depth or width?

The theorem speaks of a wide, shallow network. In practice, networks with more layers of moderate width are used. The reason, only partly theoretical, is that many interesting functions have a hierarchical structure (pixels → edges → shapes → objects), and the later layers can reuse the representations built by the earlier ones.

There are functions that a deep network represents with few neurons and that a shallow network would represent only with an exponentially larger number: Telgarsky proves it for a class of nodes that includes ReLU (in the sources). The price of depth is more delicate training: gradients that vanish or explode (seen above), more hyperparameters, more sensitivity to initialization. The techniques you will see in the last lessons (correct initialization, ReLU, batch normalization, residual connections) exist to make deep networks trainable.

Shape traps, seen in NumPy

Shape errors are the most common in working with networks. We reuse the Network class to list the shapes of its parameters and try to pass it a batch with the wrong number of features, to see how NumPy reports the problem.

print("parameter shapes:", [tuple(a.shape) for W, b in zip(network.W, network.b) for a in (W, b)])

try:
    network.forward(np.ones((5, 3)))              # the batch has 3 features instead of 4
except ValueError as e:
    print("shape error:", e)

logit = (np.maximum(0, np.ones((5, 4)) @ network.W[0] + network.b[0])) @ network.W[1] + network.b[1]
print("output without softmax: shape", logit.shape, " - these are 'logits', not probabilities")

parameter shapes: [(4, 8), (8,), (8, 3), (3,)]
shape error: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 4 is different from 3)
output without softmax: shape (5, 3)  - these are 'logits', not probabilities

What the code does. The first line confirms the parameters of the 4 → 8 → 3 network: two matrices (4, 8) and (8, 3), two biases (8,) and (3,), 67 numbers in all, like the 8·5 + 3·9 computed above. The try block reproduces the most frequent error: a batch with 3 features cannot be multiplied with a weight matrix that expects 4, and the matrix product cannot "guess".

The message, unlike broadcasting, is explicit. The last line computes by hand the forward pass up to the scores of the last layer: they are logits, and to get probabilities you need the softmax. Libraries with "with logits" loss functions include the softmax in the loss and want the raw scores; if you apply it twice the model trains worse without reporting anything.

The same network in PyTorch

Whoever moves to PyTorch will find the same shape traps and, in addition, one about the data type: NumPy normally uses float64, while PyTorch initializes the weights in float32. The listing was run with PyTorch 2.12.0 on CPU (CUDA_VISIBLE_DEVICES="", 2 threads), twice, with identical output. The fixed seed concerns only the initial weights: the three published outputs do not depend on the random values.

import torch

torch.manual_seed(0)
net_t = torch.nn.Sequential(torch.nn.Linear(4, 8), torch.nn.ReLU(), torch.nn.Linear(8, 3))
print("parameter shapes:", [tuple(p.shape) for p in net_t.parameters()])
print("total number of parameters:", sum(p.numel() for p in net_t.parameters()))

# Trap 1: data type. NumPy uses float64 by default, PyTorch float32.
try:
    net_t(torch.tensor(np.random.rand(5, 4)))
except RuntimeError as e:
    print("type error:", str(e)[:70])

# Trap 2: shape. The batch has 3 features instead of 4.
try:
    net_t(torch.rand(5, 3))
except RuntimeError as e:
    print("shape error:", str(e)[:70])

logit = net_t(torch.rand(5, 4))
print("correct output: shape", tuple(logit.shape), " - before the softmax these are 'logits', not probabilities")

parameter shapes: [(8, 4), (8,), (3, 8), (3,)]
total number of parameters: 67
type error: mat1 and mat2 must have the same dtype, but got Double and Float
shape error: mat1 and mat2 shapes cannot be multiplied (5x3 and 4x8)
correct output: shape (5, 3)  - before the softmax these are 'logits', not probabilities

What the code does. nn.Sequential builds the same 4 → 8 → 3 network as the NumPy block. The parameters have shapes (8, 4), (8,), (3, 8), (3,): Linear keeps them as (n_out, n_in), transposed with respect to our convention, and the total is 67, as in NumPy. The first try reproduces the type error: np.random.rand produces float64, the weights are float32, and the product between different types gives an explicit error; it is fixed with torch.tensor(x, dtype=torch.float32). The second is a shape error, the same as in the NumPy version, with a differently worded message. The correct output has shape (5, 3) and is made of logits.

Trade-offs and real-world traps

• Wrong shapes that pass silently: broadcasting accepts (m,) and (m, 1) in combinations that produce an (m, m) matrix instead of a vector. A classic example: subtracting a vector y of shape (m,) from a prediction of shape (m, 1) gives an (m, m) matrix, and the mean of the error has a plausible but wrong value. Always check the shapes after every non-obvious operation.
• Transposed weight matrices: (n_in, n_out) in the convention of this article, (n_out, n_in) in other libraries (see the note above). When you load or compare weights between the two conventions, remember to transpose.
• Initialization at zero or too small/large: a signal that vanishes or explodes (see above). Use He for ReLU, Xavier for tanh and sigmoid; many libraries already choose a default initialization (check which), but if you write by hand you have to do it yourself.
• Non-standardized inputs: the activations of the first layer depend on the scale of the features; with inputs in the order of thousands, even a good initialization will saturate tanh and sigmoid.
• Applying softmax twice (in the network and in the loss).
• Confusing approximation with generalization: the universal approximation theorem says nothing about new data. A large enough network also approximates the noise in the training set.
• Too much capacity: 109 386 parameters for a small classifier on a few thousand examples mean almost certain overfitting without regularization.

The seven traps are of two kinds. The first four (shapes, transpositions, initialization, input scale) produce a wrong number without giving any error: they are discovered only with explicit checks, such as printing the shapes and looking at the initial loss. The last three (double softmax, approximation confused with generalization, excessive capacity) concern judgment about what the model does: they are prevented with validation and with a careful reading of what a theorem really states.

Three worked exercises, to check that you understood

Exercise 1. How many parameters does a 20 → 64 → 32 → 1 network have? Apply the formula n_l·(n_{l−1} + 1) to each layer: 64·21 = 1 344, then 32·65 = 2 080, then 1·33 = 33; total 3 457. Note that almost two thirds of the parameters are in the second layer: it is the connections between the two widest layers that weigh the most.

Exercise 2. What shapes go through a 10 → 16 → 4 network with a batch of 32 examples? The input is (32, 10); the first layer has W¹ (10, 16) and produces (32, 16); the second has W² (16, 4) and produces (32, 4). The biases have shapes (16,) and (4,) and are added to every row.

Exercise 3. Why does a mean squared error seem to "work" but give a wrong value? The classic case is a prediction of shape (m, 1) compared with a label vector of shape (m,). Let us check it:

pred = np.array([[1.0], [2.0], [3.0], [4.0], [5.0]])         # (5, 1): network output
y = np.array([1.0, 2.0, 3.0, 4.0, 5.0])                      # (5,): perfect labels, perfect prediction

print("shape of the difference:", (pred - y).shape, " (expected: (5, 1) or (5,))")
print("wrong MSE (broadcast to 5x5): %.3f" % ((pred - y) ** 2).mean())
print("correct MSE with pred.ravel():   %.3f" % ((pred.ravel() - y) ** 2).mean())

n_20_64_32_1 = sum(n_out * (n_in + 1) for n_in, n_out in zip([20, 64, 32], [64, 32, 1]))
print("parameters of 20-64-32-1:", n_20_64_32_1)

shape of the difference: (5, 5)  (expected: (5, 1) or (5,))
wrong MSE (broadcast to 5x5): 4.000
correct MSE with pred.ravel():   0.000
parameters of 20-64-32-1: 3457

With the perfect prediction the error must be zero, but the calculation with the wrong shapes returns 4.0: pred - y broadcasts (5, 1) against (5,) into a (5, 5) matrix that compares every prediction with every label. No error, no warning, a plausible and wrong number. The cure is discipline: give predictions and labels the same shape, and print .shape before every loss computation. The last line also verifies the first exercise: 3 457.

Frequently asked questions

How many layers and how many neurons should I use? There is no formula. You start from a small, simple architecture, check that it can overfit a tiny subset of the data (a sanity check: if it cannot, there is a bug), and then increase the capacity until validation improves. For tabular data one or two hidden layers are often enough.

Why does the batch size appear in all the shapes? Because the GPU (and NumPy) are much more efficient at multiplying large matrices than at repeating vector multiplications. Processing m examples together is mathematically identical to processing them one at a time, but much faster.

Is the bias the same for all the examples? Yes: it is a parameter of the neuron, not of the example. This is why it has shape (n_l,) and is added to every row of the batch.

Does a network with more hidden neurons always learn better? No: the extra capacity is converted into a better fit to the training data, which can mean overfitting. It is always evaluated on validation.

What you can do now

You can write a complete forward pass with NumPy and check by eye that the shapes fit together, count the parameters of a network with the formula n_l·(n_{l−1} + 1) and check the initial loss before training. You can choose the weight scale for tanh and for ReLU and recognize from the experiment when the signal vanishes or explodes, and you know what the universal approximation theorem says and what it does not say, having measured the error of the piecewise linear network. The previous lesson is the perceptron and activation functions; in the next one, backpropagation derived in full, you will see how a layered network is trained.

Other lessons in the series that link to this one: A neural network from scratch in NumPy, then in PyTorch.