Neural Networks & Backpropagation Calculus | AI & Machine Learning Engineering

Neural networks extend linear models by interweaving linear matrix transformations with non-linear activation functions. The Universal Approximation Theorem guarantees that a feedforward network with a single hidden layer and non-linear activations can approximate any continuous function on compact subsets of $\mathbb{R}^d$ to arbitrary precision.

To optimize millions of parameters simultaneously, deep learning relies on Backpropagation, an elegant application of the multivariate chain rule across directed computation graphs.

Weight GradientMatrix Calculus

delta * a^T

The outer product of the incoming error vector and outgoing activation yields the exact weight Jacobian.

1. Activation Functions and Exact Derivatives#

Non-linear activations prevent deep networks from collapsing into a single trivial linear transformation $W_2(W_1 x) = W_{\text{eff}} x$ .

Architecture Diagram

Rendering diagram...

Analytical Derivatives#

Sigmoid:

\sigma(z) = \frac{1}{1 + e^{-z}} \implies \frac{d\sigma(z)}{dz} = \sigma(z)(1 - \sigma(z))

Hyperbolic Tangent (Tanh):

\tanh(z) = \frac{e^z - e^{-z}}{e^z + e^{-z}} \implies \frac{d\tanh(z)}{dz} = 1 - \tanh^2(z)

Rectified Linear Unit (ReLU):

\text{ReLU}(z) = \max(0, z) \implies \frac{d\text{ReLU}(z)}{dz} = \begin{cases} 1 & \text{if } z > 0 \\ 0 & \text{if } z < 0 \end{cases}

Leaky ReLU ( $\alpha \approx 0.01$ ):

\text{LReLU}(z) = \max(\alpha z, z) \implies \frac{d\text{LReLU}(z)}{dz} = \begin{cases} 1 & \text{if } z > 0 \\ \alpha & \text{if } z < 0 \end{cases}

The Vanishing Gradient Trap in Deep Networks

For Sigmoid, the maximum derivative occurs at

z = 0

, where

\sigma'(0) = 0.25

. In a 10-layer network, backpropagating through 10 sigmoid layers scales the gradient by

(0.25)^{10} \approx 9.5 \times 10^{-7}

, causing early layer weights to freeze. ReLU circumvents this failure mode by maintaining a constant derivative of

1.0

for all positive inputs.

2. Step-by-Step Matrix Calculus Derivation of Backpropagation#

Consider an $L$ -layer network. For layer $l \in \{1, \dots, L\}$ :

$W^{(l)} \in \mathbb{R}^{n_l \times n_{l-1}}$ : Weight matrix
$b^{(l)} \in \mathbb{R}^{n_l}$ : Bias vector
$z^{(l)} = W^{(l)} a^{(l-1)} + b^{(l)}$ : Linear pre-activation vector
$a^{(l)} = \sigma(z^{(l)})$ : Non-linear activation vector (with $a^{(0)} = x$ )

Let $J$ denote the scalar loss function.

The Error Vector Definition#

We define the error vector $\delta^{(l)}$ of layer $l$ as the partial derivative of scalar loss $J$ with respect to the linear pre-activation $z^{(l)}$ :

\delta^{(l)} \equiv \frac{\partial J}{\partial z^{(l)}} \in \mathbb{R}^{n_l}

Step 1: Base Case (Output Layer Error $\delta^{(L)}$)

By the chain rule:

\delta_i^{(L)} = \frac{\partial J}{\partial z_i^{(L)}} = \sum_k \frac{\partial J}{\partial a_k^{(L)}} \frac{\partial a_k^{(L)}}{\partial z_i^{(L)}}

Because element-wise activations depend only on their own coordinate ( $a_i^{(L)} = \sigma(z_i^{(L)})$ ):

\frac{\partial a_k^{(L)}}{\partial z_i^{(L)}} = \begin{cases} \sigma'(z_i^{(L)}) & \text{if } k = i \\ 0 & \text{if } k \ne i \end{cases}

Therefore, in vector notation:

\delta^{(L)} = \nabla_{a^{(L)}} J \odot \sigma'(z^{(L)})

Where $\odot$ represents the Hadamard element-wise vector product.

Step 2: Inductive Step (Hidden Layer Error Propagation)

Now, relate the error at hidden layer $l$ to the error at the subsequent layer $l+1$ :

\delta_j^{(l)} = \frac{\partial J}{\partial z_j^{(l)}} = \sum_{k=1}^{n_{l+1}} \frac{\partial J}{\partial z_k^{(l+1)}} \frac{\partial z_k^{(l+1)}}{\partial z_j^{(l)}}

By definition, $\frac{\partial J}{\partial z_k^{(l+1)}} = \delta_k^{(l+1)}$ . Next, evaluate $\frac{\partial z_k^{(l+1)}}{\partial z_j^{(l)}}$ :

z_k^{(l+1)} = \sum_{m=1}^{n_l} W_{k,m}^{(l+1)} a_m^{(l)} + b_k^{(l+1)} = \sum_{m=1}^{n_l} W_{k,m}^{(l+1)} \sigma(z_m^{(l)}) + b_k^{(l+1)}

Differentiating with respect to $z_j^{(l)}$ :

\frac{\partial z_k^{(l+1)}}{\partial z_j^{(l)}} = W_{k,j}^{(l+1)} \sigma'(z_j^{(l)})

Substitute this back into the summation:

\delta_j^{(l)} = \sum_{k=1}^{n_{l+1}} \delta_k^{(l+1)} W_{k,j}^{(l+1)} \sigma'(z_j^{(l)}) = \left( \sum_{k=1}^{n_{l+1}} W_{k,j}^{(l+1)} \delta_k^{(l+1)} \right) \sigma'(z_j^{(l)})

In matrix-vector notation, this is the canonical Backpropagation Equation:

\delta^{(l)} = \Big( (W^{(l+1)})^T \delta^{(l+1)} \Big) \odot \sigma'(z^{(l)})

Step 3: Compute Gradients for Weights and Biases

Now that $\delta^{(l)}$ is known, we compute the gradient with respect to weight matrix $W^{(l)}$ :

\frac{\partial J}{\partial W_{j,k}^{(l)}} = \frac{\partial J}{\partial z_j^{(l)}} \frac{\partial z_j^{(l)}}{\partial W_{j,k}^{(l)}} = \delta_j^{(l)} a_k^{(l-1)}

In outer product matrix notation:

\frac{\partial J}{\partial W^{(l)}} = \delta^{(l)} (a^{(l-1)})^T

For the bias vector:

\frac{\partial J}{\partial b^{(l)}} = \delta^{(l)}

3. Adam Optimization: Bias-Correction Derivation#

The Adam (Adaptive Moment Estimation) optimizer maintains exponentially decaying running averages of past gradients ( $m_t$ ) and past squared gradients ( $v_t$ ):

m_t = \beta_1 m_{t-1} + (1 - \beta_1) g_t

v_t = \beta_2 v_{t-1} + (1 - \beta_2) g_t^2

Mathematical Proof of Moment Bias Correction#

Because moments are initialized to zero ( $m_0 = 0$ ), the running average is heavily biased toward zero in initial iterations.

Expanding $m_t$ recursively:

m_t = (1 - \beta_1) \sum_{i=1}^t \beta_1^{t-i} g_i

Taking expectations on both sides, assuming gradients $g_i$ come from a stationary distribution with true mean $\mathbb{E}[g_i] = \mathbb{E}[g]$ :

\mathbb{E}[m_t] = \mathbb{E}\left[ (1 - \beta_1) \sum_{i=1}^t \beta_1^{t-i} g_i \right] = \mathbb{E}[g] (1 - \beta_1) \sum_{i=1}^t \beta_1^{t-i}

The summation is a finite geometric series:

\sum_{i=1}^t \beta_1^{t-i} = \sum_{k=0}^{t-1} \beta_1^k = \frac{1 - \beta_1^t}{1 - \beta_1}

Substituting this sum back:

\mathbb{E}[m_t] = \mathbb{E}[g] (1 - \beta_1) \left( \frac{1 - \beta_1^t}{1 - \beta_1} \right) = \mathbb{E}[g] (1 - \beta_1^t)

To obtain an unbiased estimator $\hat{m}_t$ such that $\mathbb{E}[\hat{m}_t] = \mathbb{E}[g]$ , we divide by $(1 - \beta_1^t)$ :

\hat{m}_t = \frac{m_t}{1 - \beta_1^t}, \qquad \hat{v}_t = \frac{v_t}{1 - \beta_2^t}

The final parameter update is:

\theta_t = \theta_{t-1} - \frac{\eta}{\sqrt{\hat{v}_t} + \epsilon} \hat{m}_t

4. From-Scratch Implementation: Zero Dependencies#

The following implementations build a 2-layer Multi-Layer Perceptron (MLP) trained on the classic non-linear XOR problem from scratch in pure Python and pure Rust.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
import math
import random
from typing import List

class MultiLayerPerceptronFromScratch:
    def __init__(self, input_dim: int, hidden_dim: int, output_dim: int, lr: float = 0.5):
        self.lr = lr
        self.hidden_dim = hidden_dim
        self.output_dim = output_dim

        # Xavier/He uniform weight initialization
        random.seed(42)
        scale1 = 1.0 / math.sqrt(input_dim)
        self.W1 = [[random.uniform(-scale1, scale1) for _ in range(input_dim)] for _ in range(hidden_dim)]
        self.b1 = [0.0] * hidden_dim

        scale2 = 1.0 / math.sqrt(hidden_dim)
        self.W2 = [[random.uniform(-scale2, scale2) for _ in range(hidden_dim)] for _ in range(output_dim)]
        self.b2 = [0.0] * output_dim

    @staticmethod
    def _sigmoid(z: float) -> float:
        if z < -700.0: return 0.0
        if z > 700.0: return 1.0
        return 1.0 / (1.0 + math.exp(-z))

    @staticmethod
    def _sigmoid_derivative(a: float) -> float:
        return a * (1.0 - a)

    def forward(self, x: List[float]):
        # Hidden layer: z1 = W1 * x + b1, a1 = sigmoid(z1)
        z1 = [sum(self.W1[j][k] * x[k] for k in range(len(x))) + self.b1[j] for j in range(self.hidden_dim)]
        a1 = [self._sigmoid(val) for val in z1]

        # Output layer: z2 = W2 * a1 + b2, a2 = sigmoid(z2)
        z2 = [sum(self.W2[j][k] * a1[k] for k in range(self.hidden_dim)) + self.b2[j] for j in range(self.output_dim)]
        a2 = [self._sigmoid(val) for val in z2]

        return a1, a2

    def train_step(self, x: List[float], y: List[float]):
        # 1. Forward propagation
        a1, a2 = self.forward(x)

        # 2. Output layer error: delta2 = (a2 - y) * sigma'(z2)
        delta2 = [(a2[j] - y[j]) * self._sigmoid_derivative(a2[j]) for j in range(self.output_dim)]

        # 3. Hidden layer error: delta1 = (W2^T * delta2) * sigma'(z1)
        delta1 = [0.0] * self.hidden_dim
        for k in range(self.hidden_dim):
            backprop_sum = sum(self.W2[j][k] * delta2[j] for j in range(self.output_dim))
            delta1[k] = backprop_sum * self._sigmoid_derivative(a1[k])

        # 4. Update W2 and b2: dW2 = delta2 * a1^T
        for j in range(self.output_dim):
            for k in range(self.hidden_dim):
                self.W2[j][k] -= self.lr * delta2[j] * a1[k]
            self.b2[j] -= self.lr * delta2[j]

        # 5. Update W1 and b1: dW1 = delta1 * x^T
        for j in range(self.hidden_dim):
            for k in range(len(x)):
                self.W1[j][k] -= self.lr * delta1[j] * x[k]
            self.b1[j] -= self.lr * delta1[j]

    def predict(self, x: List[float]) -> List[float]:
        _, a2 = self.forward(x)
        return a2

if __name__ == "__main__":
    # Non-linear XOR problem
    X = [[0.0, 0.0], [0.0, 1.0], [1.0, 0.0], [1.0, 1.0]]
    Y = [[0.0], [1.0], [1.0], [0.0]]

    mlp = MultiLayerPerceptronFromScratch(input_dim=2, hidden_dim=4, output_dim=1, lr=0.8)
    for epoch in range(10000):
        for sample_x, sample_y in zip(X, Y):
            mlp.train_step(sample_x, sample_y)

    print("XOR Predictions after 10k epochs:")
    for sample_x in X:
        pred = mlp.predict(sample_x)[0]
        print(f"Input: {sample_x} -> Prediction: {pred:.4f} (Class {round(pred)})")

5. Active Knowledge Checks#

Concept check

In the matrix calculus derivation of backpropagation, why does the error vector delta^(l) transform by the transpose of the next layer's weight matrix (W^(l+1))^T?

Concept check

Why does the Adam optimizer require bias-correction terms (dividing m_t by 1 - beta_1^t and v_t by 1 - beta_2^t) during initial training steps?