Easy 5 Steps Derivation of Neural Network (Math to Numpy code)

🎯 Bridging the gap between multivariable calculus and clean, vectorised implementation from scratch.

Background

Neural networks can feel daunting when presented as a massive wall of summations and matrices. In this guide, we break down the exact mathematical derivation of an Artificial Neural Network (ANN) into 5 simple steps, connecting every single equation directly to raw NumPy code.

Core Idea: A neural network extends logistic regression. In logistic regression, we compute a single weighted sum followed by an activation function (like sigmoid). A neural network stacks multiple layers of these computations, passing activation outputs sequentially from layer to layer.

Single Layer (Logistic Regression) $$a = \sigma(W x + b)$$
→
Multi-Layer (Neural Network) $$a^{[l]} = g^{[l]}\left(W^{[l]} a^{[l-1]} + b^{[l]}\right)$$

Derivation

Suppose we have $m$ number of houses and $n$ features describing each house represented by the data as $X_1, X_2, \dots, X_n$. We need to predict the $y$ price of each house.

To do that we want to apply a neural network having $L$ layers, where the first is input $X$ and the last is output layer $L$. $Z$ is the weighted sum, and $g$ denotes the activation function such as sigmoid, tanh, or ReLU. $A$ is the result of the activation function.

Step 1: The Forward Pass

For each layer $l = 1, 2, \dots, L$, we compute the linear combination followed by the non-linear activation:

Forward Pass (X ⟶ ŷ) Input (X) Hidden 1 Hidden l Hidden L-1 Output (ŷ) X₁ X₂ X₃ ⋮ Xₙ ⋮ ⋮ ⋮ ŷ
Figure 1: Architecture of a Deep Neural Network showing forward information flow from input features $X_1, X_2, \dots, X_n$ through hidden layer $l$ to output prediction ŷ.

$$\text{Linear sum: } Z^{[l]} = W^{[l]} A^{[l-1]} + b^{[l]}$$

$$\text{Activation: } A^{[l]} = g^{[l]}(Z^{[l]})$$

$$\text{where, } A^{[0]} = X \quad \text{and} \quad A^{[L]} = ŷ$$

From Scratch with NumPy:

def forward_pass(self, X):
    cache = {"A0": X}
    A_prev = X

    for l in range(1, self.L + 1):
        Z = A_prev @ self.weights[l - 1] + self.biases[l - 1]

        # use sigmoid activation for all layers
        A = self.sigmoid(Z)

        cache[f"Z{l}"] = Z
        cache[f"A{l}"] = A
        A_prev = A

    return cache

Step 2: The Loss Function

After computing the predictions, we obviously get some deviation from the true values $y$. So, we need to calculate error using the cost function $J$ (average loss over $m$ houses):

$$J = \frac{1}{m} \sum_{i=1}^{m} L(A^{[L]}, y)$$

$$\text{where Mean Squared Error Loss: } L = \frac{1}{2} (A^{[L]} - y)^2$$

$$\text{Overall Cost: } J = \frac{1}{2m} \sum_{i=1}^{m} (A^{[L]} - y)^2$$

From Scratch with NumPy:

def compute_loss(self, A_L, y):
    return np.mean((A_L - y) ** 2)

Step 3: The Goal of Backpropagation

Based on cost function $J$, we need to update our weights and biases, for that we need $\frac{\partial J}{\partial W^{[l]}}$ and $\frac{\partial J}{\partial b^{[l]}}$. To compute this efficiently across deep layers, we define the error term $\delta^{[l]}$ for layer $l$ as the derivative of the cost $J$ with respect to the linear pre-activation $Z^{[l]}$:

$$\delta^{[l]} = \frac{\partial J}{\partial Z^{[l]}}$$

During backpropagation, error flows in reverse. Each neuron in layer $l$ receives a weighted error contribution from every connected neuron in layer $l+1$. These individual contributions are summed up ($\Sigma$) and multiplied by the local activation derivative $g'(Z^{[l]})$ to form $\delta^{[l]}$:

Input (X) Layer l Layer l+1 Output (ŷ) W₁ⱼ · δ₁⁽ˡ⁺¹⁾ W₂ⱼ · δ₂⁽ˡ⁺¹⁾ Wₖⱼ · δₖ⁽ˡ⁺¹⁾ X₁ X₂ ⋮ Xₙ Σ ⋮ δ₁ δ₂ ⋮ δₖ δ⁽ᴸ⁾ Neuron j: δⱼ⁽ˡ⁾
Figure 2: Backpropagation in a full network. The error $\delta^{(L)}$ originates at the output layer and flows backward. For a neuron $j$ in Layer $l$, error contributions $(W_{kj} \cdot \delta_k^{(l+1)})$ from all connected neurons in Layer $l+1$ converge and are summed ($\Sigma$) to form $\delta_j^{(l)}$.

Step 4: Step-by-Step Gradient Calculation

Step 4A: Output Layer Error ($\delta^{[L]}$)

Using the chain rule at the final output layer $L$:

$$\delta^{[L]} = \frac{\partial J}{\partial A^{[L]}} \cdot \frac{\partial A^{[L]}}{\partial Z^{[L]}}$$

$$\delta^{[L]} = (A^{[L]} - y) \odot {g^{[L]}}'(Z^{[L]})$$

Step 4B: Hidden Layer Error ($\delta^{[l]}$)

For layers $l = L-1, L-2, \dots, 1$, a change in $Z^{[l]}$ affects the cost through all neurons in the next layer ($l+1$). Using the chain rule:

$$\delta^{[l]} = \frac{\partial J}{\partial Z^{[l]}} = \left[ \sum_k \frac{\partial J}{\partial Z_k^{[l+1]}} \cdot \frac{\partial Z_k^{[l+1]}}{\partial A^{[l]}} \right] \cdot \frac{\partial A^{[l]}}{\partial Z^{[l]}}$$

$$\text{In matrix form: } \delta^{[l]} = \left( (W^{[l+1]})^T \delta^{[l+1]} \right) \odot {g^{[l]}}'(Z^{[l]})$$

Step 4C: Gradients for Weights and Biases

(i) Weight Gradient ($dW^{[l]}$):

$$\frac{\partial J}{\partial W^{[l]}} = \frac{\partial J}{\partial Z^{[l]}} \cdot \frac{\partial Z^{[l]}}{\partial W^{[l]}} = \frac{1}{m} \delta^{[l]} (A^{[l-1]})^T$$

(ii) Bias Gradient ($db^{[l]}$):

$$\frac{\partial J}{\partial b^{[l]}} = \frac{\partial J}{\partial Z^{[l]}} \cdot \frac{\partial Z^{[l]}}{\partial b^{[l]}} = \frac{1}{m} \sum \delta^{[l]}$$

From Scratch with NumPy:

m = y.shape[0]
A_L = cache[f"A{self.L}"]

deltas = {}

# use derivative of sigmoid
deltas[f"delta{self.L}"] = (A_L - y) * A_L * (1 - A_L)

for l in reversed(range(1, self.L)):
    delta_next = deltas[f"delta{l + 1}"]
    deltas[f"delta{l}"] = (delta_next @ self.weights[l].T) * cache[f"A{l}"] * (1 - cache[f"A{l}"])

Step 5: Weight Update Rules

We update the weights using Gradient Descent (learning rate = $\alpha$):

$$W^{[l]}_{\text{new}} = W^{[l]}_{\text{old}} - \alpha \cdot \frac{\partial J}{\partial W^{[l]}}$$

$$\text{and}$$

$$b^{[l]}_{\text{new}} = b^{[l]}_{\text{old}} - \alpha \cdot \frac{\partial J}{\partial b^{[l]}}$$

From Scratch with NumPy:

for l in range(1, self.L + 1):
    A_prev = cache[f"A{l - 1}"]
    delta = deltas[f"delta{l}"]

    dW = A_prev.T @ delta / m
    db = np.sum(delta, axis=0, keepdims=True) / m

    self.weights[l - 1] -= learning_rate * dW
    self.biases[l - 1] -= learning_rate * db

Complete Scratch Implementation

From Scratch with NumPy:

import numpy as np

class SimpleNeuralNetwork:
    def __init__(self, layer_sizes):
        self.layer_sizes = layer_sizes
        self.L = len(layer_sizes) - 1
        self.weights = []
        self.biases = []

        for l in range(1, self.L + 1):
            self.weights.append(np.random.randn(layer_sizes[l - 1], layer_sizes[l]) * 0.01)
            self.biases.append(np.zeros((1, layer_sizes[l])))

    def sigmoid(self, z):
        return 1 / (1 + np.exp(-z))

    def forward_pass(self, X):
        cache = {"A0": X}
        A_prev = X

        for l in range(1, self.L + 1):
            Z = A_prev @ self.weights[l - 1] + self.biases[l - 1]
            A = self.sigmoid(Z)

            cache[f"Z{l}"] = Z
            cache[f"A{l}"] = A
            A_prev = A

        return cache

    def compute_loss(self, A_L, y):
        return np.mean((A_L - y) ** 2)

    def backward_pass(self, cache, y):
        m = y.shape[0]
        grads = {}
        A_L = cache[f"A{self.L}"]

        deltas = {}
        deltas[f"delta{self.L}"] = (A_L - y) * A_L * (1 - A_L)

        for l in reversed(range(1, self.L)):
            delta_next = deltas[f"delta{l + 1}"]
            deltas[f"delta{l}"] = (delta_next @ self.weights[l].T) * cache[f"A{l}"] * (1 - cache[f"A{l}"])

        for l in range(1, self.L + 1):
            A_prev = cache[f"A{l - 1}"]
            delta = deltas[f"delta{l}"]

            grads[f"dW{l}"] = (A_prev.T @ delta) / m
            grads[f"db{l}"] = np.sum(delta, axis=0, keepdims=True) / m

        return grads

    def update_parameters(self, grads, learning_rate):
        for l in range(1, self.L + 1):
            self.weights[l - 1] -= learning_rate * grads[f"dW{l}"]
            self.biases[l - 1] -= learning_rate * grads[f"db{l}"]

    def fit(self, X, y, epochs=1000, learning_rate=0.1):
        for epoch in range(epochs):
            cache = self.forward_pass(X)
            A_L = cache[f"A{self.L}"]
            loss = self.compute_loss(A_L, y)
            grads = self.backward_pass(cache, y)
            self.update_parameters(grads, learning_rate)

            if epoch % 100 == 0:
                print(f"Epoch {epoch}: loss = {loss:.4f}")

Ready to build more algorithms from scratch? Check out Articles for more math-to-code tutorials!