From a slope to backpropagation: build a tiny autograd engine

In our article on how neural networks learn, we followed a model through its training loop: make a prediction, measure its error with a loss function, calculate gradients, and update the weights. To adjust the weights, we used gradients—derivatives that tell us how the loss changes with respect to each weight. Gradient descent uses these derivatives to choose the direction of an update.

In this article, we want to focus on how those gradients are calculated. If a model’s loss is $9$, how do we work out that increasing one weight would reduce it? Which function are we differentiating, and where do its derivative rules come from?

To answer that, we need to follow the computation that produced the loss. The prediction and loss together form a function of the model’s parameters. We’ll build up from $f(x)=x^2$, calculate derivatives through individual operations, and implement a small automatic differentiation engine in Python. Then we’ll use it to train a model.

Our engine uses reverse-mode automatic differentiation, or autodiff: it records the operations in a forward evaluation at the current parameter values, then applies derivative rules backward. We will also use finite differences for two separate purposes: visualizing what a slope means and independently checking the engine’s answers. Those nearby-point estimates belong to the explanation and verification; the training loop gets its gradients from autodiff.

You only need basic Python and algebra. We will introduce the calculus as we use it. The implementation follows the scalar-computation approach taught in Andrej Karpathy’s micrograd lecture, using an explicit table of local derivatives on graph edges. It requires no machine-learning libraries.

Start with a function and its derivative

Consider:

At $x=5$, the output is $25$. The power rule gives its derivative, written $f'(x)$ or $df/dx$:

The derivative is the slope of the tangent at that point. It tells us how sensitive the output is to a small input change:

The derivative is the quantity we want. Several methods can calculate or approximate that quantity.

Three ways to obtain the same derivative

Keep the function $f(x)=x^2$ and evaluation point $x=5$ fixed:

Symbolic differentiation produces a derivative expression that we can evaluate later. Autodiff evaluates derivative rules using the values recorded during the computation. For this smooth function, both calculate the derivative at the chosen point, subject to floating-point arithmetic in code. Finite differences estimate it by comparing outputs at different points.

Neural-network backpropagation uses the third approach. We supply derivative rules for basic operations, and the engine combines them through the chain rule. It can compute the gradient at the current parameter values without generating one large symbolic expression for the whole network.

Autodiff through one multiplication

Write the same function as a program:

def f(x):
    return x * x

We already know the derivative rule for multiplication. If $z=ab$, its local derivatives are:

For x * x, both operands refer to the same input. At $x=5$, each operand is $5$, so each local derivative is also $5$. A change in $x$ affects the result through both operand paths. Their contributions add:

Select Evaluate x × x, then Seed output gradient. Follow each operand in turn to see the accumulated gradient reach $10$.

The widget starts the backward pass with the output’s derivative with respect to itself: $dq/dq=1$. Each path multiplies that incoming gradient by its local derivative. In ordinary Python, the arithmetic is:

x = 5.0
left, right = x, x
q = left * right             # forward: 25.0

output_grad = 1.0            # dq/dq
x_grad = 0.0
x_grad += output_grad * right  # left operand contributes 5.0
x_grad += output_grad * left   # right operand contributes 5.0
print(x_grad)                 # 10.0

Every value needed for this derivative came from the computation at $x=5$. We used the multiplication rule and the saved operands. Evaluating $f(5.001)$ belongs to the finite-difference method in the table; the autodiff trace above is complete with the single evaluation point.

The automation comes from recording these operand paths and applying their derivative rules for us. Later, our Value class will make the same calculation happen when we write q = x * x and call q.backward().

Knowing the operations and understanding the shape

Suppose we had only a table:

Those finite observations do not uniquely determine the function between the listed points, or its exact derivative at a listed point. Multiple functions can agree with every row and have different slopes. Computing a finite-difference estimate from available nearby samples is one option.

With the program return x * x, we have more information: we know the operations that generate the output, so we can apply their derivative rules. A neural network likewise has known computations for its predictions and loss. The training examples may arrive as a table, but the model and loss are executable mathematical functions.

Knowing those operations also does not automatically give us an understanding of the entire curve. Consider:

We can evaluate the value and derivative at any chosen point using these formulas. Can you picture every rise, dip, and turning point just from reading them? Move Known evaluation point x below, then select Reveal full curve. The local value and slope are available before the full plot is visible.

The formula specifies the function completely. Understanding its global geometry requires additional evaluation or analysis. A plot gives us more of that information over a chosen interval, while a derivative describes local sensitivity at one point.

For a neural network, the loss may depend on millions of weights and biases. We know the program that computes $L(w_1,w_2,\ldots)$ at a chosen parameter setting. Its gradient gives local sensitivities there. Those computations do not locate every minimum elsewhere in that high-dimensional space. Training uses the local information to choose an update, then calculates a new gradient at the updated parameters.

Use nearby points to visualize a slope

Now return to $f(x)=x^2$ for the first of our two uses of finite differences: building slope intuition. We have already calculated its derivative through autodiff. This experiment illustrates what that derivative means geometrically.

Compare $f(5)=25$ and $f(5.001)=25.010001$. Their average slope is:

For a general nonzero input change $h$:

As $h$ approaches zero, this ratio approaches the derivative $2x$. At $x=5$, the finite-difference estimate approaches $10$.

Reduce Small change h to $0.001$ and watch the orange secant approach the purple tangent. Then move Position x to a negative value: a negative slope means that a small increase in the input decreases the output.

This nearby-point calculation is a finite difference. Later, we will use another finite difference to independently verify the engine. In the training experiments, gradients come from recorded operations and their local derivative rules.

For $x^2$, the exact output change is $2x\Delta x+(\Delta x)^2$. The derivative captures the first-order part; it does not predict arbitrary finite changes exactly.

A prediction gives us something to evaluate

Now suppose a model receives the input $x=2$ and should predict the target $y=5$. Give it one weight $w$ and one bias $b$:

The hat in $\hat y$ marks a prediction. At $w=1$ and $b=0$, the prediction is $2$, which is $3$ below the target.

What should training minimize? Minimizing the prediction itself would drive it downward, even though we need it to increase. We need a function that measures the discrepancy between prediction and target.

Choose squared error:

The current loss is $(2-5)^2=9$. Zero would mean the prediction matches the target for this example.

Substitute the model’s computation into the loss:

This is the connection we need. The loss is now a function of the parameters we can change. For the fixed training example $(x,y)=(2,5)$:

The number $9$ alone does not tell us which way to move. With $b=0$, both $w=1$ and $w=4$ produce loss $9$, but they need to move in opposite directions. We differentiate the function that produced the loss, evaluated at the current parameters.

The same calculus, with a weight as the variable

Hold the training input $x=2$, target $y=5$, and bias $b=0$ fixed. The loss becomes an ordinary function of one variable:

We can differentiate this small expression by hand and substitute the current weight, exactly as we did for $x^2$:

The letter naming the variable has changed from $t$ to $w$. The training input and target are fixed numbers inside $L(w)$; we are asking how the loss responds to the weight. The same distinction applies when we use a whole batch: the batch’s inputs and targets stay fixed while we differentiate with respect to the parameters.

This model is small enough to expand into a polynomial. For a larger network, autodiff records the forward operations at the current parameter setting and evaluates their local derivative rules backward. It computes the numerical values of the parameter derivatives without first constructing one expanded symbolic formula for the network’s gradient.

Which derivative do we take?

For the weight, we want:

Read this as: how quickly would the loss change if we changed the weight, holding the other independent quantities fixed? The symbol $\partial$ is used because $L$ depends on several variables. For this calculation, the input, target, and bias stay fixed.

The model and loss form a composition. First compute the prediction, then subtract the target, then square the error:

Each operation has a local derivative:

The chain rule multiplies these sensitivities along the path:

At our starting values, $e=-3$ and $x=2$:

A weight increase of $0.01$ should therefore change the loss by approximately $-12\times0.01=-0.12$. The exact new loss is $(2.02-5)^2=8.8804$, a change of $-0.1196$. The agreement improves as the step gets smaller.

For the bias, the last factor changes: $\partial\hat y/\partial b=1$. So:

The derivatives with respect to all trainable parameters together form the gradient. Here that vector is $\nabla L=(-12,-6)$, in the order $(w,b)$.

The choice of loss changes the learning signal

The model tells us how its prediction responds to a weight. The loss tells us how much that prediction change matters to the objective.

For the same prediction, compare two losses:

At our starting point, the prediction gradient is $-6$ for squared error and $-1$ for absolute error. Multiplying by the same model derivative, $x=2$, gives weight gradients $-12$ and $-2$.

Switch between Squared error and Absolute error below. Move Weight w while watching both the prediction gradient and the weight gradient. The input stays fixed, so only the first factor changes.

Squared error produces a larger prediction gradient when the error’s magnitude is larger. Absolute error has constant gradient magnitude with respect to the prediction whenever the error is nonzero. Its derivative is undefined at zero; the widget selects a subgradient of zero there.

Neither loss supplies the rules for differentiation. Those rules come from calculus. The choice of loss determines the function those rules act on, and therefore the initial signal propagated through the model.

A derivative becomes an update

Gradient descent subtracts a scaled gradient:

The learning rate $\eta$ sets the step size. With squared error, $w=1$, and $\eta=0.1$:

Holding the bias at zero, the new prediction is $4.4$ and the new loss is $0.36$. Use Reset and Take one step in the loss experiment to reproduce this calculation.

Here we updated only the weight. If we also update the bias from the same gradient, it becomes $0.6$, giving prediction $2.2\times2+0.6=5$ and loss zero for this single example. These are different updates, so they produce different results.

A negative derivative suggests increasing the parameter; a positive derivative suggests decreasing it. That is a local statement. A large step can pass a minimum and increase the loss. Backpropagation calculates the gradient; the optimizer and its learning rate determine what to do with it.

Backpropagation is the chain rule through a graph

We can record the computation as a graph. Each box holds a number, and each operation connects its result to the values used to produce it.

Start at the rightmost box. The derivative of the loss with respect to itself is one:

Then work backward. A box receives a gradient telling us how the final loss responds to that box’s value. To pass a contribution to one of its inputs, multiply by the operation’s local derivative:

Select Seed gradient, then Backpropagate one operation until you reach $w$ and $x$. Every gradient displayed is a derivative of the final loss, even when the box represents an intermediate result.

The graph also gives derivatives with respect to the input and target. That does not make them trainable. We choose which values the optimizer updates: here, only $w$ and $b$. In a deeper network, a layer’s “input” may be a previous layer’s output, so its gradient becomes the route back to earlier parameters.

An important distinction appears when a value has multiple paths to the loss. Contributions from those paths add. For example:

At $a=3$, that derivative is $7$. A backward pass must keep all three contributions, including both uses of $a$ in the multiplication. This is why our engine will use += when propagating gradients.

Build a number that remembers its computation

An ordinary Python float stores a value. Our Value will also store:

• grad: the derivative of the chosen final output with respect to this value.
• edges: a sequence of (parent, local_derivative) pairs recording how this value was computed.

Addition and multiplication become graph-building operations. For $z=a+b$, both edge derivatives are $1$. For $z=ab$, they are $b$ and $a$.

Here is the first part of the engine, including integer powers and the activation we will use later:

import math


class Value:
    def __init__(self, data, edges=()):
        self.data = float(data)
        self.grad = 0.0
        self.edges = tuple(edges)  # (parent, local derivative) pairs

    @staticmethod
    def wrap(value):
        return value if isinstance(value, Value) else Value(value)

    def __add__(self, other):
        other = Value.wrap(other)
        return Value(self.data + other.data, ((self, 1.0), (other, 1.0)))

    def __mul__(self, other):
        other = Value.wrap(other)
        return Value(self.data * other.data,
                     ((self, other.data), (other, self.data)))

    def __pow__(self, exponent):
        if not isinstance(exponent, int) or exponent < 1:
            raise ValueError("This engine supports positive integer powers only")
        return Value(self.data ** exponent,
                     ((self, exponent * self.data ** (exponent - 1)),))

    def tanh(self):
        output = math.tanh(self.data)
        return Value(output, ((self, 1.0 - output * output),))

When we multiply two Value objects, __mul__ calculates the output immediately and saves two edges with their local derivatives. These derivative values are evaluated using the operands’ current data. Nothing is learned yet; we are recording a forward computation.

We implement powers directly, so (prediction - target) ** 2 records a square operation. The engine intentionally supports only positive integer powers. This covers our examples without introducing domains such as fractional powers of negative numbers.

The hyperbolic tangent, tanh, is a smooth activation function with derivative $1-\tanh^2(z)$. It lets a network express nonlinear relationships.

Python’s convenience operators make the engine pleasant to use. Reverse operators allow expressions such as 2 * w and sum(...):

    def __neg__(self):
        return self * -1.0

    def __sub__(self, other):
        return self + -Value.wrap(other)

    def __rsub__(self, other):
        return Value.wrap(other) + -self

    __radd__ = __add__
    __rmul__ = __mul__

These methods all belong inside Value. The complete Python file includes the engine, the models, and runnable examples. You can also download tinygrad.py.

Visit the graph in the right order

To calculate gradients, we need to process each result before its parents, after all of that result’s downstream contributions have arrived.

First, a depth-first traversal visits a node’s parents before appending the node itself. This creates a topological ordering. A seen set ensures a shared node is appended only once. Reversing the order puts the final output first.

Add this method to Value:

    def backward(self):
        order, seen = [], set()

        def visit(node):
            if node in seen:
                return
            seen.add(node)
            for parent, _ in node.edges:
                visit(parent)
            order.append(node)

        visit(self)
        for node in order:
            node.grad = 0.0
        self.grad = 1.0
        for node in reversed(order):
            for parent, local_derivative in node.edges:
                parent.grad += node.grad * local_derivative

There are three separate steps: clear old gradients in the reachable graph, seed the final output’s gradient with one, and traverse backward to accumulate contributions.

The seen set deduplicates nodes, while edges preserves every operand occurrence. In a * a, two edges point to the same parent. Both execute their += update. Replacing edges with a set of parents without retaining the corresponding contributions would lose that distinction.

Our backward() deliberately clears reachable gradients on each call. This makes repeated calls on the same unchanged graph reproducible. It differs from APIs that accumulate leaf gradients across backward calls and require a separate zeroing step. Within one backward pass, contributions still accumulate from all graph paths.

This implementation saves numeric local derivatives during the forward pass. Change parameter data after backward, then build a fresh graph for the next iteration. Rerunning backward on an old graph after changing data would use stale local derivatives.

The recursive graph traversal keeps this teaching engine short. Long computation chains, such as a loss summed over thousands of examples, can exceed Python’s recursion limit; an iterative traversal avoids that limit.

Check the engine against the arithmetic

Save the complete file as tinygrad.py. In another Python file in the same directory:

from tinygrad import Value

w, b = Value(1), Value(0)
prediction = w * 2 + b
loss = (prediction - 5) ** 2
loss.backward()

print(loss.data)  # 9.0
print(w.grad)     # -12.0
print(b.grad)     # -6.0

a = Value(3)
q = a * a + a
q.backward()
print(a.grad)    # 7.0: all three paths contribute

The engine never received a special formula for this model’s weight gradient. It recorded multiplication, addition, subtraction, and a square, then combined their local derivatives.

Use finite differences to check the result

This is the second purpose announced at the start: independent verification. The engine has already obtained w.grad from the recorded operations at $w=1$. To check its implementation, evaluate the loss separately on both sides of that weight:

def numeric_loss(weight):
    return (weight * 2 - 5) ** 2

h = 1e-5
finite_difference = (numeric_loss(1 + h) - numeric_loss(1 - h)) / (2 * h)
assert abs(w.grad - finite_difference) < 1e-6

This central difference is a numerical check, not the engine’s differentiation algorithm. Checking every parameter this way needs two evaluations per parameter. Reverse-mode automatic differentiation instead reuses a single recorded graph to compute all parameter gradients for a scalar output, with work proportional to the graph’s size.

Extremely small values of $h$ eventually introduce floating-point cancellation, so smaller is not always better for this numerical check. The training loops below obtain their gradients by calling backward() on the current graph.

Learn from several examples

One training example cannot determine a unique line: many $(w,b)$ pairs satisfy $2w+b=5$. Add more examples:

These points lie on $y=2x+1$, so we can check whether training finds $w\approx2$ and $b\approx1$.

Define the loss as the mean of the four squared errors:

The shared weight appears in all four predictions. During backward, each path contributes to the same w.grad. Multiplication by $1/4$ scales all contributions to match the mean loss. Summing instead of averaging would multiply both gradients by four.

Use Step to inspect one update, then Run 100 steps. Try learning rate $0.3$ from the reset state to see loss grow even though the gradients are calculated correctly.

The widget uses a JavaScript counterpart of the scalar engine above. Its gradients are calculated through a new graph on each step. Here is the same training loop in Python:

def train_line(steps=1000, learning_rate=0.05):
    data = [(0, 1), (1, 3), (2, 5), (3, 7)]
    w, b = Value(-1), Value(0)
    for _ in range(steps):
        loss = sum((w * x + b - y) ** 2 for x, y in data) * (1 / len(data))
        loss.backward()  # also clears this graph's old gradients
        for parameter in (w, b):
            parameter.data -= learning_rate * parameter.grad
    final_loss = sum((w.data * x + b.data - y) ** 2 for x, y in data) / len(data)
    return w.data, b.data, final_loss

Both gradients come from the same forward pass before either parameter changes. The next iteration recomputes predictions, loss, and derivatives at the updated parameters. The returned loss is evaluated after the last update, so it describes the returned weights.

Use the same engine for a neural network

So far, the model has been affine. A small neural network needs no new differentiation machinery: its forward pass is a larger composition of the operations we already support.

A neuron multiplies inputs by weights, adds a bias, and optionally applies tanh. An MLP chains layers of neurons. We keep the final layer linear for this regression example:

import random


class Neuron:
    def __init__(self, n_inputs, rng, nonlinear=True):
        self.weights = [Value(rng.uniform(-1, 1)) for _ in range(n_inputs)]
        self.bias = Value(0)
        self.nonlinear = nonlinear

    def __call__(self, inputs):
        total = sum((w * x for w, x in zip(self.weights, inputs)), self.bias)
        return total.tanh() if self.nonlinear else total

    def parameters(self):
        return self.weights + [self.bias]


class MLP:
    def __init__(self, sizes, seed=42):
        rng = random.Random(seed)
        self.layers = [
            [Neuron(n_in, rng, nonlinear=i < len(sizes) - 2)
             for _ in range(n_out)]
            for i, (n_in, n_out) in enumerate(zip(sizes, sizes[1:]))
        ]

    def __call__(self, inputs):
        for layer in self.layers:
            inputs = [neuron(inputs) for neuron in layer]
        return inputs

    def parameters(self):
        return [p for layer in self.layers for neuron in layer
                for p in neuron.parameters()]

MLP([2, 4, 1]) creates two inputs, four hidden neurons with tanh, and one linear output. Its parameters are 8 hidden weights, 4 hidden biases, 4 output weights, and 1 output bias: 17 trainable scalars.

Train it on XOR, encoding equal inputs as $-1$ and different inputs as $+1$. No affine function of two inputs can fit all four targets exactly; the hidden nonlinear activations allow the model to fit them.

def train_xor(steps=1000, learning_rate=0.1):
    # Same inputs twice equal -1; different inputs equal +1.
    examples = [([0, 0], -1), ([0, 1], 1), ([1, 0], 1), ([1, 1], -1)]
    model = MLP([2, 4, 1], seed=42)
    parameters = model.parameters()
    for _ in range(steps):
        loss = sum((model(x)[0] - y) ** 2 for x, y in examples) * 0.25
        loss.backward()
        for parameter in parameters:
            parameter.data -= learning_rate * parameter.grad
    predictions = [model(x)[0].data for x, _ in examples]
    final_loss = sum((p - y) ** 2 for p, (_, y) in zip(predictions, examples)) / 4
    return predictions, final_loss

The training structure is unchanged: build predictions and a scalar loss, call backward(), and update parameters. The engine sees operations and edges. It does not need to know which operations belong to a neuron or a layer.

Run the complete file:

python3 tinygrad.py

With the supplied seed, 1,000 steps, and learning rate $0.1$, the XOR predictions round to [-1.0, 1.0, 1.0, -1.0] at four decimal places. This is a small demonstration with fixed data and initialization, rather than a guarantee for arbitrary networks or learning rates.

What the loss contributed

We started with a derivative of $x^2$, which needed no notion of a model or training. Then we chose a prediction function and attached a loss. Substitution made the loss a function of the parameters, giving us a precise question to ask: how would this objective change if one parameter changed?

The loss defines the objective and its sensitivity to the prediction. The model supplies the paths from parameters to that prediction. Backpropagation combines local derivatives along those paths, adding contributions where paths meet. Gradient descent uses the resulting parameter gradients to choose an update.

For a model built from operations with usable derivative rules, a scalar loss gives our differentiation engine an output to start from, and an objective for the optimizer to reduce.

For the broader neural-network treatment, see How neural networks learn: backpropagation and gradient descent. For a complementary implementation that attaches backward functions to operations, see Karpathy’s micrograd engine and his Neural Networks: Zero to Hero materials.