Calculus and Gradients: How Networks Learn

Training a neural network means nudging millions or billions of numbers so a loss goes down. To know which way to nudge each number, you need its derivative: how much the loss changes when that number changes a little. Backpropagation, the algorithm behind every deep-learning framework, is nothing more than the chain rule from first-year calculus applied systematically to a graph of operations. This lesson builds that machinery from the ground up, works a full backward pass by hand, and shows how to check gradients the way professionals do when they write custom layers.

Intuition: sensitivity and blame

Imagine a factory line where each station transforms a part and passes it on, and a quality inspector at the end gives a score. If the score is bad, each station wants to know: "if I had changed my setting slightly, how much would the final score have moved?" The last station can answer directly. The station before it answers by multiplying its own local sensitivity by the sensitivity of everything downstream. That relay of blame, from the end of the line back to the start, is backpropagation. The derivative is the sensitivity; the chain rule is how sensitivities multiply along the line.

Derivatives

For a function $f: \mathbb{R} \to \mathbb{R}$, the derivative at $x$ is

It is the slope of the tangent line, and it gives the best linear approximation near $x$: $f(x + h) \approx f(x) + f'(x)\,h$ for small $h$. That approximation is the entire justification for gradient descent: if you step against the slope, the function decreases, at least for small steps.

Rules you will use constantly (with $c$ a constant):

ReLU is not differentiable at exactly $0$; frameworks simply pick $0$ there (a valid subgradient), and it causes no problems in practice. The sum rule $(f + g)' = f' + g'$ and product rule $(fg)' = f'g + fg'$ round out the toolkit.

Example: $f(x) = x^2$ at $x = 3$ has $f'(3) = 6$. Check with $h = 0.01$: $(3.01^2 - 9)/0.01 = 0.0601/0.01 = 6.01$, close to 6.

Partial derivatives and the gradient

Models have many parameters, so losses are functions $f: \mathbb{R}^n \to \mathbb{R}$. The partial derivative $\partial f / \partial x_i$ is the derivative with respect to $x_i$ while every other input is held fixed. Collecting all of them gives the gradient:

Worked example. $f(x, y) = x^2 y + 3y$.

• $\partial f/\partial x = 2xy$ (treat $y$ as a constant).
• $\partial f/\partial y = x^2 + 3$ (treat $x$ as a constant).
• At $(x, y) = (1, 2)$: $f = 2 + 6 = 8$ and $\nabla f = (4, 4)$.

Three facts make the gradient the central object of optimization:

1. Linear approximation: $f(\mathbf{x} + \boldsymbol{\delta}) \approx f(\mathbf{x}) + \nabla f(\mathbf{x})^\top \boldsymbol{\delta}$.
2. Steepest ascent: among all directions of unit length, $\nabla f$ points the way $f$ increases fastest, and $-\nabla f$ is the steepest descent. (By Cauchy–Schwarz, $\nabla f^\top \mathbf{u}$ over unit $\mathbf{u}$ is maximized when $\mathbf{u}$ is parallel to $\nabla f$.)
3. Stationary points: at a minimum, maximum or saddle point of a smooth function, $\nabla f = \mathbf{0}$.

The gradient has the same shape as the parameters: the gradient of a loss with respect to a $768 \times 3072$ weight matrix is itself a $768 \times 3072$ matrix. That is what param.grad holds in PyTorch.

The chain rule

If $y = f(u)$ and $u = g(x)$, then

With several intermediate variables, you sum over every path: if $L$ depends on $x$ through $u_1, \dots, u_k$, then

Worked example: one logistic neuron. $\hat{y} = \sigma(z)$ with $z = wx + b$, and $w = 2$, $x = 1$, $b = -1$.

1. Forward: $z = 2 \cdot 1 - 1 = 1$, $\hat{y} = \sigma(1) = 1/(1 + e^{-1}) \approx 0.7311$.
2. Local derivatives: $d\hat{y}/dz = \sigma(1)(1 - \sigma(1)) \approx 0.7311 \times 0.2689 \approx 0.1966$; $\partial z/\partial w = x = 1$; $\partial z / \partial b = 1$.
3. Chain: $\partial \hat{y}/\partial w = 0.1966 \times 1 = 0.1966$ and $\partial \hat y / \partial b = 0.1966$.

Notice that $\sigma'(z) \le 0.25$ everywhere. Chain twenty sigmoids and the product of local derivatives is at most $0.25^{20} \approx 9 \times 10^{-13}$: this is the vanishing gradient problem, and one reason ReLU (local derivative 1 when active) and residual connections replaced deep sigmoid stacks.

Computational graphs

Any computation can be drawn as a directed acyclic graph (DAG): nodes are intermediate values, edges show which values feed which operations. The forward pass evaluates nodes in order from inputs to output. The backward pass visits nodes in reverse order and, at each node $u$, computes the adjoint $\bar{u} = \partial L / \partial u$ by summing over the node's consumers $v$:

Each operation only needs to know its own local derivative; the graph takes care of the rest. That locality is why you can add a new layer type to a framework by writing just its forward function and its local backward rule.

Jacobians, briefly

For a vector-valued function $\mathbf{f}: \mathbb{R}^n \to \mathbb{R}^m$, the derivative is the Jacobian matrix $J \in \mathbb{R}^{m \times n}$ with $J_{ij} = \partial f_i / \partial x_j$. The chain rule becomes matrix multiplication: if $\mathbf{y} = \mathbf{f}(\mathbf{g}(\mathbf{x}))$, then $J_{\mathbf{y}, \mathbf{x}} = J_{\mathbf{f}} \, J_{\mathbf{g}}$.

Examples: for $\mathbf{y} = W\mathbf{x}$, the Jacobian is $W$ itself. For an elementwise activation, it is diagonal with entries $\phi'(x_i)$, which is why frameworks implement it as an elementwise multiply rather than building a matrix. For softmax, $\partial p_i / \partial z_j = p_i(\delta_{ij} - p_j)$ (with $\delta_{ij} = 1$ if $i = j$, else $0$).

Frameworks never build full Jacobians for large layers. Reverse mode computes vector–Jacobian products $\bar{\mathbf{y}}^\top J$ directly. For a batched linear layer $Y = XW$ with upstream gradient $G = \partial L / \partial Y$ (same shape as $Y$), the results are

Check the shapes: $X$ is $(B, d_{\text{in}})$, $G$ is $(B, d_{\text{out}})$, so $X^\top G$ is $(d_{\text{in}}, d_{\text{out}})$, matching $W$. When in doubt, derive matrix gradients by matching shapes, then confirm numerically.

Worked backpropagation through a 2-layer network

Network: input $\mathbf{x} \in \mathbb{R}^2$, hidden layer of 2 ReLU units, one linear output, squared-error loss.

Values: $\mathbf{x} = (1, 2)$, target $y = 1$,

Forward pass.

1. $z_1 = 0.1 \cdot 1 + 0.2 \cdot 2 = 0.5$; $z_2 = 0.3 \cdot 1 - 0.4 \cdot 2 = -0.5$.
2. $\mathbf{h} = (\max(0, 0.5), \max(0, -0.5)) = (0.5, 0)$.
3. $\hat{y} = 0.5 \cdot 0.5 + (-0.6) \cdot 0 + 0.1 = 0.35$.
4. $L = \tfrac12 (0.35 - 1)^2 = \tfrac12 (0.4225) = 0.21125$.

Backward pass (each line uses only the line above and a local derivative).

1. $\bar{\hat{y}} = \partial L / \partial \hat{y} = \hat{y} - y = -0.65$.
2. Output layer: $\partial L/\partial \mathbf{w}_2 = \bar{\hat y}\, \mathbf{h} = -0.65 \cdot (0.5, 0) = (-0.325, 0)$; $\partial L / \partial b_2 = -0.65$.
3. Into the hidden layer: $\bar{\mathbf{h}} = \bar{\hat y}\, \mathbf{w}_2 = -0.65 \cdot (0.5, -0.6) = (-0.325, 0.39)$.
4. Through ReLU (local derivative 1 where $z > 0$, 0 where $z < 0$): $\bar{\mathbf{z}} = (-0.325 \cdot 1, \; 0.39 \cdot 0) = (-0.325, 0)$.
5. First layer: $\partial L / \partial W_1 = \bar{\mathbf{z}}\, \mathbf{x}^\top$ (an outer product):

Interpretation: the second hidden unit is "dead" for this input (its pre-activation is negative), so neither its incoming weights nor $w_{2,2}$ receive any gradient. All gradients are negative where non-zero, so a gradient-descent step will increase those weights, raising $\hat{y}$ toward the target. With learning rate $\eta = 0.1$, the new $b_2$ would be $0.1 - 0.1 \cdot (-0.65) = 0.165$.

Automatic differentiation: forward vs reverse mode

There are three ways to get derivatives by computer:

• Symbolic differentiation manipulates formulas; it is exact but expressions can blow up in size.
• Numerical differentiation uses finite differences; it is easy but approximate and needs one or two function evaluations per input.
• Automatic differentiation (AD) applies the chain rule to the actual sequence of elementary operations your program executes, giving derivatives exact to floating-point precision.

AD has two modes. For $\mathbf{f}: \mathbb{R}^n \to \mathbb{R}^m$:

A training loss is a single scalar ($m = 1$) depending on billions of parameters ($n$ huge), so reverse mode computes the entire gradient for roughly the cost of a small constant number of forward passes (the usual rule of thumb is that a backward pass costs about twice the forward pass). Backpropagation is reverse-mode AD. The price is memory: activations from the forward pass must be kept for the backward pass, which is why activation memory dominates training and why activation checkpointing (recomputing some activations during the backward pass) trades compute for memory. Forward mode is used for Jacobian–vector products, for example in some second-order methods.

A reverse-mode autodiff engine in 30 lines

import math

class Value:
    def __init__(self, data, parents=(), local_grads=()):
        self.data, self.grad = data, 0.0
        self.parents, self.local_grads = parents, local_grads

    def __add__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        return Value(self.data + other.data, (self, other), (1.0, 1.0))

    def __mul__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        return Value(self.data * other.data, (self, other), (other.data, self.data))

    def relu(self):
        return Value(max(0.0, self.data), (self,), (1.0 if self.data > 0 else 0.0,))

    def sigmoid(self):
        s = 1 / (1 + math.exp(-self.data))
        return Value(s, (self,), (s * (1 - s),))

    def backward(self):
        order, seen = [], set()
        def visit(v):                      # topological sort of the graph
            if id(v) not in seen:
                seen.add(id(v))
                for p in v.parents: visit(p)
                order.append(v)
        visit(self)
        self.grad = 1.0
        for v in reversed(order):          # adjoints flow output -> inputs
            for p, g in zip(v.parents, v.local_grads):
                p.grad += v.grad * g

w, x, b = Value(2.0), Value(1.0), Value(-1.0)
y = (w * x + b).sigmoid()
y.backward()
print(round(y.data, 4), round(w.grad, 4), round(b.grad, 4))  # 0.7311 0.1966 0.1966

Each operation records its parents and its local derivatives during the forward pass; backward sorts the graph topologically and accumulates adjoints with += (the sum over paths in the chain rule). This reproduces the logistic-neuron numbers from earlier. PyTorch, JAX and TensorFlow do the same thing on tensors, with hand-written vector–Jacobian products for each operation. In PyTorch the equivalent is a tensor created with requires_grad=True, a forward computation, then loss.backward(), after which each leaf tensor's .grad holds its gradient.

Checking gradients numerically

When you implement a custom layer or a hand-derived gradient, verify it against the central difference

where $\mathbf{e}_i$ is the $i$-th unit vector. Its error shrinks like $O(h^2)$, versus $O(h)$ for the one-sided difference. Choose $h \approx 10^{-5}$ to $10^{-6}$ in float64: too large and the approximation error dominates; too small and floating-point cancellation does. Compare with the relative error

where $\epsilon$ is a tiny constant that avoids division by zero. In float64, values around $10^{-7}$ or below indicate a correct gradient; $10^{-2}$ or above means a bug.

import numpy as np

x, y = np.array([1., 2.]), 1.0
params = {"W1": np.array([[.1, .2], [.3, -.4]]), "b1": np.zeros(2),
          "w2": np.array([.5, -.6]), "b2": np.array(.1)}

def loss_and_grads(p):
    z = p["W1"] @ x + p["b1"]; h = np.maximum(z, 0)
    yhat = p["w2"] @ h + p["b2"]; L = 0.5 * (yhat - y) ** 2
    d_yhat = yhat - y
    d_z = d_yhat * p["w2"] * (z > 0)
    grads = {"W1": np.outer(d_z, x), "b1": d_z,
             "w2": d_yhat * h, "b2": np.array(d_yhat)}
    return L, grads

L, g = loss_and_grads(params)
print(L, g["W1"])                              # 0.21125 [[-0.325 -0.65] [0 0]]

h = 1e-6
for name, P in params.items():
    num = np.zeros_like(P)
    for i in np.ndindex(P.shape):
        old = P[i]
        P[i] = old + h; Lp, _ = loss_and_grads(params)
        P[i] = old - h; Lm, _ = loss_and_grads(params)
        P[i] = old
        num[i] = (Lp - Lm) / (2 * h)
    rel = np.abs(num - g[name]) / np.maximum(1e-12, np.abs(num) + np.abs(g[name]))
    print(name, rel.max())                      # all tiny (~1e-11)

The loop perturbs each parameter in place, recomputes the loss, and compares. Gradient checking is slow (two forward passes per parameter) so you run it on tiny inputs, once, as a unit test. Avoid checking exactly at ReLU kinks (pre-activations equal to zero), where the numerical derivative is meaningless; and disable dropout or other randomness while checking.

Common pitfalls

• Forgetting to zero gradients. Frameworks accumulate into .grad (the += above). In PyTorch, call optimizer.zero_grad() each step or gradients from previous batches leak in.
• Shape-mismatched matrix gradients. If your derived $\partial L/\partial W$ does not have $W$'s shape, a transpose is missing.
• Breaking the graph. Converting to NumPy, using .item(), or in-place edits on tensors needed for backward can silently stop gradients or raise errors.
• Vanishing and exploding gradients. Long chains of small (or large) local derivatives shrink (or blow up) the product; monitor per-layer gradient norms.
• Dead ReLUs. Units whose pre-activation is negative for all inputs get zero gradient forever, as unit 2 did in the worked example.
• Gradient checks in float32 or at kinks. Use float64 and avoid non-differentiable points, or you will chase false alarms.

In practice

• Let autodiff compute gradients; derive by hand only to understand, to write fused kernels, or to debug. Always back hand-derived gradients with a numerical check (PyTorch provides torch.autograd.gradcheck for this).
• Log the global gradient norm every step. Sudden spikes precede loss spikes; clipping the global norm (often to 1.0 in LLM training recipes) is a common safeguard.
• Memory, not compute, is often the limit in training: activations stored for reverse mode scale with batch size, sequence length and depth. Activation checkpointing, smaller micro-batches and mixed precision are the standard levers.
• Use torch.no_grad() (or the equivalent) during evaluation and inference so no graph is built.