Derivatives and Computation Graphs for Neural Network Learning

A computation graph turns an expression into a sequence of small operations. Here we differentiate one scalar output, \(J\), with respect to its inputs. The forward pass evaluates the operations in dependency order. The backward pass works toward the inputs, combining local derivatives to calculate how each value affects \(J\).

When I first learned backpropagation, the final derivative formulas felt disconnected from the computation that produced them. Drawing the function as a graph made the process more concrete: each node performs a simple local operation, and the chain rule connects those local derivatives.

The derivative as a local slope

A derivative describes how an output responds to a small change in an input. The derivative of \(f\) at a point \(a\) answers one question: if I nudge the input by a tiny amount, how much does the output move, and in what proportion? Formally \(f'(a)=\lim_{h\to 0}\frac{f(a+h)-f(a)}{h}\), When the derivative exists, the local approximation is \(f(a+h)\approx f(a)+f\prime(a)h\). A derivative of \(3\) predicts an output increase of about \(0.003\) for an input increase of \(0.001\), provided that step is small enough for the local approximation to be accurate.

For \(f(a)=4a\), the slope is 4 everywhere. Constant local derivatives occur in useful network operations, including addition and scaling. Other local derivatives depend on the values produced by the forward pass.

Local derivatives can depend on the input

For \(f(a)=a^2\), the derivative is \(2a\): at \(a=2\) a small change is amplified by about 4, and at \(a=5\) by about 10. For \(f(a)=\log a\), defined here for \(a>0\), the derivative is \(1/a\), so its magnitude grows as \(a\) approaches zero from above.

The gradient of the full loss generally depends on the current parameters, so it is evaluated again after an update. Its formula can be derived once; its value changes. This also happens with squared error on a linear model. Sigmoid provides another example: its local derivative \(\sigma(z)(1-\sigma(z))\) reaches 0.25 at \(z=0\) and approaches zero as \(|z|\) grows. Saturation can therefore reduce the gradient passed through that operation; the full parameter gradient also depends on the rest of the graph.

import numpy as np

def numeric_slope(f, a, h=1e-6):
    return (f(a + h) - f(a - h)) / (2 * h)

print(round(numeric_slope(lambda a: 4 * a, 2.0), 8))
# 4.0
print(round(numeric_slope(lambda a: a ** 2, 2.0), 8))
# 4.0
print(round(numeric_slope(lambda a: a ** 2, 5.0), 8))
# 10.0
print(round(numeric_slope(np.log, 2.0), 8))
# 0.5

The examples use \(d(a^2)/da=2a\), \(d(\log a)/da=1/a\) for \(a>0\), and the local rules for addition, multiplication, and scaling. The exercises also use sigmoid and ReLU; their needed derivative rules are given there.

What is a computation graph?

A computation graph is a directed graph that represents dependencies between values. Input nodes hold variables or parameters, operation nodes apply functions, and edges show how results flow into later operations.

• Forward pass: calculate intermediate values and the final output.
• Backward pass: calculate derivatives from the output back toward earlier values.
• Stored values: reuse useful forward-pass results when evaluating local derivatives.
• Gradient accumulation: add contributions when one variable affects the output through multiple paths.

This backward calculation is reverse-mode automatic differentiation. Backpropagation is its application to neural-network computations.

Breaking a function into local operations

Consider the following function:

\[J=3(a+bc)\]

Instead of treating it as one expression, introduce two intermediate variables:

\[u=bc\]

\[v=a+u\]

\[J=3v\]

The dependency structure is:

\[(b,c)\rightarrow u,\qquad(a,u)\rightarrow v,\qquad v\rightarrow J\]

b ──┐
    × ── u ──┐
c ──┘        + ── v ── ×3 ── J
         a ──┘

The variable \(a\) enters at the addition node, while \(b\) and \(c\) first meet at the multiplication node. This detail determines the paths followed during the backward pass.

Forward pass

Let:

\[a=5,\qquad b=3,\qquad c=2\]

Following the graph from left to right gives:

\[u=bc=3\times2=6\]

\[v=a+u=5+6=11\]

\[J=3v=3\times11=33\]

The forward pass produces both the final output \(J=33\) and the intermediate values \(u=6\) and \(v=11\). In this particular graph, the backward pass needs the input values \(b\) and \(c\); the local derivatives of addition and scaling are constants, so \(u\) and \(v\) are not needed to evaluate them. Other operations do need intermediate results, such as the sigmoid output used in its derivative. A differentiation system saves the values required by its backward rules, or recomputes them when trading extra computation for less stored data.

Local derivatives and upstream gradients

Each operation has a local derivative that describes its immediate relationship to its inputs.

Here, a partial derivative such as \(\partial u/\partial b\) varies \(b\) while holding the other input \(c\) fixed. An upstream gradient is the derivative arriving from operations closer to the final output. A node multiplies that upstream gradient by its local derivative. In code, a variable such as du conventionally means \(\partial J/\partial u\), not merely the local derivative of the multiplication node.

The chain rule

In this graph, every effect of \(a\) on \(J\) passes through \(v\), so the chain rule gives:

\[\frac{\partial J}{\partial a}=\frac{\partial J}{\partial v}\frac{\partial v}{\partial a}\]

The local effects are multiplied along a path. This is the central mathematical rule used by backpropagation.

Backward pass step by step

1. Seed the output gradient

Reverse-mode differentiation begins at the output with:

\[\frac{\partial J}{\partial J}=1\]

This value is the seed gradient. It expresses the fact that a small change in \(J\) changes \(J\) by the same amount.

2. Backpropagate through the scaling node

Because \(J=3v\):

\[d_v=\frac{\partial J}{\partial v}=\frac{\partial J}{\partial J}\frac{\partial J}{\partial v}=1\times3=3\]

3. Backpropagate through the addition node

An addition node has local derivative one with respect to each input:

\[d_a=d_v\frac{\partial v}{\partial a}=3\times1=3\]

\[d_u=d_v\frac{\partial v}{\partial u}=3\times1=3\]

4. Backpropagate through the multiplication node

For \(u=bc\), each input receives the upstream gradient multiplied by the other input:

\[d_b=d_u\frac{\partial u}{\partial b}=d_uc=3\times2=6\]

\[d_c=d_u\frac{\partial u}{\partial c}=d_ub=3\times3=9\]

The complete result is:

Gradient accumulation across multiple paths

Multiplication along a path is only half of the rule. If a variable affects the output through multiple paths, the gradient contributions from those paths must be added.

Consider:

\[J=x^2+x\]

The variable \(x\) reaches \(J\) through both the square term and the direct addition term. Their contributions are:

\[\frac{dJ}{dx}=2x+1\]

At \(x=3\), the two paths contribute 6 and 1, so the accumulated gradient is 7. Neural networks reuse parameters and activations in many operations, which makes gradient accumulation essential.

In code, start a gradient accumulator at zero and add each contribution. Writing the square as \(u=x\cdot x\) exposes two input edges from the same variable. Each edge contributes \(x\), and the direct addition contributes 1:

x = 3.0
u = x * x
J = u + x

dJ = 1.0
du = dJ
dx = 0.0
dx += dJ
dx += du * x
dx += du * x
print("J:", J)
print("dx:", dx)
# J: 12.0
# dx: 7.0

The repeated du * x lines represent the two input positions of the multiplication node. Replacing += with assignment would discard earlier contributions. In a larger graph, process a node only after its downstream uses have contributed to its gradient. Reversing a valid forward dependency order gives such an order; there is no need to enumerate every complete path separately.

Numerical gradient checking

An analytical gradient can be compared with a finite-difference approximation. For a sufficiently smooth function, central differences have a smaller truncation error order than first-order one-sided differences. Floating-point rounding and the chosen step size also affect accuracy.

\[\frac{\partial J}{\partial x}\approx\frac{J(x+\epsilon)-J(x-\epsilon)}{2\epsilon}\]

The approximation is not the gradient used for ordinary training. Coordinate-by-coordinate central differences require two forward evaluations per input or parameter. This becomes expensive for a large parameter vector; a small example makes the comparison affordable.

def objective(a, b, c):
    return 3 * (a + b * c)


def central_difference(variable, a, b, c, epsilon=1e-6):
    values_plus = {"a": a, "b": b, "c": c}
    values_minus = values_plus.copy()
    values_plus[variable] += epsilon
    values_minus[variable] -= epsilon

    j_plus = objective(**values_plus)
    j_minus = objective(**values_minus)
    return (j_plus - j_minus) / (2 * epsilon)


a, b, c = 5.0, 3.0, 2.0

numerical = {
    name: central_difference(name, a, b, c)
    for name in ("a", "b", "c")
}

analytical = {
    "a": 3.0,
    "b": 3.0 * c,
    "c": 3.0 * b,
}

print("J:", objective(a, b, c))
print("Analytical:", analytical)
print("Numerical:", {name: round(value, 6) for name, value in numerical.items()})

The output is approximately:

J: 33.0
Analytical: {'a': 3.0, 'b': 6.0, 'c': 9.0}
Numerical: {'a': 3.0, 'b': 6.0, 'c': 9.0}

Finite differences subtract nearby function values, so making the step arbitrarily small can increase rounding error. Avoid nondifferentiable points and compare more than one step size when results disagree. A relative-difference measure is often useful:

\[\text{difference}=\frac{\lVert g_{\text{analytical}}-g_{\text{numerical}}\rVert_2}{\lVert g_{\text{analytical}}\rVert_2+\lVert g_{\text{numerical}}\rVert_2}\]

If both gradient norms are zero, the displayed ratio is undefined; handle that case separately and report zero absolute difference. When both gradients are very small, relative error can be large even when absolute error is tiny. Report both, and choose tolerances for the scale, numerical precision, and test points. Agreement supports the implementation at those points, rather than proving it correct for every input.

For gradient verification and numerical-stability practices, see Gradient Checking and a Systematic Debugging Method for Neural Networks.

A compact forward and backward implementation

def forward(a, b, c):
    u = b * c
    v = a + u
    J = 3 * v
    cache = {"b": b, "c": c}
    return J, cache


def backward(cache):
    # Seed gradient
    dJ = 1.0

    # J = 3v
    dv = dJ * 3.0

    # v = a + u
    da = dv * 1.0
    du = dv * 1.0

    # u = bc
    db = du * cache["c"]
    dc = du * cache["b"]
    return {"a": da, "b": db, "c": dc}


J, cache = forward(a=5.0, b=3.0, c=2.0)
gradients = backward(cache)

print("J:", J)
print("Gradients:", gradients)

J: 33.0
Gradients: {'a': 3.0, 'b': 6.0, 'c': 9.0}

The cache makes the dependency between the passes explicit. The multiplication node needs the original values of \(b\) and \(c\), so the forward pass stores them for reuse.

The earlier numerical check compared finite differences with formulas entered by hand. The following check calls backward() itself at several inputs. Distinct values and a zero input help expose errors that equal inputs might hide:

for values in ((5.0, 3.0, 2.0), (-1.0, 0.0, 4.0), (2.0, -3.0, 0.5)):
    _, saved = forward(*values)
    actual = backward(saved)
    numerical = {
        name: central_difference(name, *values)
        for name in ("a", "b", "c")
    }
    matches = all(
        np.isclose(actual[name], numerical[name], rtol=1e-7, atol=1e-8)
        for name in actual
    )
    print(values, matches)
# (5.0, 3.0, 2.0) True
# (-1.0, 0.0, 4.0) True
# (2.0, -3.0, 0.5) True

All three checks agree within the stated tolerances. They test the implemented backward rules at these inputs; the derivation explains why those rules apply to the graph.

Why reverse-mode differentiation is efficient

Machine-learning models usually have many parameters but one scalar cost. Reverse mode calculates the derivative of that one output with respect to all upstream parameters in a single backward traversal, reusing derivatives that have already been computed.

Calculating every parameter derivative independently would repeatedly traverse the same later operations. Reverse mode avoids that duplication. This makes it particularly well suited to neural-network training.

Connection to logistic regression

Logistic regression has a longer computation graph:

\[(w,b,x)\rightarrow z=w^Tx+b\rightarrow a=\sigma(z)\rightarrow\mathcal{L}(a,y)\]

The forward pass calculates the score, probability, and loss. The backward pass applies local derivatives in reverse order. For the per-example binary cross-entropy \(\mathcal{L}=-y\log a-(1-y)\log(1-a)\) with \(a=\sigma(z)\), multiplying \(\partial\mathcal{L}/\partial a\) by \(a(1-a)\) gives:

\[\frac{\partial\mathcal{L}}{\partial z}=a-y\]

For a mean over \(m\) examples, each score receives \((a-y)/m\). That averaging factor belongs in the backward calculation once. The resulting parameter gradients are then used by the optimizer. The full derivation and NumPy training implementation are in Logistic Regression: A Complete Guide to Binary Classification.

Connection to neural networks

The feedforward networks considered here repeat affine transformations and activation functions:

\[X\rightarrow Z^{[1]}\rightarrow A^{[1]}\rightarrow Z^{[2]}\rightarrow A^{[2]}\rightarrow\cdots\rightarrow J\]

The forward pass computes and caches layer values. The backward pass begins at the cost and calculates \(dW^{[l]}\), \(db^{[l]}\), and the upstream gradient for the preceding layer. For an elementwise activation, the local derivatives multiply the incoming gradient elementwise. A matrix operation instead combines contributions with matrix products: for \(Z=WA\), the affine part passes \(W^T dZ\) toward \(A\). Linear Algebra for Deep Learning derives that expression and the parameter gradients. A shape check is useful, but it cannot replace the chain rule or a value check.

Review Neural Networks for Beginners: From Architecture to Forward Propagation for the corresponding forward computations and matrix shapes.

A note about gradient descent

If variables are trainable parameters of a meaningful cost, gradient descent updates each parameter in the direction opposite its gradient. However, \(J=3(a+bc)\) is used here only to demonstrate differentiation. It is not a practical machine-learning loss and is not bounded below, so repeatedly minimizing it would not define a useful learning problem.

In logistic regression and neural networks, \(J\) is instead a carefully chosen cost such as cross-entropy. Backpropagation calculates its gradients; the optimizer decides how to update the parameters.

The autograd mechanics that replace this hand-written derivative are in PyTorch Tensors and Autograd for People Who Wrote Backprop by Hand.

Exercises

1. By hand. For \(f(a,b,c)=(a+b)\cdot\max(c,0)\) at \(a=2,b=3,c=4\): draw the graph, run the forward pass, and compute \(\partial f/\partial a\), \(\partial f/\partial b\), \(\partial f/\partial c\) by hand. Then change \(c\) to \(-1\) and redo it. Explain in one sentence why the gradients change the way they do.

You should get: three small integers for \(c=4\), with two of them equal. For \(c=-1\), all three become the same value.

# Check your answer:
def f(a, b, c): return (a + b) * max(c, 0.0)
h = 1e-6
for c in (4.0, -1.0):
    grads = []
    for i in range(3):
        p, m = [2.0, 3.0, c], [2.0, 3.0, c]
        p[i] += h; m[i] -= h
        grads.append(round((f(*p) - f(*m)) / (2 * h), 6))
    print(f"c={c}: da, db, dc =", grads)
# c=4.0: da, db, dc = [4.0, 4.0, 5.0]
# c=-1.0: da, db, dc = [0.0, 0.0, 0.0]

2. Multiple paths. Let \(y=x\cdot x + 3x\). Compute \(dy/dx\) at \(x=4\) by treating the two occurrences of \(x\) in the product as separate inputs and summing every contribution. Verify against the numerical derivative.

You should get: a two-digit number, and three separate contributions that add up to it.

3. Implement a node. Write sigmoid forward and backward functions using a cache, following the implementation above. Its local derivative is \(\sigma(z)(1-\sigma(z))\), and the backward method should reuse the cached output rather than recomputing \(\sigma(z)\). Check the derivative of \(J(z)=2\sigma(z)\) numerically, so the incoming gradient is 2 rather than 1.

You should get: close agreement at the moderate inputs tested below. At a saturated input, inspect absolute error as well as relative error before diagnosing a bug.

import numpy as np

def sigmoid_forward(z, cache):
    a = 1.0 / (1.0 + np.exp(-z))
    cache["a"] = a
    return a

def sigmoid_backward(da, cache):
    a = cache["a"]
    return da * a * (1.0 - a)

def relative_error(a, b):
    denominator = abs(a) + abs(b)
    return abs(a - b) / denominator if denominator else 0.0

def objective_sigmoid(z):
    return 2.0 * sigmoid_forward(z, {})

h = 1e-6
for z in (-2.0, 0.0, 2.0, 30.0):
    cache = {}
    sigmoid_forward(z, cache)
    analytical = sigmoid_backward(2.0, cache)
    numerical = (objective_sigmoid(z+h) - objective_sigmoid(z-h)) / (2*h)
    print(f"z={z:4.0f} abs={abs(analytical-numerical):.2e} "
          f"rel={relative_error(analytical, numerical):.2e}")

cache = {}
sigmoid_forward(2.0, cache)
local = cache["a"] * (1.0 - cache["a"])
for da in (1.0, 2.0):
    correct = sigmoid_backward(da, cache)
    print(f"omitted upstream, da={da:.0f}: rel={relative_error(correct, local):.6f}")
# z=  -2 abs=1.13e-11 rel=2.69e-11
# z=   0 abs=6.99e-11 rel=6.99e-11
# z=   2 abs=4.42e-11 rel=1.05e-10
# z=  30 abs=1.87e-13 rel=1.00e+00
# omitted upstream, da=1: rel=0.000000
# omitted upstream, da=2: rel=0.333333

4. Break it deliberately. For \(J(x,y)=xy\), deliberately return (x, x) as the gradient instead of (y, x). Compare it with central differences at equal and distinct inputs, then test one update at \((x,y)=(1,-2)\).

You should get: a missed bug when the inputs are equal, a visible discrepancy at distinct inputs, and an increase in this objective after the deliberately incorrect update.

import numpy as np

def product(x, y):
    return x * y

h = 1e-6
for x, y in ((2.0, 2.0), (1.0, -2.0)):
    numerical = np.array([
        (product(x+h, y) - product(x-h, y)) / (2*h),
        (product(x, y+h) - product(x, y-h)) / (2*h),
    ])
    wrong = np.array([x, x])
    denominator = np.linalg.norm(wrong) + np.linalg.norm(numerical)
    relative = np.linalg.norm(wrong-numerical) / denominator
    print(f"x={x:.0f}, y={y:.0f}: relative error {relative:.6f}")

x, y, learning_rate = 1.0, -2.0, 0.01
new_x, new_y = x-learning_rate*x, y-learning_rate*x
print(f"objective: {product(x,y):.4f} -> {product(new_x,new_y):.4f}")
# x=2, y=2: relative error 0.000000
# x=1, y=-2: relative error 0.821854
# objective: -2.0000 -> -1.9899

5. Reason about cost. A network has one scalar loss and \(10^{6}\) parameters. Estimate the number of forward passes needed to obtain all gradients by central differences, and compare with reverse-mode differentiation. Assume one backward pass costs about two forward passes. Use the ratio to explain why coordinate-by-coordinate central differences are impractical for training a model of this size.

You should get: a ratio around \(10^{6}\) — enough that a training step measured in seconds becomes one measured in days.