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?

The vectors and matrices lesson represented incoming observations as feature vectors and combined them with dot products. A neural network builds directly on that linear foundation. Its forward pass evaluates an output from inputs using stored weight parameters and biases. We'll hold those weights fixed first so you can inspect every floating-point addition and multiplication by hand before training adjusts 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 an illustrative zero-to-ten scale rather than raw token counts or seconds. 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$. That would mean a request with zero prompt size and zero queue wait must always produce a score of zero. The bias shifts the hyperplane away from the origin so the model can set its own baseline threshold.

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

$$z = w^\top x + b = 0.7(6) + 0.3(8) - 4 = 2.6$$

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:

$$\sigma(z) = \frac{1}{1 + e^{-z}}$$

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. Sigmoid outputs stay strictly between zero and one for any finite input.

A bounded score isn't proof of a 93.1% timeout probability. These parameters were chosen by hand for illustration. A true probability estimate requires a calibrated statistical model trained on verified request outcomes.^{[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:

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:

At a score threshold of 0.5, the high-load request crosses the alert line while the quiet request stays safely 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.

Why shouldn't you treat 0.931 as a verified timeout probability in production?

Answer

These parameters were chosen by hand for demonstration, not learned from historical data. Sigmoid bounds the output between zero and one, but calibration requires empirical proof that requests scoring near 0.93 actually time out roughly 93% of the time.

Linear maps collapse without non-linearities

One neuron scores a weighted sum. Why can't we build a deep network simply by stacking multiple linear layers?

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$:

$$h = W_1 x + b_1$$ $$y = W_2 h + b_2$$

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

$$y = W_2(W_1 x + b_1) + b_2 = (W_2 W_1)x + (W_2 b_1 + b_2)$$

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

$$y = W_{\text{merged}} x + b_{\text{merged}}$$

No curvature appeared. Because the product of two matrices is just another matrix, stacking fifty linear layers produces the exact same representational power as a single linear layer. A fifty-layer purely linear network can't even solve the classic XOR classification problem, where $(0, 1)$ and $(1, 0)$ are positive but $(0, 0)$ and $(1, 1)$ are negative. No single flat hyperplane can separate 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:

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]. Because of that clipping, the network output changes from -0.290 to -0.250. The two layers can no longer be compressed into a single matrix multiplication.

Comparing modern activation functions

Choosing the right non-linearity shapes how gradients flow during optimization. Five core functions appear throughout modern deep learning architectures:

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 works well for binary classification outputs, but it causes severe vanishing gradients in deep hidden layers because $|z| > 3$ pushes derivatives near zero. Its outputs are also strictly positive (not zero-centered), which forces all incoming weight updates in a layer to share the same sign.

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 zero-centered, meaning negative inputs produce negative outputs. That property centers hidden activations around zero and improves gradient descent dynamics compared to sigmoid. It still saturates for $|z| > 2$, though, so deep stacks still suffer from 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 eliminates vanishing gradients in the positive regime and evaluates with minimal CPU and GPU overhead. Its weakness is the dying ReLU problem: if a large negative bias or massive gradient update pushes pre-activations negative for all inputs, the unit outputs zero and receives zero gradient forever.

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 abrupt corner at zero, GELU dips slightly negative ($\approx -0.17$ near $z \approx -0.75$) with non-zero curvature everywhere. That smooth tail keeps small gradient signals alive where ReLU would cut them off completely. GELU is the standard activation across BERT, GPT-2, GPT-3, and Vision Transformers.

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)$. Discovered through neural architecture search, Swish is smooth, non-monotonic, and self-gated. It reaches a minimum of roughly $-0.278$ near $z \approx -1.28$. Modern Large Language Models rely heavily on Swish in SwiGLU feed-forward layers (used in LLaMA, Mistral, Gemma, and DeepSeek), where a linear projection is elementwise gated by its own Swish transformation.

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

At $z = -1.0$, ReLU clips to zero, while GELU returns -0.159 and Swish returns -0.269. That smooth dip lets gradients push dead units back toward the active positive zone during training.

Why do modern large language models favor GELU and Swish over standard ReLU in their feed-forward blocks?

Answer

ReLU has a sharp discontinuity at zero and completely zeroes out gradients for negative pre-activations, risking permanently dead units. GELU and Swish are smooth and non-monotonic; their slight negative dips allow small gradient signals to flow through negative inputs during backpropagation.

Multi-layer dense networks and tensor shape contracts

Multiple neurons combine into a dense layer, also known as a fully connected layer. In a dense layer, every output unit connects to every input feature.

PyTorch stores layer weights as (n_out, n_in), meaning number of output units by number of input features.^{[4]Reference 4Linearhttps://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$:

$$z = Wx + b$$

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,).

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

$$W_1 = \begin{bmatrix} 0.3 & -0.1 \\ 0.1 & 0.6 \end{bmatrix}, \quad b_1 = \begin{bmatrix} 0.0 \\ 0.0 \end{bmatrix}$$

The hidden pre-activation vector is:

$$z_1 = W_1 x + b_1 = \begin{bmatrix} 0.3(2.0) - 0.1(5.0) \\ 0.1(2.0) + 0.6(5.0) \end{bmatrix} = \begin{bmatrix} 0.10 \\ 3.20 \end{bmatrix}$$

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$:

$$z_2 = w_2^\top h + b_2 = 0.4(0.10) + 0.1(3.20) - 0.5 = 0.04 + 0.32 - 0.5 = -0.14$$ $$\text{score} = \sigma(-0.14) = \frac{1}{1 + e^{0.14}} \approx 0.465$$

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

Batched matrix multiplication and tensor shapes

Evaluating requests one by one wastes GPU compute and memory bandwidth. Stacking $B$ requests as rows of an input batch matrix $X \in \mathbb{R}^{B \times D_{\text{in}}}$ unlocks parallel matrix operations on modern tensor cores.

Under the PyTorch weight layout where $W_1$ is stored as (n_out, n_in), we multiply the row-major batch by $W_1^\top$:

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

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,). NumPy uses broadcasting to automatically copy $b_1$ across all $B$ rows.
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:

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). Neither NumPy nor PyTorch will throw a shape error. But X @ W1 pairs each neuron with the wrong transposed weights:

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:

$$\text{parameters} = n_{\text{out}} \times n_{\text{in}} + n_{\text{out}}$$

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:

The image classifier needs 235,146 parameters because a dense layer connects every pixel to every hidden unit, ignoring spatial adjacency.

Deriving backpropagation for a two-layer network

Evaluating the forward pass computes predictions. To train the network, we must calculate the partial derivative of the loss with respect to every weight and bias: $\frac{\partial \mathcal{L}}{\partial W_1}$, $\frac{\partial \mathcal{L}}{\partial b_1}$, $\frac{\partial \mathcal{L}}{\partial W_2}$, $\frac{\partial \mathcal{L}}{\partial b_2}$.

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

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)$:

$$\mathcal{L} = -\frac{1}{B} \sum_{i=1}^B \left[ y_i \log \hat{y}_i + (1 - y_i) \log(1 - \hat{y}_i) \right]$$

Step 2: Output logit gradient

Differentiate $\mathcal{L}$ with respect to the output logit $z_2$ for a single sample:

$$\frac{\partial \mathcal{L}}{\partial z_2} = \frac{\partial \mathcal{L}}{\partial \hat{y}} \cdot \frac{\partial \hat{y}}{\partial z_2}$$

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

$$\frac{\partial \mathcal{L}}{\partial \hat{y}} = -\frac{y}{\hat{y}} + \frac{1 - y}{1 - \hat{y}} = \frac{\hat{y} - y}{\hat{y}(1 - \hat{y})}$$

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

$$\frac{\partial \hat{y}}{\partial z_2} = \hat{y}(1 - \hat{y})$$

Multiplying them together produces an elegant cancellation:

$$\frac{\partial \mathcal{L}}{\partial z_2} = \left( \frac{\hat{y} - y}{\hat{y}(1 - \hat{y})} \right) \cdot \hat{y}(1 - \hat{y}) = \hat{y} - y$$

The sigmoid derivative denominator cancels the loss derivative denominator completely. For a batch of size $B$:

$$\frac{\partial \mathcal{L}}{\partial Z_2} = \frac{1}{B} (\hat{Y} - Y) \quad \text{shape: } (B, D_{\text{out}})$$

Step 3: Layer 2 parameter gradients

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

$$\frac{\partial \mathcal{L}}{\partial W_2} = \left(\frac{\partial \mathcal{L}}{\partial Z_2}\right)^\top H \quad \text{shape: } (D_{\text{out}}, D_{\text{hid}})$$ $$\frac{\partial \mathcal{L}}{\partial b_2} = \sum_{i=1}^B \left(\frac{\partial \mathcal{L}}{\partial Z_2}\right)_i \quad \text{shape: } (D_{\text{out}},)$$

Step 4: Propagating gradients to the hidden layer

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

$$\frac{\partial \mathcal{L}}{\partial H} = \frac{\partial \mathcal{L}}{\partial Z_2} W_2 \quad \text{shape: } (B, D_{\text{hid}})$$

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 \le 0$:

$$\frac{\partial \mathcal{L}}{\partial Z_1} = \frac{\partial \mathcal{L}}{\partial H} \odot \mathbf{1}_{Z_1 > 0} \quad \text{shape: } (B, D_{\text{hid}})$$

Here $\odot$ denotes elementwise multiplication.

Step 6: Layer 1 parameter gradients

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

$$\frac{\partial \mathcal{L}}{\partial W_1} = \left(\frac{\partial \mathcal{L}}{\partial Z_1}\right)^\top X \quad \text{shape: } (D_{\text{hid}}, D_{\text{in}})$$ $$\frac{\partial \mathcal{L}}{\partial b_1} = \sum_{i=1}^B \left(\frac{\partial \mathcal{L}}{\partial Z_1}\right)_i \quad \text{shape: } (D_{\text{hid}},)$$

Verifying NumPy gradients against PyTorch autograd

Implement this exact backward pass in pure NumPy and verify every single gradient against PyTorch autograd (loss.backward()):

Every hand-derived gradient matches PyTorch autograd to 12 decimal places. Notice that backpropagation requires saving forward activations ($X$, $Z_1$, and $H$) in memory: we needed $H$ to compute $\text{d}W_2$, $Z_1$ to compute the ReLU gate, and $X$ to compute $\text{d}W_1$. In deep networks, storing these forward activations is what dominates GPU training memory.

Why does backpropagation require holding forward activations in GPU memory until the backward pass completes?

Answer

Parameter gradients depend directly on incoming activations: $\frac{\partial \mathcal{L}}{\partial W_2} = (\frac{\partial \mathcal{L}}{\partial Z_2})^\top H$ and $\frac{\partial \mathcal{L}}{\partial W_1} = (\frac{\partial \mathcal{L}}{\partial Z_1})^\top X$. Without caching $H$, $X$, and pre-activation signs for ReLU, the backward pass cannot evaluate weight updates without recomputing the forward pass from scratch.

Diagnosing feature scale and numerical overflow

Even with correct matrix math and backpropagation equations, real implementations fail if numerical scale is ignored.

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:

$$z = 0.3 \cdot \text{queue wait} + 0.3 \cdot \text{prompt tokens} - 2.0$$

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:

Standardizing each feature by subtracting training-set means and dividing by training-set standard deviations brings inputs onto comparable scales: 1.5 and 2.0. Their contributions become 0.45 and 0.60. Standardization prevents wide-ranging features from dwarfing small features and warping the loss surface into elongated ravines.

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:

$$\sigma(z) = \begin{cases} \frac{1}{1 + e^{-z}} & \text{if } z \ge 0 \\ \frac{e^z}{1 + e^z} & \text{if } z < 0 \end{cases}$$

For negative logits, evaluating $e^z$ produces numbers close to zero rather than astronomical numbers, ensuring numerical stability across the entire real line.

Now combine shape assertions, the stable sigmoid branch, and routing logic into a verified scoring pipeline:

Requests CA-208 and CA-412 cross the 0.60 routing threshold, redirecting them to isolated workers with dedicated timeouts.

Universal approximation and its practical limits

The Universal Approximation Theorem, established by Cybenko (1989) and extended by Hornik et al. (1989), shows that a feed-forward network with a single hidden layer, a sufficient finite number of non-linear neurons, and a linear output can approximate any continuous function on a compact (closed and bounded) subset of $\mathbb{R}^n$ to arbitrary precision $\epsilon > 0$.^{[6]Reference 6Approximation by Superpositions of a Sigmoidal Function.https://doi.org/10.1007/BF02551274 [7]Reference 7Multilayer Feedforward Networks are Universal Approximators.https://doi.org/10.1016/0893-6080(89)90020-8}

The geometric intuition is surprisingly visual:

• Take two opposing sigmoid functions: $\sigma(x + a) - \sigma(x - a)$. Their difference forms a localized bump or ridge in the coordinate space.
• With ReLU, two shifted ramp functions can construct a triangle bump: $\text{ReLU}(x) - 2\text{ReLU}(x - 1) + \text{ReLU}(x - 2)$.
• By summing many of these localized bumps across dimensions, a wide hidden layer can tile any continuous surface, much like a Riemann sum tiles an integral or Lego blocks build a curved arch.

While theoretically reassuring, the theorem has three critical practical catches:

1. Exponential width explosion: Approximating complex or highly oscillatory functions with just a single hidden layer can demand exponentially many hidden neurons ($O(2^d)$ for $d$ dimensions). Deep architectures compose functions hierarchically: each added layer folds the input space, creating exponential expressive power with only linear parameter growth.

2. Existence doesn't mean learnability: The theorem guarantees that a set of weights exists in parameter space. It doesn't guarantee that gradient descent will find those weights from random initialization, or that optimization won't get trapped in poor local minima or flat saddle points.

3. Generalization requires inductive bias: The theorem only guarantees function fitting on the compact region covered by training data; it says nothing about unseen test points. Architectures like CNNs and Transformers succeed because their structure encodes domain symmetries: CNNs enforce translation invariance on image grids, and Transformers enforce permutation equivariance across token sequences.

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.

Compare your solutions with these worked answers.

Answer

1. Lowering the bias by 2 subtracts exactly 2.0 from every logit: quiet becomes -3.70 (score $\approx 0.024$), and high load becomes 0.60 (score $\approx 0.646$). Because sigmoid increases monotonically with its input, reducing the logit always reduces the output score.
2. With three features and two hidden units: $X$ has shape (4, 3), $W_1$ has shape (2, 3) under PyTorch layout, and $b_1$ has shape (2,). Hidden activations retain shape (4, 2), so $W_2$ (shape (1, 2)) and output logits (shape (4, 1)) remain unchanged.
3. When targets are all $1.0$, $y_i = 1$ for all samples. Since predictions $\hat{y}_i < 1$, the logit error $\text{d}Z_2 = \frac{1}{B}(\hat{Y} - Y)$ is strictly negative for every sample. Negative gradients mean that gradient descent updates ($W \leftarrow W - \alpha \nabla W$) will add positive increments, pulling logits upward to increase sigmoid scores toward 1.0.

Mastery check