Derivatives, Gradients and the Chain Rule

What this module teaches:

A derivative answers one question: if I nudge this input slightly, how much does the output change?

import numpy as np

# Loss function: L = w^2
# The derivative of w^2 with respect to w is 2w (calculus rule)
# At w=3: derivative = 2*3 = 6

def loss(w):
    return w ** 2

def derivative_of_loss(w):
    return 2 * w   # exact derivative from calculus

# Numerical approximation — the definition of a derivative
# derivative = (f(w + tiny) - f(w)) / tiny   as tiny → 0
def numerical_derivative(f, w, h=1e-5):
    return (f(w + h) - f(w)) / h

w = 3.0
print(f"Loss at w=3:          {loss(w):.4f}")          # 9.0
print(f"Exact derivative:     {derivative_of_loss(w):.4f}")   # 6.0
print(f"Numerical derivative: {numerical_derivative(loss, w):.4f}")  # ~6.0

# What the derivative tells us:
deriv = derivative_of_loss(w)
if deriv > 0:
    print(f"Derivative is positive ({deriv:.1f}) → increase w increases loss → DECREASE w")
elif deriv < 0:
    print(f"Derivative is negative ({deriv:.1f}) → increase w decreases loss → INCREASE w")
else:
    print("Derivative is 0 → we are at a minimum!")

# One gradient descent step
learning_rate = 0.1
w_new = w - learning_rate * deriv   # move opposite to derivative
print(f"w before: {w:.2f}, w after one step: {w_new:.2f}")  # 3.0 → 2.4
print(f"Loss before: {loss(w):.4f}, loss after: {loss(w_new):.4f}")   # 9.0 → 5.76

import numpy as np

# Simple linear regression: predict delivery time from distance
# Loss = Mean Squared Error = (1/n) * sum((y_pred - y_true)^2)
# y_pred = w * x + b

np.random.seed(42)
n = 100
x = np.random.uniform(0.5, 8.0, n)           # distance in km
y_true = 8.6 + 7.3 * x + np.random.randn(n) * 3  # actual delivery time

# Initialise weights randomly
w = np.random.randn()   # slope
b = np.random.randn()   # intercept

def predict(w, b, x):
    return w * x + b

def mse_loss(y_pred, y_true):
    return np.mean((y_pred - y_true) ** 2)

# ── Compute gradients manually ────────────────────────────────────────
# Loss L = (1/n) * sum((w*x + b - y)^2)
# ∂L/∂w = (2/n) * sum((w*x + b - y) * x)   ← partial derivative w.r.t. w
# ∂L/∂b = (2/n) * sum((w*x + b - y))        ← partial derivative w.r.t. b

def compute_gradients(w, b, x, y_true):
    n = len(x)
    y_pred = predict(w, b, x)
    error  = y_pred - y_true           # residuals: shape (n,)

    dL_dw = (2/n) * np.sum(error * x)  # gradient w.r.t. weight
    dL_db = (2/n) * np.sum(error)       # gradient w.r.t. bias

    return dL_dw, dL_db

# One gradient step
dL_dw, dL_db = compute_gradients(w, b, x, y_true)
print(f"Initial w={w:.3f}, b={b:.3f}")
print(f"Gradient dL/dw={dL_dw:.3f}, dL/db={dL_db:.3f}")

learning_rate = 0.01
w = w - learning_rate * dL_dw   # step opposite to gradient
b = b - learning_rate * dL_db

y_pred = predict(w, b, x)
print(f"Loss after one step: {mse_loss(y_pred, y_true):.4f}")

# The gradient of w tells us: for every 1-unit increase in w,
# loss increases by dL_dw units.
# Since we want to DECREASE loss, we move in the OPPOSITE direction.

# Complete gradient descent training loop for Swiggy delivery prediction

np.random.seed(42)
n = 1000
x = np.random.uniform(0.5, 8.0, n)
y_true = 8.6 + 7.3 * x + np.random.randn(n) * 3

w = 0.0   # start from scratch
b = 0.0

learning_rate = 0.01
batch_size    = 32
epochs        = 100
losses        = []

n_batches = n // batch_size

for epoch in range(epochs):
    # Shuffle data at the start of each epoch
    # Shuffling prevents the model from memorising the order
    indices = np.random.permutation(n)
    x_shuffled = x[indices]
    y_shuffled = y_true[indices]

    epoch_loss = 0.0

    for i in range(n_batches):
        # Extract one mini-batch
        start = i * batch_size
        end   = start + batch_size
        x_batch = x_shuffled[start:end]   # shape (32,)
        y_batch = y_shuffled[start:end]   # shape (32,)

        # Forward pass: compute predictions
        y_pred = w * x_batch + b          # shape (32,)

        # Compute loss on this batch
        batch_loss = np.mean((y_pred - y_batch) ** 2)
        epoch_loss += batch_loss

        # Backward pass: compute gradients
        error = y_pred - y_batch
        dL_dw = (2 / batch_size) * np.sum(error * x_batch)
        dL_db = (2 / batch_size) * np.sum(error)

        # Update weights — move opposite to gradient
        w = w - learning_rate * dL_dw
        b = b - learning_rate * dL_db

    losses.append(epoch_loss / n_batches)

    if epoch % 10 == 0:
        rmse = np.sqrt(losses[-1])
        print(f"Epoch {epoch:3d} | Loss: {losses[-1]:.4f} | RMSE: {rmse:.2f} min")

print(f"
Learned: w = {w:.4f} (true: 7.3), b = {b:.4f} (true: 8.6)")
# Epoch   0 | Loss: 312.4  | RMSE: 17.67 min   ← terrible, just started
# Epoch  10 | Loss: 42.8   | RMSE: 6.54 min
# Epoch  50 | Loss: 10.2   | RMSE: 3.19 min
# Epoch  90 | Loss: 9.1    | RMSE: 3.02 min    ← close to true noise level
# Learned: w = 7.2891 (true: 7.3), b = 8.6432 (true: 8.6)  ← recovered!

The chain rule in plain English: to find how much the final output changes when an early input changes, multiply the sensitivities at each step along the chain.

# Chain rule implemented manually for a tiny 2-layer network

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def sigmoid_derivative(x):
    s = sigmoid(x)
    return s * (1 - s)   # derivative of sigmoid

def relu(x):
    return np.maximum(0, x)

def relu_derivative(x):
    return (x > 0).astype(float)   # 1 where x > 0, 0 elsewhere

np.random.seed(42)
n = 32

# Network: input(4) → hidden(8) → output(1)
x = np.random.randn(n, 4)
y = np.random.randn(n, 1)

W1 = np.random.randn(4, 8) * 0.1
b1 = np.zeros((1, 8))
W2 = np.random.randn(8, 1) * 0.1
b2 = np.zeros((1, 1))

lr = 0.01

for step in range(5):
    # ── Forward pass ──────────────────────────────────────────────────
    Z1 = x @ W1 + b1          # (32,4) @ (4,8) → (32,8)
    A1 = relu(Z1)              # (32,8) — apply ReLU

    Z2 = A1 @ W2 + b2          # (32,8) @ (8,1) → (32,1)
    y_pred = Z2                # linear output (regression)

    loss = np.mean((y_pred - y) ** 2)

    # ── Backward pass — chain rule applied backwards ───────────────────
    # Start at loss, work backwards to W1

    # Gradient of loss w.r.t. output (∂L/∂y_pred)
    dL_dZ2 = 2 * (y_pred - y) / n         # (32,1)

    # Gradient w.r.t. W2 (∂L/∂W2) = A1ᵀ @ dL_dZ2
    dL_dW2 = A1.T @ dL_dZ2               # (8,32) @ (32,1) → (8,1)
    dL_db2 = dL_dZ2.sum(axis=0)

    # Gradient w.r.t. A1 — chain through W2
    dL_dA1 = dL_dZ2 @ W2.T               # (32,1) @ (1,8) → (32,8)

    # Gradient through ReLU (∂A1/∂Z1 = relu_derivative(Z1))
    dL_dZ1 = dL_dA1 * relu_derivative(Z1)  # element-wise: chain rule

    # Gradient w.r.t. W1 (∂L/∂W1) = xᵀ @ dL_dZ1
    dL_dW1 = x.T @ dL_dZ1               # (4,32) @ (32,8) → (4,8)
    dL_db1 = dL_dZ1.sum(axis=0)

    # ── Weight update ──────────────────────────────────────────────────
    W2 -= lr * dL_dW2
    b2 -= lr * dL_db2
    W1 -= lr * dL_dW1
    b1 -= lr * dL_db1

    print(f"Step {step+1}: loss = {loss:.6f}")

# This is EXACTLY what PyTorch's loss.backward() does automatically.
# The only difference: PyTorch builds a computation graph and
# differentiates it symbolically. Same chain rule. No magic.

# Partial derivatives — isolating each weight's effect

# Loss function with two weights: L = (w1*x1 + w2*x2 - y)^2
# ∂L/∂w1 = 2*(w1*x1 + w2*x2 - y) * x1   (treat w2 as constant)
# ∂L/∂w2 = 2*(w1*x1 + w2*x2 - y) * x2   (treat w1 as constant)

x1, x2 = 3.0, 1.5   # two features: distance, traffic score
y = 35.0             # actual delivery time
w1, w2 = 6.0, 2.0   # current weights

def loss_2d(w1, w2):
    y_pred = w1 * x1 + w2 * x2
    return (y_pred - y) ** 2

def partial_w1(w1, w2):
    error = (w1 * x1 + w2 * x2 - y)
    return 2 * error * x1

def partial_w2(w1, w2):
    error = (w1 * x1 + w2 * x2 - y)
    return 2 * error * x2

dL_dw1 = partial_w1(w1, w2)
dL_dw2 = partial_w2(w1, w2)

print(f"Prediction: {w1*x1 + w2*x2:.1f} min (actual: {y} min)")
print(f"∂L/∂w1 = {dL_dw1:.3f}  → w1 change contributes {abs(dL_dw1):.1f}× to loss change")
print(f"∂L/∂w2 = {dL_dw2:.3f}  → w2 change contributes {abs(dL_dw2):.1f}× to loss change")
print(f"
Gradient vector: [{dL_dw1:.3f}, {dL_dw2:.3f}]")

# Verify with numerical approximation
h = 1e-5
numerical_w1 = (loss_2d(w1 + h, w2) - loss_2d(w1, w2)) / h
numerical_w2 = (loss_2d(w1, w2 + h) - loss_2d(w1, w2)) / h
print(f"
Numerical check: [{numerical_w1:.3f}, {numerical_w2:.3f}]")
# Should match partial derivatives above ✓

# Gradient step
lr = 0.01
w1_new = w1 - lr * dL_dw1
w2_new = w2 - lr * dL_dw2
print(f"
Before: w1={w1:.2f}, w2={w2:.2f}")
print(f"After:  w1={w1_new:.2f}, w2={w2_new:.2f}")
print(f"Loss before: {loss_2d(w1, w2):.2f}, after: {loss_2d(w1_new, w2_new):.2f}")

# How autograd works conceptually — using NumPy to simulate it

# In PyTorch you would do:
# import torch
# w = torch.tensor(3.0, requires_grad=True)
# loss = w ** 2
# loss.backward()
# print(w.grad)   # tensor(6.0)  ← ∂L/∂w = 2w = 6

# We'll simulate this with a tiny autograd engine in NumPy

class Value:
    """A scalar value that remembers how it was computed."""
    def __init__(self, data, _children=(), _op='', label=''):
        self.data = data
        self.grad = 0.0            # gradient starts at 0
        self._backward = lambda: None
        self._prev = set(_children)
        self.label = label

    def __mul__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        out = Value(self.data * other.data, (self, other), '*')
        def _backward():
            # Chain rule: d(self*other)/d(self) = other
            self.grad  += other.data * out.grad
            other.grad += self.data  * out.grad
        out._backward = _backward
        return out

    def __pow__(self, n):
        out = Value(self.data ** n, (self,), f'**{n}')
        def _backward():
            self.grad += n * (self.data ** (n-1)) * out.grad
        out._backward = _backward
        return out

    def backward(self):
        # Topological sort then call _backward in reverse
        topo, visited = [], set()
        def build(v):
            if v not in visited:
                visited.add(v)
                for child in v._prev: build(child)
                topo.append(v)
        build(self)
        self.grad = 1.0   # dL/dL = 1.0 (start of chain)
        for node in reversed(topo):
            node._backward()

# Compute loss = w^2 and get gradient automatically
w    = Value(3.0, label='w')
loss = w ** 2

print(f"Loss = {loss.data}")   # 9.0
loss.backward()
print(f"∂L/∂w = {w.grad}")    # 6.0  ← correct! (2w = 2×3 = 6)

# More complex: L = (w1*x + b - y)^2
w1   = Value(6.0,  label='w1')
x    = Value(3.0,  label='x')
b    = Value(1.0,  label='b')
y    = Value(35.0, label='y')

y_pred = w1 * x * Value(1.0) * Value(1.0)  # simplified: w1*x
error  = y_pred * Value(1.0)               # just y_pred for demo
loss2  = error ** 2

loss2.backward()
print(f"w1.grad = {w1.grad:.3f}")   # ∂L/∂w1
# This is exactly what PyTorch does with arbitrary networks

# PyTorch autograd — the real thing (shown for reference)
# Run this when you install PyTorch in the Deep Learning section

# import torch
#
# x  = torch.tensor([3.2, 2.0, 15.0, 7.0])       # one Swiggy order
# W  = torch.randn(4, 1, requires_grad=True)       # weights to learn
# b  = torch.zeros(1,    requires_grad=True)       # bias to learn
# y  = torch.tensor([37.0])                        # actual delivery time
#
# # Forward pass
# y_pred = x @ W + b
# loss   = ((y_pred - y) ** 2).mean()
#
# # Backward pass — computes ALL gradients automatically
# loss.backward()
#
# print(W.grad)   # ∂L/∂W — shape (4,1)
# print(b.grad)   # ∂L/∂b — shape (1,)
#
# # Gradient descent step
# with torch.no_grad():   # don't track the update itself
#     W -= 0.01 * W.grad
#     b -= 0.01 * b.grad
#     W.grad.zero_()      # clear gradients — IMPORTANT: they accumulate by default
#     b.grad.zero_()

# The three lines every PyTorch training loop always has:
# 1. optimizer.zero_grad()   — clear gradients from last step
# 2. loss.backward()         — compute new gradients
# 3. optimizer.step()        — update weights using gradients

# Demonstrating vanishing gradients with sigmoid vs ReLU

def sigmoid(x): return 1 / (1 + np.exp(-x))
def sigmoid_grad(x):
    s = sigmoid(x)
    return s * (1 - s)   # max value is 0.25 at x=0

def relu(x): return np.maximum(0, x)
def relu_grad(x): return (x > 0).astype(float)   # 1 or 0

# Simulate gradient flowing backward through 10 layers
x = np.array([0.5])   # some activation value

sigmoid_gradient = 1.0
relu_gradient    = 1.0

print("Gradient magnitude after each layer:")
print(f"{'Layer':<8} {'Sigmoid':<20} {'ReLU'}")
print("─" * 45)

for layer in range(1, 11):
    sigmoid_gradient *= sigmoid_grad(x)[0]
    relu_gradient    *= relu_grad(x)[0]
    print(f"{layer:<8} {sigmoid_gradient:<20.8f} {relu_gradient:.8f}")

# Layer 1  : Sigmoid = 0.23500 | ReLU = 1.0
# Layer 5  : Sigmoid = 0.00072 | ReLU = 1.0
# Layer 10 : Sigmoid = 0.000000 | ReLU = 1.0
# ← Sigmoid gradient vanishes. ReLU gradient stays at 1.0.

# Solutions to vanishing gradients:
# 1. Use ReLU or GELU activations instead of sigmoid/tanh
# 2. Use residual/skip connections (ResNet) — gradient highway
# 3. Use batch normalisation — normalises activations between layers
# 4. Careful weight initialisation (Xavier/He) — keeps activations in range

# Solutions to exploding gradients:
# 1. Gradient clipping — cap gradient magnitude above a threshold
# 2. Smaller learning rate
# 3. Batch normalisation

# Gradient clipping example:
large_gradient = np.array([150.0, -89.3, 203.1, -45.6])
max_norm = 1.0

gradient_norm = np.linalg.norm(large_gradient)
if gradient_norm > max_norm:
    clipped = large_gradient * (max_norm / gradient_norm)
    print(f"
Before clipping: norm = {gradient_norm:.1f}")
    print(f"After clipping:  norm = {np.linalg.norm(clipped):.4f}")
    print(f"Direction preserved: {np.allclose(clipped / np.linalg.norm(clipped),
                                               large_gradient / gradient_norm)}")