From derivatives to backpropagation: build a tiny autograd engine

In our earlier article on how neural networks learn, we walked through the basic training loop: the model makes a prediction, we measure how wrong it is using a loss function, we tweak the weights, and then we repeat. To adjust each weight correctly, we need to know whether making it bigger or smaller will reduce the loss. One simple way to figure this out is to nudge a weight slightly up or down, recalculate the loss, and see if it improved. The problem is that real models can have millions or billions of these weights. Testing every single one individually would be incredibly slow.

This is where derivatives come in. A derivative describes how quickly a function’s output changes as we change an input, at a particular input value. In this article, we’ll mostly use the loss as our function and examine how it changes with the model’s parameters, keeping the training data fixed. The same idea applies to other functions, including the individual operations inside the model.

For each weight, the derivative tells us how sensitive the loss is to a small change in that weight, assuming all the other weights stay the same. This gives us a direction without testing each possible change: a positive derivative tells us to decrease the weight to reduce the loss locally; a negative derivative tells us to increase it.

Here’s how the ideas fit together:

• Finite differences: estimate a derivative using two function evaluations at nearby inputs.
• Limit definition: mathematically define the exact derivative as the interval shrinks toward zero.
• Symbolic differentiation: use analytical derivative rules to derive a formula for the derivative.
• Backpropagation (reverse-mode autodiff): combine numerical local derivatives backward through a computation graph to calculate gradients.

We’ll build a small Python engine for reverse-mode automatic differentiation, check its gradients with finite differences, and use it to train a model.

A familiar rate of change

Consider the relationship $y=3x$. If $x$ increases by $1$, $y$ increases by $3$. If $x$ increases by $0.1$, $y$ increases by $0.3$. Whatever the starting value, the change in $y$ is three times the change in $x$:

The symbol $\Delta$ means a change between two values. You can recognize this rate of change without knowing calculus. The derivative gives it a name and notation:

We read $dy/dx$ as the derivative of $y$ with respect to $x$. The denominator tells us which variable we change: $x$. The numerator tells us whose response we measure: $y$. So $dy/dx=3$ means that changing $x$ by a small amount changes $y$ by three times that amount.

The function $y=3x$ tells us the output value for a given input. Its derivative tells us how quickly that output changes as the input changes. Here, the rate is always $3$, so $\Delta y=3\,\Delta x$ is exact, even for large changes. The graph is a straight line with the same slope everywhere.

For $y=x^2$, the derivative is $2x$, which we’ll derive shortly. At $x=1$, the rate is $2$; at $x=5$, it is $10$. Because the rate changes as we move, using the slope at the starting point to predict an output change generally gives only a local approximation. Let’s explore this distinction using a loss function whose slope changes with its input.

Finding a direction by testing nearby values

Suppose our loss depends on one parameter $a$:

At $a=5$, the loss is $28$. Should we increase or decrease $a$ to reduce it?

We could try both directions:

Decreasing $a$ by $0.1$, from $5$ to $4.9$, reduced the loss from $28$ to $27.01$. Increasing $a$ by $0.1$, from $5$ to $5.1$, raised the loss from $28$ to $29.01$. These evaluations suggest that a small decrease from $a=5$ will reduce the loss.

The direction that reduces the loss depends on the current value of $a$. If we started at $a=-5$, the loss would also be $28$, but increasing $a$ to $-4.9$ would lower it to $27.01$. Decreasing $a$ to $-5.1$ would raise it to $29.01$. For this function, moving toward zero reduces the loss: from a positive value we should decrease $a$, and from a negative value we should increase it.

This is why we describe the direction as local: it tells us which way to move from the current parameter value for a small change. The direction can change as we move, and the size of the change matters too.

Let’s now use the same values to estimate how quickly the loss changes as we change the parameter. We can do this using a method called finite differences. A finite-difference estimate is a numerical approximation of a function’s derivative, or rate of change, obtained by dividing the difference between its values at two nearby points by the distance between those points.

There are three common types. The difference is which points we choose relative to the point where we want to estimate the derivative:

1. Forward difference: compare the current point with a point ahead of it.
2. Backward difference: compare the current point with a point behind it.
3. Central difference: compare two points equally spaced on opposite sides of the current point.

For example, we could estimate a car’s current speed by comparing its distance reading now with the reading a moment ago and dividing by the elapsed time. This is a backward difference: we use a past reading because a future one isn’t available yet.

We’ll use forward differences for our main walkthrough. Let’s start by estimating the slope of our loss at $a=5$, then compare the result with backward and central differences.

Writing the parameter change as $h$, the forward-difference formula is:

With $h=0.1$, we compare the loss at the current value $a=5$ with the loss at $a+h=5.1$:

Our forward-difference estimate of $10.1$ approximates the slope at $a=5$. It uses the table’s two loss values, $L(5)=28$ and $L(5.1)=29.01$. Because the estimate is positive, it indicates that decreasing $a$ should reduce the loss.

This ratio is the average change in loss per unit change in the parameter over that interval: we divide the total loss change by the total parameter change. The calculation uses only the two endpoint values. Those values alone don’t tell us how the function behaves between them: the loss could change more slowly in one part of the interval and more quickly in another. So $10.1$ describes the rate across the whole interval, rather than the exact rate at every point within it.

Here we know the formula, so we can plot the curve between the endpoints and see that distinction. The graph below uses a wider interval, from $a=1$ to $a=5$, to make the difference visible. The dashed line connecting the endpoints has the average slope; the green curve shows how the loss actually changes. The shaded triangle shows the total parameter change $\Delta a=4$ and loss change $\Delta L=24$, so the average slope is $24/4=6$.

We can compare the average rate over the whole interval with the rates over its two halves. Splitting at $a=3$, where the loss is $12$, gives:

The overall average of $6$ summarizes the total change; it doesn’t mean the loss rises at a constant rate throughout the interval.

Now let’s see what backward and central differences produce for the same loss $L(a)=a^2+3$, at the same point $a=5$, using the same step size $h=0.1$. Only the choice of sample points changes.

Notice how central difference divides by $2h$, while forward and backward differences divide by $h$. The denominator is the distance between the two sample points. For central difference, those points are $a-h$ and $a+h$, so:

At $a=5$ and $h=0.1$:

Each point is $0.1$ from the current point, but they are $0.2$ apart. Forward and backward differences use the current point as one endpoint, so their distances are $(a+h)-a=h$ and $a-(a-h)=h$, respectively.

Central difference gives $10$ in this example. We’ll derive that slope directly next, and return to the accuracy of central differences when we use them to check our engine.

With millions of parameters, estimating a slope this way for each one would require millions of additional loss evaluations. To avoid those evaluations, we need a mechanism that gives us the slope directly.

From finite differences to the derivative

Finite differences measure an average rate of change over an interval. To get the slope at the current parameter value, we consider what that rate approaches as the interval shrinks toward zero. For our loss function, we can work out this limiting value algebraically and obtain a formula for the slope that we can evaluate directly, without testing nearby parameter values.

Let’s first see what shrinking the interval does to our forward-difference estimate. At $a=5$, we used $h=0.1$ and got $10.1$. Reducing $h$ to $0.001$ gives:

With $h=0.001$, our estimate is $10.001$. As we shrink $h$ further toward zero, the estimate approaches $10$.

The value we approach is called the derivative. It tells us how quickly a function changes at a particular point. In our case, the function is the loss $L(a)$, and the point is $a=5$. We write the derivative as $L'(a)$ or $dL/da$. Using the forward-difference form we started with, its definition is:

Using our earlier change notation, $\Delta a=h$ and $\Delta L=L(a+h)-L(a)$. The ratio $\Delta L/\Delta a$ is the average rate over an interval; $dL/da$, written with the Latin letter $d$, is its limit as that interval shrinks toward zero.

For our loss function at $a=5$, this limit is $10$, so $L'(5)=10$. For a differentiable function, all three forms—forward, backward, and central—approach the same derivative value as $h$ approaches zero.

For $L(a)=a^2+3$, we can simplify that comparison algebraically. For nonzero $h$:

As $h$ approaches zero, the ratio approaches $2a$, giving us the derivative rule:

In this way, we build a derivative function: $L'(a)=2a$ gives the slope of $L$ at any value of $a$. When we evaluate it at a particular point, such as $a=5$, we get the slope at that point: $L'(5)=10$. We can calculate this directly, without evaluating $L(5+h)$ to discover the slope. Differentiation gives us a rule for calculating the local rate of change.

We now have a formula for the slope. Next, we’ll use it to choose a direction, then compare deriving such a formula with computing the derivative numerically through the function’s operations.

Using the derivative to choose a direction

The derivative tells us how sensitive the loss is to a small parameter change. At $a=5$, the derivative $L'(5)=10$ means that a small change in $a$ produces a change in loss approximately $10$ times as large:

The derivative gives the instantaneous rate of change at a point. Another way to understand it is as the best constant rate for approximating changes near that point. Here, $L'(5)=10$ is exact. The approximation comes from using that same rate throughout a nonzero parameter change, even though the slope changes as we move. Geometrically, this gives the tangent line: the best local linear approximation to the loss curve.

For example, increasing $a$ by $0.01$ means $\Delta a=0.01$. Substituting this into the formula gives:

The positive result means the loss increases by approximately $0.1$. Decreasing $a$ by $0.01$ means $\Delta a=-0.01$:

The negative result means the loss decreases by approximately $0.1$.

We can choose the local direction from the derivative’s sign without testing both possibilities. At $a=-5$, the derivative is $-10$, so a small increase in $a$ reduces the loss. At $a=0$, the derivative is zero; for this particular function, that is its minimum. A zero derivative does not establish a minimum for every function.

Visualizing the slope

We can visualize the difference between a finite-difference estimate and a derivative on the loss curve. Draw a straight line through two points on the curve: this is a secant, and its slope gives the average rate of change between them. Keep the first point fixed and move the second toward it. The secant approaches the tangent at the first point, whose slope is the derivative.

In the widget below:

• The blue curve is the loss function $L(a)=a^2+3$.
• The orange secant shows the finite-difference estimate, using two points on the curve.
• The purple tangent shows the exact derivative at the current value of $a$.

At the initial settings $a=2$ and $h=0.8$, the orange line’s slope is $4.8$, while the purple line’s slope is $4$. Shrinking $h$ brings the orange slope closer to the purple one.

Set Parameter a to $5$ and reduce Small change h to $0.001$ below to watch the orange secant approach the purple tangent. At $a=5$, the slope estimate becomes $10.001$, close to the derivative $10$. Move Parameter a from $5$ to $-5$ to see the direction of improvement reverse.

The derivative describes a local change. At $a=5$, decreasing all the way to $a=-6$ would increase the loss from $28$ to $39$, despite starting in the suggested direction. The size of the update matters too.

Let’s compare the derivative’s approximation with the exact change in loss. Using the original function $L(a)=a^2+3$, we can calculate the exact change directly by subtracting the original loss from the new loss:

The derivative approximation, $\Delta L\approx L'(a)\Delta a=2a\Delta a$, keeps only the first term. The exact change also includes $(\Delta a)^2$, which becomes small relative to the first term for sufficiently small changes when $a\ne0$. For the move from $a=5$ to $a=-6$, we have $\Delta a=-11$: the derivative predicts $10\times(-11)=-110$, but the exact change is $-110+121=11$. That is why the loss rises from $28$ to $39$ despite the initial direction pointing downhill.

Evaluating a derivative at one point tells us the local sensitivity; it does not locate every minimum elsewhere on the curve. For a neural network, the gradient likewise describes the loss near the current parameter values, while finding useful parameter values requires further updates and evaluations.

Three ways to calculate or estimate a derivative

So far, we’ve used two approaches to find the slope of our loss function: estimating it from nearby values and deriving a formula we can evaluate directly. These are called finite differences and symbolic differentiation.

Neural-network engines use a third approach: automatic differentiation, or autodiff. This is the approach we’ll implement in our Python engine.

Symbolic differentiation produces a derivative expression that we can then evaluate. It can use the chain rule without expanding the original function, and software can perform the symbolic work for us. Our autodiff engine instead evaluates local derivative rules using recorded values and accumulates numerical gradients. The distinction is what we compute—an expression or numerical sensitivities—not whether we use calculus rules, expand brackets, or work by hand. Finite differences provide a separate numerical approximation that we can use to check the engine.

A derivative is the quantity we want; these are three ways to estimate or calculate it. All three methods can find the derivative of the same function at the same point.

Let’s compare these three approaches using the same loss function, $L(a)=a^2+3$, at $a=5$:

Autodiff relies on knowing the operations that produced the output. A finite table of input–output observations alone does not uniquely determine the function between those points or its exact derivatives. A model, however, is a program: its operations are known even when its training data arrive as a table.

All three methods also work for functions with many inputs. To see why neural-network engines use autodiff, consider how much larger their calculations are than our example.

A neural network is a large function built from many repeated, connected mathematical operations: multiplication, addition, and activation functions. Its parameters are the weights and biases used in those operations. Results pass through successive layers, creating a nested calculation that eventually produces a prediction. The loss function then compares that prediction with the target.

For a fixed training batch, we can treat this entire prediction-and-loss calculation as one function of the parameters:

Here, the $w_i$ represent all the trainable parameters, including weights and biases. Our example $L(a)=a^2+3$ has one parameter and a few operations; a neural network’s loss can depend on millions or billions of parameters.

The current point is then the current combination of all those parameter values. We want to know how sensitive the loss is to each parameter while the others stay fixed. Each of these rates is a partial derivative; together, they form the gradient.

Estimating this gradient with finite differences would require additional loss evaluations for each parameter. Reverse-mode autodiff computes all these partial derivatives efficiently for a loss that produces a single number: it uses one forward computation followed by a backward pass. The engine combines the individual derivative rules through the chain rule, without generating one large symbolic expression for the whole network. This is the approach used in backpropagation.

To understand how that backward pass works, we’ll first look at the derivative rules for individual operations, then see how the chain rule connects them.

From local derivatives to the chain rule

We can associate mathematical operations with derivative rules that describe how their outputs respond to small changes in their inputs. These rules are the building blocks for calculating how each parameter affects the final loss.

During training, we need each parameter’s individual effect to decide how to update it. Holding the other parameters fixed isolates that effect; changing several at once would mix their effects together. This is how we calculate each derivative, but it doesn’t require us to update only one parameter at a time. We can use the individual derivatives to update all the parameters together.

Earlier, we wrote $dy/dx$ to describe how $y$ changes when we change $x$. Now our operations have two independent inputs, $x$ and $y$, and an output $z$. To measure how $z$ changes when we vary $x$ while holding $y$ fixed, we use a partial derivative, written $\partial z/\partial x$. We read this as 'the partial derivative of $z$ with respect to $x$'. The order means the same thing as before; the symbol $\partial$ indicates that we’re varying one input while holding the others fixed.

It doesn’t mean we ignore the value of $y$; it means we hold that value fixed while changing $x$. For example, if $z=xy$, the derivative with respect to $x$ is $y$, because $y$ is a constant multiplier of $x$. Its value still determines how strongly a change in $x$ affects $z$.

Let’s see how changing $x$ affects addition and multiplication at $x=3$ and $y=4$, with $y$ held fixed. The table shows each operation’s derivative rule, its value at these inputs, and a finite-difference check that increases only $x$ by $0.001$:

Now we reverse the roles: vary $y$ while keeping $x$ fixed. This gives the partial derivative with respect to $y$, written $\partial z/\partial y$. The checks now increase only $y$ by $0.001$:

The tables above demonstrate that for addition, changing either input changes the output by the same amount, so both partial derivatives are $1$. For multiplication, the change is scaled by the other input: $\partial z/\partial x=y$ and $\partial z/\partial y=x$.

The finite-difference checks recalculate $z$ after changing each input separately. The derivative rules need no such extra evaluations: using the existing inputs $x=3$ and $y=4$, we obtain both multiplication derivatives directly, $\partial z/\partial x=4$ and $\partial z/\partial y=3$.

The derivative rules for addition and multiplication shown in the tables are called local because each describes how one operation’s output changes when one of its inputs changes. To follow an input’s effect through a longer calculation, we need to account for how those operations are connected.

Rules for differentiating combined functions

We’ve examined addition and multiplication with independent inputs. But an operation’s inputs can themselves be functions of the same variable. There are three differentiation rules we’ll use to handle these combinations. Let $f$ and $g$ be differentiable functions of $x$:

The sum rule adds the rates of change of the two terms. It gives another way to differentiate $x^2+3$: the derivative of $x^2$ is $2x$, and the constant $3$ contributes zero.

The product rule accounts for changes in both factors. One contribution comes from the first factor changing, scaled by the current value of the second; the other comes from the second factor changing, scaled by the current value of the first. Multiplying two functions does not mean multiplying their derivatives.

The chain rule applies to function composition: one function’s output becomes another function’s input. It multiplies the outer function’s derivative, evaluated at the intermediate value $g(x)$, by the inner function’s derivative.

Addition, multiplication, and composition are different ways to combine functions. A larger calculation can use all three. Let’s first follow a change through one path, then see what happens when the same input affects the output through two paths.

Following one path

Suppose we square $x$ to get an intermediate value $a$, then multiply $a$ by $3$ to get $z$:

We want to know how changing $x$ affects the final output $z$. One approach is to substitute $a=x^2$ into $z=3a$, writing the output directly in terms of the input:

The chain rule gives us the same answer using the two relationships we already know. The first function’s derivative tells us how the intermediate value changes; the second function’s derivative tells us how that intermediate change affects the final output:

Multiplying these rates along the path gives:

At $x=3$, the two local slopes are $6$ and $3$, giving a combined slope of $18$. A small change $h$ in $x$ therefore changes $a$ by approximately $6h$, which changes $z$ by approximately $3\times6h=18h$.

Our autodiff engine will record these operations and combine their numerical local derivatives: at $x=3$, it will multiply $3\times6$ to obtain $18$.

Following two paths

Now let’s use $x$ in two places: first add it to $y$ to produce $a$, then multiply $a$ by $x$:

At $x=3,y=4$, the forward computation gives $a=7$ and $z=21$. Changing $x$ now affects $z$ in two ways: it changes the intermediate value $a$, and it changes the $x$ that multiplies $a$. Our derivative must account for both effects while keeping $y$ fixed. For this small expression, we can find the answer by expanding it and differentiating directly:

Using the sum rule to differentiate each term with respect to $x$, we get $2x$ from $x^2$ and $y$ from $xy$. At $x=3,y=4$, this gives $2(3)+4=10$. The product rule gives the same result without expanding: differentiating the two factors in $(x+y)x$ gives $1\cdot x+(x+y)\cdot1=3+7=10$. This tells us that, starting at $x=3$ and keeping $y=4$ fixed, increasing $x$ by a tiny amount $h$ increases $z$ by approximately $10h$.

The chain rule lets us obtain the same derivative by combining the local derivatives along each path from an input to the final output. We multiply the derivatives along each path, as in the previous example, and then add the two contributions:

• Indirect effect, $x\to a\to z$: changing $x$ affects $z$ through the intermediate value $a$. Increasing $x$ by $h$ increases $a=x+y$ by $h$. Holding the multiplication’s other input at $x=3$, that contributes $3h$ to the output change. The local derivatives along this path multiply to $3\times1=3$.
• Direct effect, $x\to z$: the multiplication uses $x$ itself. Holding its other input at $a=7$, increasing $x$ by $h$ contributes $7h$ to the output change. The local derivative along this path is $a=7$.

Changing $x$ affects both paths, so we add their contributions:

We hold one multiplication input fixed only to isolate each contribution. In the complete calculation, both inputs change, giving the combined first-order change $\Delta z\approx3h+7h=10h$. We can see how this relates to the exact change by increasing $x$ from $3$ to $3+h$, while keeping $y=4$:

The original output was $21$, so the exact change is $\Delta z=10h+h^2$. The derivative captures the $10h$ part. The extra $h^2$ comes from both multiplication inputs changing together; it explains why increasing $x$ by a full $1$ changes $z$ by $11$, even though the derivative at the starting point is $10$.

The widget below shows both paths from $x$ to $z$. Use Indirect via a and Direct to z to highlight them. Move Input x to see their contributions change, then increase Change h to $1$ to compare the predicted change of $10$ with the actual change of $11$ at $x=3$. Shrink $h$ to see the approximation improve.

For comparison, a finite difference evaluates the full function at a nearby input:

We can also calculate the sensitivity to $y$ using the same saved values. There is one path, $y\to a\to z$, giving $\partial z/\partial y=3\times1=3$. Neither derivative required another evaluation at a perturbed input.

The pattern is the same for a loss: multiply local derivatives along a path, and add contributions from different paths. Let’s now turn that arithmetic into a numerical trace of our original loss. The engine will record its operations as a graph and work backward through them.

Trace the same loss with autodiff

Return to $L(a)=a^2+3$ at $a=5$, the same function and point used in our comparison table. Symbolic differentiation gave us $L'(a)=2a$, then $L'(5)=10$. Now we’ll calculate that $10$ by following the executed operations:

def loss(a):
    q = a * a
    return q + 3

The forward calculation gives $q=25$ and $L=28$. The addition has local derivative $1$ with respect to $q$. The multiplication has two operand paths, both referring to the same $a$; each local derivative is the other operand’s value, $5$.

Start the backward pass with $dL/dL=1$. Passing through the addition gives $dL/dq=1\times1=1$. Each multiplication path then contributes $1\times5$ to the gradient of $a$, so:

Select Evaluate a × a + 3, then Seed loss gradient. Use Backpropagate addition, followed by Follow left operand and Follow right operand, to watch the gradient reach $10$.

In ordinary Python, the same arithmetic is:

a = 5.0
q = a * a                   # forward: 25.0
loss = q + 3.0              # forward: 28.0

loss_grad = 1.0             # dL/dL
q_grad = loss_grad * 1.0    # derivative of q + 3 with respect to q
a_grad = 0.0
a_grad += q_grad * a        # left operand contributes 5.0
a_grad += q_grad * a        # right operand contributes 5.0
print(a_grad)               # 10.0

Every value needed came from the computation at $a=5$. We never evaluated the loss at $5+h$ and never constructed the expression $2a$ as the output of a symbolic differentiation procedure. We applied the local rules to numbers and accumulated the result.

The automation comes from recording the operations and applying those rules for us. Our Value class will do this when we write loss = a * a + 3 and call loss.backward(). First, let’s connect the same process to a loss built from a model’s predictions and targets.

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 $a^2+3$:

The variable is now the weight $w$, rather than our earlier parameter $a$. 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.

Next, we’ll calculate the same weight derivative by following the model’s operations, as our engine will do.

Which derivative do we take?

For the weight, we follow the path from $w$ through the prediction and error to the loss:

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$.

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)$. We combined the local derivatives by hand; now let’s watch a backward pass do the same calculation on this exact model and loss.

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.

As in our two-path example, contributions add when a value has multiple paths to the output. We’ll use this compact case to check that behavior in the engine:

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.

We choose a central difference because sampling both sides cancels the leading approximation error. For sufficiently smooth functions, forward and backward differences typically have error proportional to $h$, while central difference has error proportional to $h^2$. Halving $h$ therefore roughly halves the one-sided error but reduces central-difference error to about a quarter. For our quadratic loss, central difference is exact in exact arithmetic; for general functions, it remains an approximation.

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. With the engine checked, we can examine how the choice of loss affects those gradients and then use them to update parameters.

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.

Why we recalculate gradients after an update

The update moves the weight, but the gradient was evaluated at its previous value. The chain rule combines derivative rules, and the numbers those rules produce can change after that move. Let’s reuse our loss $L(a)=a^2+3$ and let its input come from a multiplication: $a=wx$. Keep $x=2$ fixed and treat $w$ as the weight we adjust:

The multiplication gives $da/dw=2$, and the loss gives $dL/da=2a$. Combining them:

At $w=1$, we get $a=2$ and $L=7$, so the gradient is $4\times2=8$. If an update moves $w$ to $0.5$, then $a$ becomes $1$ and the gradient becomes $2\times2=4$:

The rules stay the same, but their numerical values can change after an update. The multiplication’s derivative stays $2$ because $x$ is fixed; the loss derivative changes with $a$. Reusing the old gradient of $8$ at $w=0.5$ would describe the slope at the previous weight, not the current one.

Move Weight w below to see these values change. From the initial weight of $1$, select Apply one update to move to $0.5$ and recalculate the gradient. This demonstration uses learning rate $\eta=1/16$, so each update subtracts one sixteenth of the current gradient from the weight.

This connects the two senses in which we refer to derivatives and approximations as local. Each operation supplies a local derivative rule. The resulting gradient describes the loss near the current weight and gives a local approximation of its change. After updating the weights, our engine builds a fresh graph and computes a new gradient. We can now repeat that process to train on several examples.

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