Neural Networks from Scratch

Suppose CodeAssist receives a code-completion request with a prompt-size score of 6 and an active queue-wait score of 8. A quiet request has values 2 and 3. We want a scoring function that flags the first request as high risk before worker threads time out. Raw numbers don't make an operational decision on their own, though. What computation should sit between those incoming features and a routing decision?

You already have the ingredients for a neural network: feature vectors and dot products from the vectors and matrices lesson. A network arranges those computations into layers, then uses learned weights and biases to turn inputs into outputs. We'll freeze the weights first so you can check the arithmetic before training starts changing them.

One neuron scores a single request

Start with one scoring unit, or neuron. Multiply each input feature by its own weight, then add their contributions together. Both inputs use illustrative zero-to-ten scales; raw token counts and seconds would have different units and ranges. Before reading the sum, check which feature contributes more:

The contributions total 6.6. Now add a bias, an offset of -4.0, to get 2.6. This weighted sum plus bias is an affine map. A purely linear map satisfies $f(\alpha u + \beta v) = \alpha f(u) + \beta f(v)$, which forces $f(0) = 0$. Without bias, zero input features give a weighted sum of zero before activation, which sigmoid would map to 0.5. The bias lets the model change that baseline and shift the decision boundary away from the origin.

In vector notation, the feature vector $x$ holds [6, 8], the weight vector $w$ holds [0.7, 0.3], and $b = -4.0$:

The superscript $\top$ denotes a transpose, making $w^\top x$ the dot product. An activation function transforms the resulting pre-activation $z$. For this binary risk scorer, pass $z$ through the sigmoid function, which maps 2.6 to roughly 0.931:

Here $e$ is Euler's constant, roughly 2.718. Before entering sigmoid, $z$ is called a logit. Hidden-layer pre-activations aren't generally logits because they don't feed directly into a logistic probability mapping. Mathematically, sigmoid stays strictly between zero and one for finite inputs. Floating-point implementations can round extreme outputs to zero or one; we'll address that numerical distinction below.

These hand-chosen parameters make the score 0.931. A production probability claim needs calibration evidence on representative, verified request outcomes; training alone leaves calibration unmeasured.^{[1]Reference 1On Calibration of Modern Neural Networkshttps://arxiv.org/abs/1706.04599}

Compute this scoring path once with Python's standard library. The snippet below prints the logit and sigmoid output separately so you can verify each arithmetic step:

from math import exp

prompt_size = 6.0
queue_wait = 8.0
w_prompt, w_queue, b = 0.7, 0.3, -4.0

z = prompt_size * w_prompt + queue_wait * w_queue + b
score = 1.0 / (1.0 + exp(-z))

print(f"logit:         {z:.3f}")
print(f"sigmoid score: {score:.3f}")
assert abs(z - 2.6) < 1e-9

logit:         2.600
sigmoid score: 0.931

View larger

A bias shifts the intervention threshold without altering feature values. Keep the same weights and compare a quiet request with our high-load request. Notice how each logit changes while the weights stay fixed:

from math import exp

w_prompt, w_queue, b = 0.7, 0.3, -4.0
requests = {
    "quiet": (2.0, 3.0),
    "high load": (6.0, 8.0),
}

for name, (prompt_size, queue_wait) in requests.items():
    z = prompt_size * w_prompt + queue_wait * w_queue + b
    score = 1.0 / (1.0 + exp(-z))
    print(f"{name:9s} logit={z:5.2f} score={score:.3f}")

quiet     logit=-1.70 score=0.154
high load logit= 2.60 score=0.931

At a score threshold of 0.5, the high-load request crosses the alert line while the quiet request stays below it. Equivalently, compare each logit with zero. The bias subtracts the exact same amount from both logits, shifting the decision boundary without touching feature weights.

Linear maps collapse without non-linearities

One neuron scores a weighted sum. Why can't a stack of linear layers do anything that one layer cannot?

Consider a scalar example: the first layer doubles an input and adds 1; the second triples that result and subtracts 2. Starting from 4, the first layer gives 2(4) + 1 = 9. The second gives 3(9) - 2 = 25. Notice that multiplying 4 directly by 6 and adding 1 also gives 25. For any input $x$, $3(2x + 1) - 2 = 6x + 1$. The two layers collapsed into one.

The exact same collapse happens with vectors and matrices. Let $h$ be a hidden vector. The first layer applies weights $W_1$ and biases $b_1$; the second applies $W_2$ and $b_2$:

Now substitute the expression for $h$ directly into the output equation:

Define $W_{\text{merged}} = W_2 W_1$ and $b_{\text{merged}} = W_2 b_1 + b_2$. The two layers collapse into:

The merged map remains affine: stacking affine layers preserves this form, while narrow hidden layers can further constrain its rank. An affine decision boundary cannot represent XOR, where $(0, 1)$ and $(1, 0)$ are positive but $(0, 0)$ and $(1, 1)$ are negative. No single hyperplane separates those diagonal pairs.

A non-linear activation placed between layers breaks that linear superposition, allowing the network to carve curved and piecewise-linear decision surfaces.^{[2]Reference 2Deep Learning.https://www.deeplearningbook.org/}

Test this collapse directly in NumPy. Without an activation, the two-layer calculation matches the merged affine map down to machine precision. Inserting ReLU breaks that equivalence:

import numpy as np

x = np.array([1.0, 4.0])  # small prompt, long queue
w1 = np.array([[0.3, -0.1], [0.1, 0.6]])
b1 = np.array([0.0, 0.0])
w2 = np.array([0.4, 0.1])
b2 = -0.5

z1 = w1 @ x + b1
two_affine = w2 @ z1 + b2
merged_affine = (w2 @ w1) @ x + (w2 @ b1 + b2)
with_relu = w2 @ np.maximum(0.0, z1) + b2

print("hidden pre-activation:", z1)
print("affine equals merged:", np.allclose(two_affine, merged_affine))
print(f"two affine output:    {float(two_affine):.3f}")
print(f"with ReLU output:     {float(with_relu):.3f}")

# Every affine map satisfies f(x) + f(-x) = 2 f(0).
def forward(point):
    return w2 @ np.maximum(0.0, w1 @ point + b1) + b2

print("affine symmetry holds:", np.isclose(forward(x) + forward(-x), 2 * forward(np.zeros_like(x))))

hidden pre-activation: [-0.1  2.5]
affine equals merged: True
two affine output:    -0.290
with ReLU output:     -0.250
affine symmetry holds: False

The first hidden pre-activation is -0.1. ReLU clips it to zero, changing the hidden vector from [-0.1, 2.5] to [0.0, 2.5]. The output changes from -0.290 to -0.250. One changed output only distinguishes these two calculations; the symmetry check goes further. This network violates an identity satisfied by every affine map, so no single affine map reproduces it for all inputs. Other choices of weights can still make a network with ReLU behave affinely.

Comparing modern activation functions

The activation you choose changes both what a layer can represent and how gradients flow through it during training. Compare these five on the same pre-activation:

1. Sigmoid: $\sigma(z) = 1 / (1 + e^{-z})$, mapping $(-\infty, \infty)$ to $(0, 1)$. Its derivative is $\sigma'(z) = \sigma(z)(1 - \sigma(z))$, reaching a peak of only $0.25$ at $z = 0$. Sigmoid is useful for binary classification outputs. Its local derivative is at most 0.25 and shrinks in the tails, which can contribute to vanishing gradients in deep hidden stacks. For one neuron on one example, $\partial\mathcal{L}/\partial w_j = \delta x_j$: nonnegative inputs make these weight gradients share the sign of $\delta$. That isn't a rule for every weight in a layer or a batch update.

2. Tanh (Hyperbolic Tangent): $\tanh(z) = (e^z - e^{-z}) / (e^z + e^{-z}) = 2\sigma(2z) - 1$, mapping $(-\infty, \infty)$ to $(-1, 1)$. Its derivative is $\tanh'(z) = 1 - \tanh^2(z)$, peaking at $1.0$ at $z = 0$. Tanh is odd and permits both negative and positive outputs. The mean of actual hidden activations depends on their input distribution. Its derivative shrinks as $|z|$ grows, so it can also contribute to vanishing gradients.

3. ReLU (Rectified Linear Unit): $\text{ReLU}(z) = \max(0, z)$, mapping $(-\infty, \infty)$ to $[0, \infty)$. Its derivative is exactly $1$ for $z > 0$ and $0$ for $z < 0$. ReLU's positive branch doesn't shrink the gradient locally, but products of weights and other network operations can still make gradients vanish or explode. A unit that is negative on every training input receives no gradient through its ReLU branch and can remain inactive. Changes in upstream features or optimization dynamics can revive it, so inactivity isn't an unconditional permanent state. ReLU is continuous at zero, where its ordinary derivative doesn't exist; implementations commonly choose a zero backward value.

4. GELU (Gaussian Error Linear Unit): $\text{GELU}(z) = z \Phi(z) = z \cdot \frac{1}{2}[1 + \text{erf}(z / \sqrt{2})]$, mapping $(-\infty, \infty)$ to roughly $[-0.17, \infty)$. Hendrycks and Gimpel introduced GELU to smoothly gate values by their probabilistic magnitude under standard Gaussian noise.^{[3]Reference 3Gaussian Error Linear Units (GELUs).https://arxiv.org/abs/1606.08415} Unlike ReLU's corner, GELU is smooth and dips slightly negative ($\approx -0.17$ near $z \approx -0.75$). Its derivative is zero at its minimum and becomes very small far into the negative tail, so gradient flow still depends on the inputs and surrounding layers. GELU appears in BERT and several GPT and Vision Transformer architectures; exact or approximate variants depend on the checkpoint.

5. Swish / SiLU (Sigmoid Linear Unit): $\text{Swish}(z) = z \cdot \sigma(\beta z)$, and when $\beta = 1$, it's called SiLU: $\text{SiLU}(z) = z \cdot \sigma(z)$, mapping to roughly $[-0.278, \infty)$. Swish was studied through automated activation search; its $\beta = 1$ form also appears as SiLU in the earlier GELU paper.^{[3]Reference 3Gaussian Error Linear Units (GELUs).https://arxiv.org/abs/1606.08415 [4]Reference 4Searching for Activation Functionshttps://arxiv.org/abs/1710.05941} SiLU is smooth, non-monotonic, and self-gated, with a minimum near $-0.2785$ at $z \approx -1.28$. SwiGLU uses two separate projections, multiplying $\operatorname{SiLU}(xW)$ elementwise by $xV$, then applying an output projection.^{[5]Reference 5GLU Variants Improve Transformerhttps://arxiv.org/abs/2002.05202} Activation choice varies across models: the published Gemma 4 26B-A4B configuration uses a GELU approximation, so check a checkpoint's published configuration before assuming it uses SwiGLU.^{[6]Reference 6Gemma 4 26B-A4B IT Configurationhttps://huggingface.co/google/gemma-4-26B-A4B-it/blob/main/config.json}

View larger

Evaluate all five activations across sample pre-activations to see how their numerical behaviors diverge:

from math import erf, exp, sqrt, tanh

def sigmoid(z: float) -> float:
    return 1.0 / (1.0 + exp(-z)) if z >= 0 else exp(z) / (1.0 + exp(z))

def gelu(z: float) -> float:
    return z * 0.5 * (1.0 + erf(z / sqrt(2.0)))

def swish(z: float) -> float:
    return z * sigmoid(z)

zs = [-2.0, -1.0, 0.0, 1.0, 2.0]
print(f"{'z':>5}  {'sigmoid':>7}  {'tanh':>7}  {'relu':>5}  {'gelu':>6}  {'swish':>6}")
for z in zs:
    print(f"{z:5.1f}  {sigmoid(z):7.3f}  {tanh(z):7.3f}  {max(0.0, z):5.3f}  {gelu(z):6.3f}  {swish(z):6.3f}")

z  sigmoid     tanh   relu    gelu   swish
 -2.0    0.119   -0.964  0.000  -0.046  -0.238
 -1.0    0.269   -0.762  0.000  -0.159  -0.269
  0.0    0.500    0.000  0.000   0.000   0.000
  1.0    0.731    0.762  1.000   0.841   0.731
  2.0    0.881    0.964  2.000   1.954   1.762

At $z = -1.0$, ReLU clips to zero, while GELU returns -0.159 and SiLU returns -0.269. These smooth functions don't cut off the entire negative branch. The direction and size of an update still depend on the loss gradient and the activation's local derivative, which can itself be negative or zero. One unit's curve is settled. A layer with several units still has to say which weight meets which feature when a batch of requests arrives together.

Multi-layer dense networks and tensor shape contracts

A dense, or fully connected, layer connects every output unit to every input feature. The shape of each tensor is the contract that keeps those weights aligned.

PyTorch stores layer weights as (n_out, n_in), meaning number of output units by number of input features.^{[7]Reference 7Linearhttps://docs.pytorch.org/docs/stable/generated/torch.nn.Linear.html} Each row of the weight matrix stores the parameters for one output neuron. For a single column vector $x$:

If a layer takes $D_{\text{in}} = 2$ input features and computes $D_{\text{hid}} = 2$ hidden units, $W$ has shape (2, 2), $x$ has shape (2,), and $z$ has shape (2,).

View larger

Trace our complete 2 -> 2 -> 1 network by hand on request $x = [2.0, 5.0]$:

The hidden pre-activation vector is:

Both values are positive, so ReLU leaves them untouched: $h = [0.10, 3.20]$. The second layer has weights $w_2 = [0.4, 0.1]$ and bias $b_2 = -0.5$:

Run this single-request calculation with explicit row loops to confirm every intermediate number:

from math import exp

x = [2.0, 5.0]
w1 = [[0.3, -0.1], [0.1, 0.6]]
w2 = [0.4, 0.1]
b2 = -0.5

z1 = [sum(row[j] * x[j] for j in range(2)) for row in w1]
h = [max(0.0, value) for value in z1]
z2 = sum(w2[i] * h[i] for i in range(2)) + b2
score = 1.0 / (1.0 + exp(-z2))

print("hidden pre-activation:", [round(value, 2) for value in z1])
print("hidden after ReLU:    ", [round(value, 2) for value in h])
print(f"output logit:          {z2:.3f}")
print(f"sigmoid score:         {score:.3f}")
assert abs(z1[0] - 0.1) < 1e-12
assert abs(z1[1] - 3.2) < 1e-12
assert abs(z2 - (-0.14)) < 1e-12

hidden pre-activation: [0.1, 3.2]
hidden after ReLU:     [0.1, 3.2]
output logit:          -0.140
sigmoid score:         0.465

Batched matrix multiplication and tensor shapes

The same calculation works for many requests at once. Stack $B$ requests as rows of $X \in \mathbb{R}^{B \times D_{\text{in}}}$ and one matrix multiplication evaluates them together. Larger operations can improve GPU utilization, depending on shapes, hardware, and precision. The NumPy example below runs on the CPU.

Under the PyTorch weight layout where $W_1$ is stored as (n_out, n_in), we multiply the batch, whose rows are examples, by $W_1^\top$:

Follow the dimensional contract step by step:

1. Input batch $X$ has shape (B, Din) = (B, 2).
2. Weight matrix $W_1$ has shape (Dhid, Din) = (2, 2). Its transpose $W_1^\top$ has shape (Din, Dhid) = (2, 2).
3. Matrix product $X W_1^\top$ has shape (B, Dhid) = (B, 2).
4. Bias vector $b_1$ has shape (Dhid,) = (2,). Broadcasting reuses that vector for every row without constructing a repeated bias matrix.
5. Activation $H = \text{ReLU}(Z_1)$ preserves shape (B, Dhid) = (B, 2).
6. Output layer $Z_2 = H W_2^\top + b_2$ maps (B, Dhid) to (B, Dout) = (B, 1).
7. Sigmoid $\hat{Y} = \sigma(Z_2)$ maps logits to probability scores of shape (B, Dout).

Execute this batched forward pass on four requests simultaneously:

import numpy as np

x = np.array([
    [2.0, 5.0],  # hand-traced request
    [6.0, 8.0],  # high load
    [1.0, 1.0],  # tiny prompt, short queue
    [8.0, 2.0],  # large prompt, short queue
])
w1 = np.array([[0.3, -0.1], [0.1, 0.6]])
b1 = np.array([0.0, 0.0])
w2 = np.array([[0.4, 0.1]])
b2 = np.array([-0.5])

h = np.maximum(0.0, x @ w1.T + b1)
logits = h @ w2.T + b2
scores = 1.0 / (1.0 + np.exp(-logits))

print("input shape: ", x.shape)
print("hidden shape:", h.shape)
print("score shape: ", scores.shape)
print("scores:", np.round(scores[:, 0], 3))
assert scores.shape == (4, 1)
assert np.isclose(scores[0, 0], 0.465, atol=5e-4)

input shape:  (4, 2)
hidden shape: (4, 2)
score shape:  (4, 1)
scores: [0.465 0.608 0.413 0.641]

The first row matches our scalar hand trace of 0.465.

Beware the silent transpose trap: when $D_{\text{in}} == D_{\text{hid}}$, both X @ W1 and X @ W1.T produce valid matrices of shape (B, Dhid). The operation X @ W1 pairs each neuron with the wrong transposed weights:

import numpy as np

x = np.array([[2.0, 5.0]])
w1 = np.array([[0.3, -0.1], [0.1, 0.6]])

right = x @ w1.T
wrong = x @ w1

print("W1 @ x convention:   ", right[0])
print("X @ W1 other layout: ", wrong[0])
print("same numbers:", np.allclose(right, wrong))
assert np.allclose(right[0], [0.1, 3.2])
assert np.allclose(wrong[0], [1.1, 2.8])

W1 @ x convention:    [0.1 3.2]
X @ W1 other layout:  [1.1 2.8]
same numbers: False

The wrong multiply produces [1.1, 2.8] instead of [0.1, 3.2]. Dimensional shape checks alone won't catch this bug when layers are square. Always verify your weight layout convention against a single-sample test vector.

Counting trainable parameters

Each weight matrix entry and bias element is a trainable parameter that updates during training. In a dense layer with $n_{\text{in}}$ inputs and $n_{\text{out}}$ outputs, the parameter count is:

The product counts one connection weight per input-output pair; the final term adds one bias per output unit. In our 2 -> 2 -> 1 network:

• Layer 1 (2 -> 2): $2 \times 2 + 2 = 6$ parameters.
• Layer 2 (2 -> 1): $1 \times 2 + 1 = 3$ parameters.
• Total parameters: $6 + 3 = 9$.

Compare that tiny network to a flattened image classifier processing $28 \times 28 = 784$ pixels:

def dense_parameters(n_in: int, n_out: int) -> int:
    if n_in <= 0 or n_out <= 0:
        raise ValueError("dense dimensions must be positive")
    return n_in * n_out + n_out

tiny = dense_parameters(2, 2) + dense_parameters(2, 1)
image_classifier = (
    dense_parameters(784, 256)
    + dense_parameters(256, 128)
    + dense_parameters(128, 10)
)

print("2 -> 2 -> 1 parameters:", tiny)
print("784 -> 256 -> 128 -> 10 parameters:", f"{image_classifier:,}")
try:
    dense_parameters(0, 2)
except ValueError as error:
    print(error)

2 -> 2 -> 1 parameters: 9
784 -> 256 -> 128 -> 10 parameters: 235,146
dense dimensions must be positive

The image classifier needs 235,146 parameters because a dense layer connects every pixel to every hidden unit, ignoring spatial adjacency. The count says how large the fit is. It doesn't say which of the nine weights in the tiny network should move when a score is wrong.

Deriving backpropagation for a two-layer network

The forward pass gives us predictions. Training asks how each weight and bias contributed to the loss, so we calculate the partial derivatives with respect to all four parameter arrays: $\frac{\partial \mathcal{L}}{\partial W_1}$, $\frac{\partial \mathcal{L}}{\partial b_1}$, $\frac{\partial \mathcal{L}}{\partial W_2}$, and $\frac{\partial \mathcal{L}}{\partial b_2}$.

Backpropagation applies the multivariate chain rule recursively, moving backward from the final loss through each intermediate activation.^{[8]Reference 8Learning Representations by Back-Propagating Errors.https://doi.org/10.1038/323533a0}

View larger

Here is the step-by-step backward derivation for our 2 -> 2 -> 1 network with Binary Cross-Entropy loss.

Step 1: The loss function

For a batch of $B$ requests with true labels $y \in \{0, 1\}^B$ and sigmoid predictions $\hat{y} = \sigma(z_2)$:

Step 2: Output logit gradient

First take $B = 1$, so the loss contains one sample. Differentiate it with respect to the output logit $z_2$:

Differentiating the cross-entropy loss with respect to $\hat{y}$ gives:

Differentiating sigmoid with respect to its logit $z_2$ gives:

Multiplying them together cancels the shared factor:

The sigmoid derivative's factor $\hat{y}(1 - \hat{y})$ cancels the reciprocal factor in the loss derivative. For a batch of size $B$:

Here $D_{\text{out}}=1$, the loss is unweighted, and the mean covers all $B$ scalar predictions. If you extend this to $C$ independent binary outputs with the default BCEWithLogitsLoss(reduction="mean"), the divisor becomes $BC$: PyTorch averages across the batch and output axes.^{[9]Reference 9BCEWithLogitsLosshttps://docs.pytorch.org/docs/2.14/generated/torch.nn.BCEWithLogitsLoss.html} Using (y_hat - y) / B in that extension makes the gradients $C$ times too large. Summing each request's output losses and then averaging requests defines a different objective. For pos_weight, use the target-weighted gradient formula.

Step 3: Layer 2 parameter gradients

Since $Z_2 = H W_2^\top + b_2$:

Step 4: Propagating gradients to the hidden layer

By the chain rule, the gradient flowing back to hidden activations $H$ is:

Step 5: Passing through the ReLU gate

Because $H = \max(0, Z_1)$, the local derivative is $1$ if $Z_1 > 0$ and $0$ if $Z_1 < 0$. At zero the ordinary derivative doesn't exist; we choose zero for the backward calculation, matching PyTorch ReLU:

Here $\odot$ denotes elementwise multiplication.

Step 6: Layer 1 parameter gradients

Since $Z_1 = X W_1^\top + b_1$:

Verifying NumPy gradients against PyTorch autograd

Implement the backward pass in NumPy, check both weight matrices with finite differences, and, when PyTorch is installed, compare all four parameter gradients with autograd (loss.backward()). Compute loss directly from logits as $\log(1 + e^z) - yz$, using logaddexp to avoid overflowing the exponential. BCEWithLogitsLoss also computes the equivalent loss from logits for numerical stability.^{[9]Reference 9BCEWithLogitsLosshttps://docs.pytorch.org/docs/2.14/generated/torch.nn.BCEWithLogitsLoss.html}

import numpy as np

# Input batch: B=4 requests, Din=2 features
X = np.array([
    [2.0, 5.0],
    [6.0, 8.0],
    [1.0, 1.0],
    [8.0, 2.0],
])
# Ground-truth binary targets (1 = timeout, 0 = success)
y = np.array([[0.0], [1.0], [0.0], [1.0]])

# Network parameters: 2 -> 2 -> 1
W1 = np.array([[0.3, -0.1], [0.1, 0.6]])
b1 = np.array([0.0, 0.0])
W2 = np.array([[0.4, 0.1]])
b2 = np.array([-0.5])

# Forward pass in NumPy
Z1 = X @ W1.T + b1
H = np.maximum(0.0, Z1)
Z2 = H @ W2.T + b2
y_hat = 1.0 / (1.0 + np.exp(-Z2))
loss = np.mean(np.logaddexp(0.0, Z2) - y * Z2)

# Backward pass via chain rule
B = X.shape[0]
dZ2 = (y_hat - y) / B
dW2 = dZ2.T @ H
db2 = np.sum(dZ2, axis=0)

dH = dZ2 @ W2
dZ1 = dH * (Z1 > 0.0)
dW1 = dZ1.T @ X
db1 = np.sum(dZ1, axis=0)

# Numerical gradient check using two-sided finite differences
def compute_loss(w1_val, b1_val, w2_val, b2_val):
    z1_val = X @ w1_val.T + b1_val
    h_val = np.maximum(0.0, z1_val)
    z2_val = h_val @ w2_val.T + b2_val
    return np.mean(np.logaddexp(0.0, z2_val) - y * z2_val)

eps = 1e-6
num_dW2 = np.zeros_like(W2)
for i in range(W2.shape[0]):
    for j in range(W2.shape[1]):
        wp, wm = W2.copy(), W2.copy()
        wp[i, j] += eps
        wm[i, j] -= eps
        num_dW2[i, j] = (compute_loss(W1, b1, wp, b2) - compute_loss(W1, b1, wm, b2)) / (2 * eps)

num_dW1 = np.zeros_like(W1)
for i in range(W1.shape[0]):
    for j in range(W1.shape[1]):
        wp, wm = W1.copy(), W1.copy()
        wp[i, j] += eps
        wm[i, j] -= eps
        num_dW1[i, j] = (compute_loss(wp, b1, W2, b2) - compute_loss(wm, b1, W2, b2)) / (2 * eps)

assert np.allclose(dW2, num_dW2, rtol=0.0, atol=1e-7)
assert np.allclose(dW1, num_dW1, rtol=0.0, atol=1e-7)

# PyTorch autograd verification if installed
torch_checked = False
try:
    import torch
    t_W1 = torch.tensor(W1, dtype=torch.float64, requires_grad=True)
    t_b1 = torch.tensor(b1, dtype=torch.float64, requires_grad=True)
    t_W2 = torch.tensor(W2, dtype=torch.float64, requires_grad=True)
    t_b2 = torch.tensor(b2, dtype=torch.float64, requires_grad=True)
    t_loss = torch.nn.functional.binary_cross_entropy_with_logits(
        torch.relu(torch.nn.functional.linear(torch.tensor(X, dtype=torch.float64), t_W1, t_b1)) @ t_W2.T + t_b2,
        torch.tensor(y, dtype=torch.float64),
    )
    t_loss.backward()
    for manual, parameter in [(dW1, t_W1), (db1, t_b1), (dW2, t_W2), (db2, t_b2)]:
        assert np.allclose(manual, parameter.grad.numpy(), rtol=1e-12, atol=1e-12)
    torch_checked = True
except ImportError:
    pass

print(f"NumPy BCE loss: {loss:.4f}")
print("dW2 gradient:", np.round(dW2, 4))
print("db2 gradient:", np.round(db2, 4))
print("dW1 gradient:", np.round(dW1, 4))
print("db1 gradient:", np.round(db1, 4))
print("Gradient check passed:", True)
print("PyTorch autograd checked:", torch_checked)

# Use the gradients for one optimization step on this same batch.
learning_rate = 0.1
next_loss = compute_loss(
    W1 - learning_rate * dW1, b1 - learning_rate * db1,
    W2 - learning_rate * dW2, b2 - learning_rate * db2,
)
print(f"BCE after one gradient step: {next_loss:.4f}")
assert next_loss < loss

NumPy BCE loss: 0.5252
dW2 gradient: [[-0.2631 -0.2639]]
db2 gradient: [0.0319]
dW1 gradient: [[-0.3878 -0.1113]
 [-0.097  -0.0278]]
db1 gradient: [0.0128 0.0032]
Gradient check passed: True
PyTorch autograd checked: True
BCE after one gradient step: 0.4980

The shown run has PyTorch installed, so both weight and both bias gradients pass the tight comparison. With NumPy alone, the finite-difference checks still run, but PyTorch autograd checked prints False.

The checks also depend on where we evaluate them. Every hidden pre-activation here is positive, and the perturbations of size 1e-6 keep them positive. At a ReLU corner, a central difference can disagree with the chosen backward value because there is no ordinary derivative there.^{[10]Reference 10torch.autograd.gradcheck.gradcheckhttps://docs.pytorch.org/docs/main/generated/torch.autograd.gradcheck.gradcheck.html} Compute the same check directly at zero and on either side:

import torch

eps = 1e-6
for value in [-0.1, 0.0, 0.1]:
    z = torch.tensor(value, dtype=torch.float64, requires_grad=True)
    torch.relu(z).backward()
    central = (max(0.0, value + eps) - max(0.0, value - eps)) / (2 * eps)
    print(f"z={value:4.1f}: backward={z.grad.item():.1f}, central={central:.1f}")

z=-0.1: backward=0.0, central=0.0
z= 0.0: backward=0.0, central=0.5
z= 0.1: backward=1.0, central=1.0

At ReLU's nondifferentiable corner, the central difference gives 0.5 while PyTorch selects its backward value 0.0; this is expected. Check away from corners in double precision, keep perturbations within the same active branches, and compare several perturbation sizes: very small differences can lose precision, while large ones can cross branches. The smooth-region agreement validates this fixture on the tested inputs and weights; evaluate other configurations separately.

The last step uses those gradients: loss on the same four requests falls from 0.5252 to 0.4980. This verifies the update on this batch at this learning rate. Evaluate held-out requests and other step sizes to assess broader behavior.

Our backward calculation reuses $H$ for $\text{d}W_2$, the signs of $Z_1$ for the ReLU gate, and $X$ for $\text{d}W_1$. These values can be saved or recomputed when needed. Activation memory can be large at long sequence lengths or large batches, but weights, gradients, and optimizer state also contribute to training memory. Which term dominates depends on the workload. The formulas can still be right while the run breaks. What happens when one feature is a token count and another is a zero-to-ten score?

Diagnosing feature scale and numerical overflow

Shapes and gradient formulas can all be correct while scale still breaks a run. Two common causes are features measured in very different units, and exponentials that overflow.

Feature scale mismatch

Suppose CodeAssist changes units: queue wait is measured on a 0 to 10 scale, while prompt length is measured in raw tokens ranging from 10 to 1,000:

A queue wait of 8 contributes $0.3(8) = 2.4$. A prompt length of 900 tokens contributes $0.3(900) = 270.0$. Even though both features share the exact same coefficient of 0.3, the token count dominates the affine sum by two orders of magnitude:

raw = [8.0, 900.0]  # queue-wait score, prompt token count
raw_contributions = [0.3 * value for value in raw]

mean = [5.0, 500.0]
scale = [2.0, 200.0]
standardized = [(value - m) / s for value, m, s in zip(raw, mean, scale)]
scaled_contributions = [0.3 * value for value in standardized]

print("raw contributions:       ", raw_contributions)
print("standardized features:    ", [round(value, 2) for value in standardized])
print("scaled contributions:     ", [round(value, 2) for value in scaled_contributions])
assert raw_contributions == [2.4, 270.0]

raw contributions:        [2.4, 270.0]
standardized features:     [1.5, 2.0]
scaled contributions:      [0.45, 0.6]

Standardizing each feature by subtracting training-set means and dividing by training-set standard deviations brings their typical scales closer: 1.5 and 2.0 here. Reusing the illustrative coefficients gives contributions 0.45 and 0.60; a real model's coefficients must be trained or transformed for the chosen units. Standardization can help optimization, while feature correlations and outliers still need separate attention.

Fit the means and scales on the training split, then reuse them for validation, test data, and inference. Using test data to fit the scaler leaks information. A constant training feature needs a policy too: scikit-learn's StandardScaler uses a scale factor of one.^{[11]Reference 11StandardScalerhttps://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html}

Exponential overflow in activations

The standard sigmoid formula 1.0 / (1.0 + exp(-z)) is numerically safe for $z = 2.6$, but it crashes for extreme negative logits. At $z = -1000$, exp(-z) evaluates exp(1000), exceeding the floating-point range of 64-bit floats and triggering an OverflowError.

Use an algebraically equivalent branch that evaluates $e^z$ for negative inputs:

from math import exp

def stable_sigmoid(z: float) -> float:
    if z >= 0:
        return 1.0 / (1.0 + exp(-z))
    ez = exp(z)
    return ez / (1.0 + ez)

for logit in [-1000.0, -2.0, 0.0, 2.0, 1000.0]:
    print(f"logit={logit:7.1f} score={stable_sigmoid(logit):.6f}")

logit=-1000.0 score=0.000000
logit=   -2.0 score=0.119203
logit=    0.0 score=0.500000
logit=    2.0 score=0.880797
logit= 1000.0 score=1.000000

For finite negative logits, $e^z$ produces small numbers. Extreme outputs can still round to zero or one, making log(score) unsafe; compute the loss directly from logits.

Compare two extremely confident wrong predictions using the loss directly from logits:

import numpy as np

logits = np.array([-1000.0, 1000.0])
targets = np.array([1.0, 0.0])
losses = np.logaddexp(0.0, logits) - targets * logits
print("per-request losses:", losses)
print("mean BCE loss:", losses.mean())
assert np.all(np.isfinite(losses))
assert np.allclose(losses, [1000.0, 1000.0])

per-request losses: [1000. 1000.]
mean BCE loss: 1000.0

Now combine shape assertions, the stable sigmoid branch, and an illustrative routing rule in a checked toy scoring report:

import numpy as np

def stable_sigmoid(values: np.ndarray) -> np.ndarray:
    scores = np.empty_like(values, dtype=float)
    nonnegative = values >= 0
    scores[nonnegative] = 1.0 / (1.0 + np.exp(-values[nonnegative]))
    exp_values = np.exp(values[~nonnegative])
    scores[~nonnegative] = exp_values / (1.0 + exp_values)
    return scores

feature_names = ["prompt-size score", "queue-wait score"]
request_ids = ["CA-104", "CA-208", "CA-311", "CA-412"]
x = np.array([
    [2.0, 5.0],
    [6.0, 8.0],
    [1.0, 1.0],
    [8.0, 2.0],
])
w1 = np.array([[0.3, -0.1], [0.1, 0.6]])
b1 = np.array([0.0, 0.0])
w2 = np.array([[0.4, 0.1]])
b2 = np.array([-0.5])

assert x.shape[1] == len(feature_names) == w1.shape[1]
assert w1.shape[0] == b1.shape[0] == w2.shape[1]
assert w2.shape[0] == b2.shape[0] == 1

hidden = np.maximum(0.0, x @ w1.T + b1)
logits = hidden @ w2.T + b2
scores = stable_sigmoid(logits)

assert scores.shape == (len(request_ids), 1)
assert np.isclose(scores[0, 0], 0.465, atol=5e-4)
assert np.all(np.isfinite(stable_sigmoid(np.array([-1000.0, 1000.0]))))

print("Timeout-risk scoring report")
for request_id, score in zip(request_ids, scores[:, 0], strict=True):
    route = "slow_path" if score >= 0.60 else "fast_path"
    print(f"{request_id}: score={score:.3f} -> {route}")

Timeout-risk scoring report
CA-104: score=0.465 -> fast_path
CA-208: score=0.608 -> slow_path
CA-311: score=0.413 -> fast_path
CA-412: score=0.641 -> slow_path

Requests CA-208 and CA-412 cross the illustrative 0.60 threshold and receive slow_path labels. Hand-picked weights produce these four scores, with no training or calibration. What functions can a wide enough network represent, and what other conditions shape whether training learns one?

Universal approximation and its practical limits

The four-request report used hand-picked weights. Universal approximation theorems ask what functions a wide enough network can represent, and they come with important limits. Cybenko (1989) proved a result for sigmoidal activations, and Hornik et al. (1989) established related results.^{[12]Reference 12Approximation by Superpositions of a Sigmoidal Function.https://doi.org/10.1007/BF02551274 [13]Reference 13Multilayer Feedforward Networks are Universal Approximators.https://doi.org/10.1016/0893-6080(89)90020-8} A useful version allows a continuous non-polynomial activation, such as sigmoid or ReLU: with biases, a linear output, and enough hidden units, one hidden layer can approximate any continuous function uniformly on a compact (closed and bounded) domain to arbitrary accuracy $\epsilon > 0$.^{[14]Reference 14Multilayer Feedforward Networks with a Non-Polynomial Activation Function Can Approximate Any Functionhttps://archive.nyu.edu/handle/2451/14329} Non-linearity alone isn't sufficient; a fixed polynomial activation is a counterexample.

For one-dimensional intuition, take $a > 0$. The difference $\sigma(x + a) - \sigma(x - a)$ forms a smooth bump. Three shifted ReLU units form an exact triangle:

It rises from zero to one on $[0, 1]$, falls to zero on $[1, 2]$, and is zero outside that interval. Shifted and scaled triangles can approximate a continuous curve by increasingly fine piecewise-linear segments. This example builds one-dimensional intuition; higher-dimensional approximation requires additional arguments.

To use this result in practice, ask three further questions:

1. How much width is enough? Required width depends on the target function, activation, input dimension, and desired error. Depth-separation results construct particular function families that small deep networks represent much more efficiently than shallow ones.^{[15]Reference 15Benefits of Depth in Neural Networkshttps://arxiv.org/abs/1602.04485} These examples show where depth helps; gains for a particular task depend on its function and architecture.

2. Representability and learnability: The theorem establishes that a suitable set of weights exists in parameter space. Training still depends on whether gradient descent can reach suitable weights from its initialization; optimization can encounter poor local minima or flat saddle points.

3. Can training data identify the function? The mathematical approximation covers the chosen compact domain, including unsampled points. The theorem assumes a target function and proves suitable weights exist; finite training examples may not identify those weights or yield good test performance. Architectural assumptions can help: CNNs share local filters across an image, while attention, positions, and masks determine how a Transformer combines tokens. These choices need to fit the task. A CNN's translation behavior depends on its architecture and boundary handling; position-aware or causally masked attention can break permutation equivariance.

Practice problems

1. In bias-shifts-threshold.py, change the bias from -4.0 to -6.0. Calculate the new logits and explain why both requests receive lower risk scores even though their input features didn't change.
2. Extend batch-forward-pass.py to support three input features (add a server CPU load score) while keeping two hidden units. Write out the required shapes of $X$, $W_1$, and $b_1$ before running the code.
3. In backpropagation-two-layer-mlp.py, change the ground-truth target vector from [0, 1, 0, 1] to all ones [1, 1, 1, 1]. Predict the sign of $\text{d}Z_2$ and explain which direction the weights must move to reduce the loss.

Mastery check