I built a small autograd engine, used it to differentiate one expression, and compared it with PyTorch and a numerical estimate:

forward: L = -0.9051482536
my tape:  dL/da = -0.1807066389
torch64:  dL/da = -0.1807066389
central difference (h=0.1):  dL/da = -0.1815845310

The tape and PyTorch print the same ten decimal places. The central difference with h=0.1 differs by 0.0008778921. Calling the third row wrong would already assume the first two are the oracle. It is more useful to ask what each method computed. The tape applies stored local derivative rules. PyTorch executes its own automatic-differentiation graph. The central difference measures a secant across an interval of width 0.2. Only the limit of that measurement as the interval approaches zero defines the derivative, and binary floating point prevents taking that limit literally.

The discrepancy is therefore the entrance, not a failed check. Following dL/da backward exposes why a reused value needs accumulated contributions, why graph order matters, why one scalar loss makes reverse mode attractive, what a matrix multiply saves for its backward operation, and which tensors a framework must keep alive. The numerical path reaches fp16 only after those mechanics make a tiny gradient necessary.

the first nudge was ordinary subtraction

Start without derivative notation. Let

y = x²

At x=3, the output is 9. Increase x by 0.1 and the output becomes 3.1² = 9.61. The input changed by 0.1 and the output changed by 0.61, so the measured change per unit input is

(9.61 - 9) / (3.1 - 3) = 6.1

Repeating with a smaller nudge moves that ratio toward 6:

nudge       measured output change per input change
 0.1                         6.1
 0.01                        6.01
 0.001                       6.001

For x², expanding (x+h)² shows why:

((x+h)² - x²) / h
= (x² + 2xh + h² - x²) / h
= 2x + h

As h approaches zero, the extra h approaches zero and the ratio approaches 2x. That limiting rate is the derivative. At x=3, it is 6. The notation dy/dx = 2x records which output changes with respect to which input.

The opening expression has two inputs:

z = a*b + a
L = tanh(z)

Changing a while holding b fixed asks for a partial derivative, written dL/da here to keep the typography compact. Changing b while holding a fixed asks for dL/db. The vector containing all partial derivatives of one scalar output with respect to its inputs is its gradient:

gradient of L with respect to [a,b] = [dL/da, dL/db]

The operations can be separated into named intermediate values:

m = a*b
z = m+a
L = tanh(z)

At a=1.5 and b=-2:

m = -3
z = -1.5
L = -0.9051482536

Each line has a local rate. For multiplication m=a*b, changing a by a small amount changes m by approximately b times that amount, so dm/da=b. Changing b gives dm/db=a. For addition z=m+a, either input passes its change through with factor 1. For tanh, the local derivative is 1 - tanh(z)². At z=-1.5, that factor is:

dL/dz = 1 - (-0.9051482536)²
       = 0.1807066389

To reach L from b, a nudge travels through b -> m -> z -> L. Each stage scales the incoming nudge, so the total scale is the product:

dL/db = dL/dz * dz/dm * dm/db
       = 0.1807066389 * 1 * 1.5
       = 0.2710599584

That multiplication of local rates is the chain rule. It is not specific to neural networks. It is what happens when the output of one numeric function becomes the input of another.

The route from a is different because a appears twice. One route is a -> m -> z -> L. The other is the direct a -> z -> L. Their contributions add:

route through m:
  dL/dz * dz/dm * dm/da
  = 0.1807066389 * 1 * -2
  = -0.3614132778

direct route:
  dL/dz * dz/da
  = 0.1807066389 * 1
  = 0.1807066389

sum:
  dL/da = -0.1807066389

This is why the opening derivative is negative even though the last +a route contributes positively. The multiplication route has twice the magnitude and the opposite sign.

The intermediate values and arrows form a directed graph. Each operation creates an output node and records edges to the values it consumed. The forward evaluation computes m, then z, then L. A reverse traversal starts with dL/dL=1 and distributes that quantity toward the leaves. An operation needs its saved forward values and its own local derivative. It does not need to understand the whole expression.

The automatic differentiation survey by Baydin and colleagues distinguishes this numeric graph transformation from symbolic differentiation and finite differences. The program executes ordinary operations while recording enough structure to apply the chain rule. The result is a derivative evaluated at the actual input values, subject to the same floating point arithmetic as the forward program.

the number remembered where it came from

Every scalar is wrapped in a Value. It stores the numeric value, its accumulated gradient, the parent nodes that produced it, and a function that applies one local backward rule. Python’s __add__ and __mul__ method names let the ordinary + and * operators build these nodes:

class Value:
    def __init__(self, x, parents=()):
        self.x = x
        self.grad = 0.0
        self.parents = parents
        self.back = lambda: None

    def __add__(self, other):
        result = Value(self.x + other.x, (self, other))
        def backward():
            self.grad += result.grad
            other.grad += result.grad
        result.back = backward
        return result

    def __mul__(self, other):
        result = Value(self.x * other.x, (self, other))
        def backward():
            self.grad += other.x * result.grad
            other.grad += self.x * result.grad
        result.back = backward
        return result

    def tanh(self):
        value = math.tanh(self.x)
        result = Value(value, (self,))
        def backward():
            self.grad += (1 - value*value) * result.grad
        result.back = backward
        return result

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

        def visit(node):
            if id(node) in seen:
                return
            seen.add(id(node))
            for parent in node.parents:
                visit(parent)
            order.append(node)

        visit(self)
        self.grad = 1.0
        for node in reversed(order):
            node.back()

The nested backward functions are closures. A closure retains access to variables from the operation call after that call returns. The multiplication closure retains self, other, and result, including the two forward values needed by the local rule. A tensor framework saves tensors for the same reason, although it uses native objects and generated operation nodes rather than thousands of Python closures.

visit performs a depth-first traversal. A node enters order only after its parents. That produces a topological order: every forward dependency appears before the node that consumes it. Reversing the list makes every consumer run before its inputs. The seen set matters when a node is reused. Without it, the traversal could insert the same node once per incoming path and apply its backward rule repeatedly.

Setting the final node’s gradient to 1 expresses dL/dL=1. Every earlier gradient begins at zero. A local rule reads the gradient that reached its output, scales it, and accumulates the result in its parents.

The engine runs the opening expression directly:

a = Value(1.5)
b = Value(-2.0)
loss = (a*b + a).tanh()
loss.backward()
my tape:  dL/da = -0.1807066389   dL/db = 0.2710599584

The number agrees with the analytic calculation from the two paths. It is not an “exact analytic gradient” in the mathematical sense. math.tanh and every multiplication are binary64 evaluations. The result is an automatic derivative of the floating point values produced along this path.

PyTorch 2.9.1 produces the same printed values in torch.float64:

torch64:  dL/da = -0.1807066389   dL/db = 0.2710599584

The PyTorch 2.9 backward documentation says the engine differentiates the graph with the chain rule and accumulates results into leaf .grad fields. Agreement here checks one graph and two leaf gradients. It does not validate the tape’s behavior for every operation, mutation, view, device, dtype, or reused graph.

two harmless edits broke the graph

I broke two lines independently and reran the same expression. The first change replaced every += with =. The second traversed the graph in forward topological order rather than reverse order.

The expression L = tanh(a·b + a) is chosen to make this bite: a appears twice, once inside a·b and once in the bare + a. When gradient flows backward, it arrives at a along both of those paths, and the two contributions have to add. That’s what += does. Here’s the correct engine, the overwrite bug (= instead of +=), and the analytic answer, side by side:

exact:            dL/da = -0.1807066389   dL/db = 0.2710599584
correct +=:       dL/da = -0.1807066389   dL/db = 0.2710599584
BUG =(overwrite): dL/da = -0.3614132778   dL/db = 0.2710599584

The program still completes. It leaves dL/db correct because b has one outgoing path. For a, the multiply contribution lands last and overwrites the direct addition contribution. The remaining value is -0.3614132778, exactly the route-through-m result derived above. The factor of two relative to the correct total is a property of b=-2 and the direct +a route in this expression, not a general signature of an overwrite bug.

The second load-bearing line is the ordering. The backward pass has to visit a node only after every node that consumes its output, because a node’s back reads out.grad, and out.grad isn’t finished accumulating until all its children have run. Reverse topological order guarantees that. To see what “guarantees” is buying, I kept the correct += rules and just called the closures in forward build order instead of reversed:

BUG wrong order:  dL/da = 0.0000000000   dL/db = 0.0000000000

Every leaf remains zero because its consumer runs before a gradient reaches that consumer’s output. The root runs last and deposits a gradient into z, but z.back has already run. No second traversal propagates the late value. Neither deliberate bug raises an exception because all operations remain valid floating point assignments. Correctness depends on the graph invariant, not on type checking.

the interval could not become zero

The central difference evaluates the expression at two nearby points:

(L(a+h) - L(a-h)) / (2h)

It is the slope of the chord between those points. The derivative is the limit as h approaches zero. A program must choose a finite, representable h, so the method is always an approximation. I swept powers of ten against the analytic binary64 calculation:

h=1e-01  err=8.78e-04     <- the fourth-decimal disagreement from the top
h=1e-02  err=8.78e-06
h=1e-04  err=8.78e-10
h=1e-05  err=6.68e-12     <- smallest observed error in this sweep
h=1e-06  err=4.88e-11
h=1e-08  err=3.82e-10
h=1e-10  err=7.81e-08
h=1e-12  err=1.78e-05     <- getting worse again

The left side of the valley is truncation error. The function curves across the interval, so the chord differs from the tangent. Reducing h shrinks that error. The right side comes from binary64 evaluation. L(a+h) and L(a-h) become nearly equal. Their leading bits cancel during subtraction, leaving a small difference with fewer trustworthy bits. Division by 2h magnifies the remaining error. The floating point investigation followed the same cancellation mechanism in ordinary arithmetic.

There is no step size that removes both errors. In this declared sweep, 1e-5 has the smallest observed error. A different function, input, or dtype can put the bottom elsewhere. Finite differences remain useful for detecting a sign error, a missing path, or a large factor. They are not a source of exact final bits.

I fitted the slope of log(error) against log(h) over h=1e-1 through 1e-4, before cancellation dominates. Central difference measures slope 2.000. Forward difference, (L(a+h)-L(a))/h, measures 0.993:

forward-diff log-log slope (h=1e-1..1e-4): 0.993  (theory 1)
central-diff log-log slope (h=1e-1..1e-4): 2.000  (theory 2)

A slope near one means reducing h by ten reduces truncation error by about ten. A slope near two means the error falls by about one hundred. The central formula cancels the first-order truncation term by sampling both sides:

forward difference (L(a+h)-L(a))/h :
  h=1e-04  err=1.636e-05
  h=1e-06  err=1.636e-07
  h=1e-08  err=3.819e-10     <- forward's best
  h=1e-10  err=7.810e-08

The forward estimate’s smallest error in this sweep is 3.819e-10 at h=1e-8. The central estimate reaches 6.679e-12 at h=1e-5. Neither minimum is a general recommendation. The table records the error surface of one function in one dtype.

central difference (L(a+h)-L(a-h))/2h :
  h=1e-05  est=-0.180706638930  err=6.679e-12
  h=1e-06  est=-0.180706638875  err=4.883e-11
  h=1e-08  est=-0.180706638542  err=3.819e-10

the subtraction could be removed

For analytic functions, a complex perturbation avoids subtracting two nearby real outputs. Evaluate at a+ih, take the imaginary component, and divide by h:

For a function with a valid analytic continuation, its Taylor series contains:

f(a+ih) = f(a) + ih f'(a) - h² f''(a)/2 - ih³ f'''(a)/6 + ...

The imaginary component divided by h is f'(a) plus an error proportional to h². There is no subtraction of two nearly equal values. I changed math.tanh to cmath.tanh and swept the step:

complex-step Im(L(a+ih))/h :
  h=1e-001  est=-0.1798282148870546  err=8.784e-04
  h=1e-002  est=-0.1806978572767289  err=8.782e-06
  h=1e-006  est=-0.1807066389235607  err=8.785e-14
  h=1e-008  est=-0.1807066389236485  err=2.776e-17
  h=1e-100  est=-0.1807066389236485  err=2.776e-17
  h=1e-200  est=-0.1807066389236485  err=2.776e-17
  h=1e-250  est=-0.1807066389236486  err=0.000e+00

Below h=1e-8, the comparison with the analytic binary64 expression stays within 2.776e-17, and one row happens to subtract to zero. That means both computations produced the same binary64 value in that row. It does not prove either is the exact real-number derivative.

The method has a structural limit. A function must accept complex input and preserve an analytic continuation. I captured three controls:

complex step:
  tanh at 1.5:          0.18070663892364855
  ReLU branch at +2:   1.0        expected local real slope 1
  ReLU branch at -2:   0.0        expected local real slope 0
  ReLU branch at 0:    0.0        real derivative undefined
  abs at +2:           0.0        real derivative is 1

Python cannot compare a complex number with zero, so max(0, 2+ih) raises TypeError. I supplied a ReLU-like branch based on the real component. It reproduces the local slopes away from zero but returns one arbitrary branch value at the kink, where the real derivative does not exist. abs is the silent failure: abs(2+ih) is real, so its imaginary component is zero and the complex-step estimate is zero even though the real derivative at 2 is one.

The method is a strong check for the smooth scalar expression. It is not a drop-in gradient checker for a program containing non-analytic operations, comparisons, casts, indexing decisions, or complex-unsupported kernels.

The dtype comparison remains separate. The tape uses Python binary64 values. PyTorch float32 gives dL/da = -0.1807066202, differing by 1.87e-8. Both computations apply the same derivative expression to rounded intermediate values with different precision.

one checked coordinate left seven places to be wrong

Checking a in the scalar expression says nothing about the tape’s derivative with respect to b. The XOR check perturbs one of seventeen parameters. Selecting a few individual coordinates becomes weaker as a model grows. A single attention matrix in GPT-2 has hundreds of thousands of entries, and a bug can be confined to a slice the checker never touches.

A finite difference can check a whole direction at once. Let the parameters be a vector x, choose another vector v, and evaluate:

L(x + h*v) - L(x - h*v)
-------------------------
           2h

This estimates the rate at which the loss changes when every parameter moves in direction v. If the analytic gradient is g, the chain rule predicts the same directional derivative as a dot product:

g dot v = sum over i of g[i] * v[i]

The dot product works because parameter x[i] moves by v[i] units for one unit of movement along the chosen direction. Every coordinate contribution adds into the one scalar rate.

I built an eight-parameter float64 function:

L(x) = sum(tanh(A*x)²)

The fixed matrix A has six rows. PyTorch computes the eight-entry gradient. Five random direction vectors are normalized to length one, then checked with central differences at h=1e-5:

direction   numerical slope       gradient dot v       absolute error
   0       -0.053887421725        -0.053887421730        5.501e-12
   1        0.019459103751         0.019459103748        2.887e-12
   2        0.790708188125         0.790708188117        7.882e-12
   3        0.748836225295         0.748836225108        1.871e-10
   4       -0.031304150250        -0.031304150274        2.380e-11

The largest error is 1.871e-10. To make the negative control precise, I add 0.1 to gradient coordinate three while leaving coordinate zero untouched. A coordinate-zero check reports zero error. The five directional checks report errors between 0.004327 and 0.06596, so each one exposes the altered coordinate through its nonzero component in that direction.

Random directions do not prove a gradient correct. A bad gradient can be orthogonal to one chosen direction, and a small bad component can fall below finite-difference error. Several independent directions make a large, structured discrepancy harder to miss while requiring two forward evaluations per direction instead of two per parameter.

The checker also needs a deterministic function. If dropout draws a different mask at x+h*v and x-h*v, the numerator contains the effect of both the parameter change and the mask change. Batch statistics, data augmentation, concurrent reductions, and nondeterministic kernels can create the same problem. Holding random state and inputs fixed is part of the measurement, not an optional cleanup.

A point at a nonsmooth boundary requires another interpretation. ReLU at zero has no ordinary derivative, yet an autograd framework must choose a backward value for execution. A finite difference crossing the boundary can disagree with that chosen convention even when the implementation follows its documented rule. Moving the check point away from the boundary tests the local branch. Testing the boundary checks the framework’s convention rather than a unique mathematical derivative.

Gradient checking therefore asks a narrower question than training: does the recorded local backward program agree with independent local measurements at this fixed state? It can catch a missing transpose, an overwrite instead of an accumulation, a sign error, or a wrong scale. It cannot show that the loss surface is easy to optimize or that the gradient remains useful after a large step.

one output changed the useful direction

The tape uses reverse mode, but automatic differentiation does not have to run backward. Forward mode carries a value and its derivative together. A pair

(value, tangent)

records the ordinary value and how that value changes for one selected input direction. Addition adds both components. Multiplication follows the product rule:

(a, da) * (b, db)
= (a*b, da*b + a*db)

To get the derivative with respect to the first input, seed its tangent with 1 and every other input with 0. Every operation propagates that one nudge alongside its value. This is often implemented with dual numbers, but the pair above contains the mechanism.

I applied forward and reverse mode to the same expression with four scalar inputs and one scalar output:

L = tanh(x0*x1 + x2*x2 + x3*x3)

Both modes return:

dL/d[x0,x1,x2,x3]
= [-0.068031297262,
    0.051023472946,
    0.017007824315,
   -0.051023472946]

One forward-mode seed traverses the six graph operations once and gives one directional derivative. Obtaining all four partial derivatives takes four seeds and 24 graph-operation evaluations. Reverse mode runs the six forward operations once, seeds the one output, then invokes six backward closures to produce all four leaf gradients.

forward graph operations per seed:       [6, 6, 6, 6]
forward operations for full gradient:    24
reverse backward closures for gradient:   6
maximum gradient difference:               0

These counts exclude the arithmetic inside each tangent or closure and do not predict wall time in a tensor framework. They expose the scaling choice. Forward mode is attractive when there are few input directions and many outputs. Reverse mode is attractive when many parameters feed a small number of outputs. Training usually reduces a batch to one scalar loss, while the model can have millions or billions of parameters. One reverse seed can reach every parameter.

For a vector output y and vector input x, all pairwise partial derivatives form the Jacobian matrix J, where J[i,j] is dy[i]/dx[j]. Materializing that whole matrix can be enormous. Reverse mode instead computes a vector-Jacobian product: an upstream vector v weights the output rows and the engine returns v^T J. A scalar loss is the special case with one output and seed 1.

This is why PyTorch requires an explicit gradient argument when backward starts from a non-scalar tensor. The argument is the upstream vector, not a request to build the full Jacobian. The same object appeared as out.grad in the scalar tape. Every local backward rule consumes an upstream gradient shaped like its output and returns contributions shaped like its inputs.

the transpose appeared because indices had to match

The model in the previous investigation repeatedly multiplies matrices. Let

Y = X W

with X shaped (2,3) and W shaped (3,2):

X = [[1,2,3],       W = [[  1, -1],
     [4,5,6]]            [  2,  0],
                          [0.5,  3]]

Y = [[ 6.5,  8],
     [17.0, 14]]

Backward does not begin with dY=1 unless Y is itself the scalar loss. Suppose later operations send this upstream gradient:

G = dL/dY = [[ 1, 2],
             [-1, 0.5]]

Take one input value, X[0,1], which is 2. It contributes to both outputs in row 0:

Y[0,0] includes X[0,1] * W[1,0]
Y[0,1] includes X[0,1] * W[1,1]

The two paths contribute:

dL/dX[0,1]
= G[0,0] * W[1,0] + G[0,1] * W[1,1]
= 1*2 + 2*0
= 2

Repeating that path sum for every input produces:

dL/dX = G W^T
      = [[-1.0,  2.0, 6.5],
         [-1.5, -2.0, 1.0]]

The transpose makes the inner dimensions match and aligns each output coordinate with the weight that connected it to an input coordinate.

One weight, W[0,1], is reused by both input rows:

Y[0,1] includes X[0,0] * W[0,1]
Y[1,1] includes X[1,0] * W[0,1]

Its contributions add:

dL/dW[0,1]
= X[0,0] * G[0,1] + X[1,0] * G[1,1]
= 1*2 + 4*0.5
= 4

All weight gradients can be written:

dL/dW = X^T G
      = [[-3.0, 4.0],
         [-3.0, 6.5],
         [-3.0, 9.0]]

The probe computes both formulas from the arrays and compares them with PyTorch float64 gradients. Every entry is equal, with maximum difference 0. The two transposes are not backward conventions to memorize. They fall out of which row and column participated in each forward output.

The operation also reveals what backward needs from forward. Computing dL/dX needs W. Computing dL/dW needs X. If those values are not otherwise available, the autograd node must retain them. For a transformer MLP, saved activation rows can be far larger than the small example. That retained state becomes the memory problem later in the article.

the target logit moved every other logit

The GPT-2 loss used below starts from one target token, but its first backward value is not confined to that token. The model produces one logit for each possible next token. A logit is an unrestricted score. Softmax turns the scores into positive probabilities that sum to one:

              exp(z[i])
p[i] = -----------------------
        sum over k of exp(z[k])

For three logits [2.0, 0.5, -1.0], the direct calculation gives:

p = [0.785597034589,
     0.175290392140,
     0.039112573271]

The first class is the target. Its negative log likelihood is:

L = -log(p[0])
  = 0.241311296657

The logarithm changes a product of probabilities across several tokens into a sum of token losses. It also assigns a large penalty when the probability of the observed target is close to zero. The derivative can be found without differentiating a quotient for every pair of logits. Substitute the softmax definition into the loss:

L = -log(exp(z[target]) / sum(exp(z[k])))

  = -z[target] + log(sum(exp(z[k])))

The derivative of the first term with respect to logit z[j] is negative one when j is the target and zero otherwise. For the second term, let S = sum(exp(z[k])). Its derivative is:

d log(S) / dz[j]
= (1 / S) * dS/dz[j]
= exp(z[j]) / S
= p[j]

Combining the two paths gives one compact rule:

dL/dz[j] = p[j] - 1        when j is the target
          = p[j]            otherwise

For the three logits:

dL/dz = [-0.214402965411,
          0.175290392140,
          0.039112573271]

PyTorch float64 agrees with each entry to within 6.939e-18 in the probe. The gradient entries sum to 2.776e-17, which is zero at the rounding scale of this calculation. That near-zero sum is not an accident. Adding the same constant to every logit leaves all softmax probabilities unchanged because the new exponential factor cancels between numerator and denominator. A movement that raises every score equally cannot change the loss, so the gradient has no component in that direction.

The signs make the update concrete. Gradient descent subtracts the gradient. It raises the target logit because its gradient is negative. It lowers every other logit because their gradients are positive. The size of each non-target push equals that token’s current probability. A wrong token receiving 0.17 probability is pushed down more strongly than one receiving 0.039.

Softmax also creates a numerical trap before backward begins. Evaluating exp(1000) in binary64 overflows, even though the probabilities for logits [1000, 999, 998] are ordinary:

stable probabilities:
[0.665240955775, 0.244728471055, 0.090030573170]

Subtracting the largest logit before exponentiation changes the vector to [0, -1, -2]. Softmax is unchanged by that common shift, and every exponential is now at most one. The probe’s naive Python calculation raises OverflowError; the shifted calculation succeeds. A fused cross entropy implementation can compute the equivalent logsumexp expression without forming an unstable probability followed by a separate logarithm.

The upstream vector then meets the output matrix. If the last hidden row is h with width D, and each vocabulary row W[j] produces one logit,

z[j] = sum over d of h[d] * W[j,d]

then the output weight gradient is an outer product:

dL/dW[j,d] = dL/dz[j] * h[d]

Every vocabulary row receives a gradient, not only the target row. The hidden state receives a weighted sum of all output rows:

dL/dh[d] = sum over j of dL/dz[j] * W[j,d]

That vector is the first upstream gradient entering the final normalization and the transformer blocks. It has exactly the shape of the final hidden row. The vocabulary dimension disappears because contributions from all logits add into it.

For a batch or a sequence, the program must also declare a reduction. Summing token losses makes the gradient scale with token count. Taking their mean divides every token contribution by that count. Ignored padding tokens change the denominator if the implementation defines the mean over non-ignored tokens. This is why two runs can use the same examples and still have different gradient magnitudes if one averages per token and another averages per sequence.

The small loss probe uses one row, one target, and the default mean, which is identical to the single token loss. The GPT-2 probe below also uses one target position. That keeps the first upstream vector traceable. It does not represent the weighting used by a complete language model training batch.

the softmax sent the error down three roads

The output gradient reaches an attention block as one tensor, then splits through values, attention weights, scores, queries, and keys. Writing those paths for one head exposes more accumulation than the scalar expression did.

The forward operations are:

S = Q K^T / sqrt(D)
A = softmax each row of S
O = A V

Q, K, and V each have one row per token and D columns for one head. Row i of S contains the scores from query position i to every key position. Row i of A contains probabilities that sum to one. Row i of O is a weighted sum of all value rows.

I used two positions and width two:

Q = [[1.0, 0.0],    K = [[1.0, 1.0],    V = [[2.0, -1.0],
     [0.5, 1.0]]         [0.0, 2.0]]         [0.0,  3.0]]

The score scale is sqrt(2). Forward produces:

S = [[0.7071, 0.0000],
     [1.0607, 1.4142]]

A = [[0.6698, 0.3302],
     [0.4125, 0.5875]]

O = [[1.3395, 0.3210],
     [0.8250, 1.3499]]

Suppose the rest of the model returns this upstream gradient:

G = dL/dO = [[ 1.0, -0.5],
             [-1.0,  2.0]]

The last forward operation was O=A V, so the matrix rule already derived gives:

dL/dV = A^T G
dL/dA = G V^T

The numbers are:

dL/dV = [[ 0.2572, 0.4902],
         [-0.2572, 1.0098]]

dL/dA = [[ 2.5, -1.5],
         [-4.0,  6.0]]

One value row is reused by every query row. Its gradient adds contributions from both rows of G, weighted by how much each query attended to that value. One attention weight affects every output feature in its row. Its gradient is the dot product between that row of G and the corresponding value row.

Softmax is the part where changing one input changes every output in the same row. For one row with probabilities a and upstream gradient g, the pairwise derivative is:

da[j]/ds[k] = a[j] * (1 - a[j])   when j equals k
            = -a[j] * a[k]         otherwise

A score raises its own probability and lowers the others because the row must still sum to one. Multiplying that Jacobian by the upstream vector can be written without constructing a square matrix:

dot = sum over j of g[j] * a[j]
dL/ds[k] = a[k] * (g[k] - dot)

The shared dot subtracts the component that would move all logits together. Applying it to both rows gives:

dL/dS = [[ 0.8847, -0.8847],
         [-2.4235,  2.4235]]

Each row sums to zero within floating point error. The measured row sums are 1.110e-16 and -4.441e-16. This is the same shift invariance seen in the vocabulary softmax. Adding one constant to every score in a row cannot change its probabilities, so backward removes that direction.

The score matrix came from Q K^T / sqrt(D). Applying the matrix rule and the scale gives:

dL/dQ = (dL/dS) K / sqrt(D)
dL/dK = (dL/dS)^T Q / sqrt(D)

For the small arrays:

dL/dQ = [[ 0.6256, -0.6256],
         [-1.7137,  1.7137]]

dL/dK = [[-0.2312, -1.7137],
         [ 0.2312,  1.7137]]

The hand formulas match PyTorch float64 with maximum differences 2.220e-16 for queries, 2.776e-17 for keys, and zero for values. Those differences are rounding in operation order. They are not a tolerance chosen after a failed float32 check.

The division by sqrt(D) appears in backward too. It scales both query and key gradients because it scaled every score in forward. Omitting it would make those two gradients too large by sqrt(D) while leaving the value gradient correct. A test that inspected only V would miss the bug.

GPT-2 adds a causal mask before softmax. A future position receives a value that acts like negative infinity, so its softmax probability is zero. In exact arithmetic, its local contribution remains zero. Implementations use a large negative finite value or a specialized masked kernel, and the backward definition must follow that representation without producing an invalid infinity subtraction. The two-position probe deliberately has no mask, so it checks the dense formulas rather than claiming to exercise that kernel path.

Queries, keys, and values are not independent leaves in the model. They are three projections of the same residual input X:

Q = X Wq
K = X Wk
V = X Wv

Backward produces three contributions to the residual:

dL/dX = dL/dQ Wq^T
       + dL/dK Wk^T
       + dL/dV Wv^T

That plus sign is the original reused a problem at matrix scale. Overwriting one route would erase information from an entire projection. The three weight gradients similarly reuse X:

dL/dWq = X^T dL/dQ
dL/dWk = X^T dL/dK
dL/dWv = X^T dL/dV

Multiple heads repeat these operations on slices, concatenate their outputs, and pass them through another projection. Backward splits the concatenated gradient into head slices, runs each head’s paths, then adds their input contributions. A fused kernel may never materialize these named matrices in memory. It must still produce a result equivalent to these dependencies, subject to its declared dtype and numerical tolerance.

Layer normalization couples features in another direction. For one token row with D features, it computes:

mean = sum(x[i]) / D
variance = sum((x[i] - mean)²) / D
xhat[i] = (x[i] - mean) / sqrt(variance + epsilon)
y[i] = gamma[i] * xhat[i] + beta[i]

Changing one feature changes the mean and variance used by every other feature. Backward cannot treat this as D independent scalar divisions. Let g be dL/dy, let u = g*gamma, and let r = 1/sqrt(variance+epsilon). Combining the direct normalization path with the mean and variance paths gives:

dL/dx[i] = r/D * (
    D*u[i]
    - sum(u)
    - xhat[i] * sum(u*xhat)
)

The two sums are reductions across the whole feature row. They remove the components caused only by a common shift and a common rescaling of the row. The parameter gradients are simpler:

dL/dgamma[i] = g[i] * xhat[i]
dL/dbeta[i]  = g[i]

The float64 probe uses x=[1,2,4], gamma=[0.5,1.5,-1], beta=[0.25,-0.5,1], and upstream gradient [1,-2,0.5]. Its forward statistics are:

mean:                        2.333333333333
variance:                    1.555555555556
inverse standard deviation: 0.801781148588
xhat: [-1.069041531450, -0.267260382863, 1.336301914313]

The hand formula and PyTorch return:

dL/dx:
[ 1.088132295123, -1.632197154115, 0.544064858993]

dL/dgamma:
[-1.069041531450, 0.534520765725, 0.668150957156]

dL/dbeta:
[1.0, -2.0, 0.5]

Maximum differences are 2.220e-16, zero, and zero. The three input gradients sum to -2.220e-16. Adding a constant to every input leaves the normalized row unchanged, so the backward result has no meaningful component in that common-shift direction.

epsilon prevents division by zero when every feature is equal, but it also changes the derivative. If the variance is much larger than epsilon, its effect is small. If the variance is comparable to epsilon, the inverse standard deviation and every input gradient depend materially on the chosen constant. The probe uses 1e-5, matching the pinned GPT-2 configuration. Changing it would define a different forward function, so a gradient check must use the same value on both paths.

The formula also explains why accumulating its reductions in a wider dtype can matter. Mean, variance, sum(u), and sum(u*xhat) combine many features. Rounding in those reductions is then broadcast back into every element of the row. The float64 probe tests the algebra with little rounding pressure. It does not establish the error of a fused float16, bfloat16, or float32 kernel.

GPT-2 applies normalization before each attention or MLP sublayer and adds the sublayer output back to the residual. If one block is written:

y = x + F(layer_norm(x))

then the upstream gradient reaches x by two routes. One is the direct residual route. The other runs through F and the coupled normalization formula. The direct route is why the derivative contains an identity term, but the second route can reinforce or oppose it. The measured residual norms below are not expected to fall monotonically with depth.

Layer normalization backward needs the normalization result or enough saved statistics to reconstruct it. A fused implementation may save the mean and inverse standard deviation rather than every named intermediate. This is another place where the mathematical graph specifies dependencies while the implementation chooses which values to retain and which to recompute.

the tape reached the model from the previous run

The scalar rules and matrix transposes should reach the exact GPT-2 artifact from one forward pass. I loaded the same revision and made one artificial training target: increase the next-token probability of Paris after The capital of France is. The starting probability was 0.0322432, so its negative log likelihood was 3.434449.

The run used gpt2 revision 607a30d783dfa663caf39e06633721c8d4cfcd7e, Transformers 4.57.1, PyTorch 2.9.1, eager attention, float32 model weights, evaluation mode, four Torch CPU threads, and the Apple M4 in the current machine. The model reports 124,439,808 unique parameter elements. Evaluation mode disables its dropout; use_cache=False prevents inference cache state from changing this backward path. No GPU participated in the measurement.

This is not a claim that the prompt came from GPT-2’s training data with Paris as its observed target. It is a controlled loss chosen to produce one gradient through the pinned forward path.

I retained the thirteen residual states, called backward once, and inspected the last-token gradient norm at each state:

state  0: 2.18991232    state  5: 1.05256379    state 10: 0.29460910
state  1: 0.41888338    state  6: 0.95679313    state 11: 0.25041935
state  2: 0.74908227    state  7: 0.74968779    state 12: 2.86155653
state  3: 0.98113734    state  8: 0.50454092
state  4: 1.03733277    state  9: 0.36300051

The norms are not monotonic. A residual state participates in normalization, attention, MLP, residual, and final projection paths. Norm alone does not tell which parameter update is important or whether a layer is vanishing. It establishes that an upstream gradient reached every retained state.

The combined attention projection stores query, key, and value weights beside one another. Splitting its gradient on the output dimension gives:

layer       ||dL/dWq||    ||dL/dWk||    ||dL/dWv||
  0          1.309296      1.447122      15.353137
  5          9.655268      6.891427      14.710448
 11          3.499980      3.177557       5.239401

These are Frobenius norms over whole slices, not per-weight averages. The value slice has the largest norm in the three sampled layers except neither conclusion generalizes to other prompts, losses, or layers without a sweep.

I checked the direction on layer zero’s complete query slice. Moving the slice by L2 distance 0.01 against its normalized gradient lowered the loss from 3.434449 to 3.421283. Moving the same distance with the gradient raised it to 3.447407. Copying the original weights back restored 3.434449 exactly in the rerun.

That local intervention gives the gradient its operational meaning. For a sufficiently small step, negative gradient is a descent direction for this loss. It does not say that taking repeated fixed steps of 0.01 would train the model, preserve other prompts, or remain in the local linear region.

The scalar tape still creates one Python object, parent tuple, and closure for every elementary operation. I compared a length-1,000 dot product followed by tanh, forward and backward, with one Torch thread. After warmup, I randomized eleven rounds and used the median:

length-1000 dot + tanh, forward+backward
three process runs, eleven randomized rounds within each:

  Python tape median range:  6474.8 to 6557.5 us
  Torch median range:           8.7 to   10.2 us
  per-process ratio range:    642.6 to  741.8 times

The range belongs to this microbenchmark, Python version, Torch build, and M4 CPU. It is not a framework throughput claim. The spread between process runs is also why one exact slowdown would imply more stability than the measurement has. The tape expresses the dot as thousands of scalar object operations. Torch expresses it as tensor operations over contiguous storage and sends the work to compiled kernels. Both use reverse accumulation, but that shared algorithm does not make their execution paths comparable at the instruction level.

the old gradient was still in the leaf

The tape initializes grad to zero, but it never clears gradients before a second backward call. PyTorch deliberately behaves the same way for leaf tensors. I ran backward twice on one retained graph:

maximum |gradient after two calls - 2*gradient after one call| = 0

This is accumulation, not stale memory by accident. A parameter can receive contributions from several graph paths in one backward call and from several backward calls before an optimizer step. The same += rule supports both.

It also creates a common failure. If an optimizer step is supposed to use one batch and the program forgets to clear .grad, the next step includes the previous batch’s gradient. PyTorch’s documentation explicitly says backward accumulates into leaf fields. Calling

optimizer.zero_grad(set_to_none=True)

sets those fields to None. A later backward can allocate them on first contribution. Setting them to numeric zero also works, but can require touching every gradient buffer. In the probe, every parameter gradient is None after the call.

Intentional accumulation is useful when one full batch does not fit in memory. Split seven examples into microbatches of three and four. The full loss is the mean squared error over all seven. If each microbatch computes its own mean, its backward contribution must be weighted by its fraction of the full batch:

for start, stop in [(0, 3), (3, 7)]:
    micro_loss = mse(
        model(inputs[start:stop]),
        targets[start:stop],
        reduction="mean",
    )
    fraction = (stop - start) / 7
    (fraction * micro_loss).backward()

The fractions are 3/7 and 4/7. Their weighted sum is the seven-example mean. On a float64 linear layer:

full-batch gradient norm:                 1.2828833995
weighted microbatch maximum difference:  1.11e-16

The remaining difference is at binary64 rounding scale. The negative control calls backward on each microbatch mean without weighting:

unweighted microbatch maximum difference: 0.6625104282
unweighted/full gradient-norm ratio:       2.0097860545

Adding two means gives each microbatch equal weight even though one has three examples and the other has four. It also roughly doubles the scale because two mean losses are summed. Dividing by the number of microbatches would fix equal-size microbatches, but it is wrong for the unequal split used here. Weighting by example count states the invariant directly.

Gradient accumulation can reduce the number of examples whose activations are live at once. It does not shrink model parameters, their gradient buffers, or optimizer state. It can also differ from a true large-batch step when the forward path depends on batch composition, random masks, running statistics, or a loss that couples examples. The exact equality above belongs to a deterministic linear layer with a separable mean loss.

The same scale problem appears when examples are split across processes instead of across backward calls. I simulated two data-parallel ranks holding three and four examples. Each rank has an identical two-output linear model and computes its local gradient.

If both ranks take a local mean and the system averages the two rank gradients equally, each rank receives half the influence:

naive global gradient = (rank0_mean + rank1_mean) / 2

That gives each of the three examples on rank zero weight 1/6. Each of the four examples on rank one receives weight 1/8. A global seven-example mean should give every example weight 1/7. The measured maximum gradient error is 0.16242961.

Weighting by local example count restores the intended mean:

global gradient = (3*rank0_mean + 4*rank1_mean) / 7

The maximum difference from a single-process seven-example gradient falls to 2.220e-16. An equivalent convention has each rank compute a sum, adds those sums across ranks, then divides by the global number of loss elements. Because the probe’s mean squared error has two output elements per example, that denominator is 7*2. It also agrees within 2.220e-16.

This arithmetic precedes any choice of communication algorithm. An all-reduce can return a sum or a framework can divide by the world size around it. Either is correct only when the loss reduction and local batch sizes make the final scale equal to the intended objective. Equal local batch sizes make an average of rank means correct, which can hide the assumption until one rank receives a short final batch, a filtered example, or a variable number of tokens.

Language model batches make the denominator especially easy to lose. One rank may have the same number of sequences but fewer non-padding target tokens. Averaging rank-level mean token losses equally weights ranks, not tokens. A global token mean requires the summed loss numerator and valid-token count to be reduced consistently.

The simulator runs in one process and performs ordinary tensor additions. It does not measure all-reduce latency, overlap, bucket scheduling, or floating point differences caused by another reduction tree. It isolates the invariant that communication must preserve: every local contribution must arrive once, and the final divisor must describe the objective actually being optimized.

the leaf was a different kind of node

The small tape stores a grad field on every Value, so inspecting an intermediate after backward works by default. PyTorch separates the gradient needed while executing the graph from the gradient retained for the user. A tensor produced by an operation has a grad_fn that points into the backward graph. A parameter created directly by the user is a leaf. It has no creating operation, but it can own a persistent .grad buffer.

This program makes the distinction visible:

x = torch.tensor([2.0, 3.0], requires_grad=True)
y = x * x

print(x.is_leaf, x.grad_fn)  # True, None
print(y.is_leaf, y.grad_fn)  # False, MulBackward...

y.sum().backward()
print(x.grad)                # tensor([4., 6.])
print(y.grad)                # not retained by default

The backward operation needs dL/dy in order to compute dL/dx, but that does not require keeping y.grad after its consumer has run. Releasing intermediate gradient buffers shortens their lifetime. Calling y.retain_grad() asks PyTorch to retain this one:

leaf grad:             [4.0, 6.0]
retained nonleaf grad: [1.0, 1.0]

The nonleaf gradient is one because the scalar loss was y.sum(). Each element of y contributes to the sum with local derivative one.

The PyTorch 2.9 autograd mechanics note states that the forward pass builds a graph of Function objects and that only leaf tensors with requires_grad=True accumulate into .grad by default. It also says the graph is recreated on each iteration. An if statement can select a different operation path on the next example because eager autograd records the path that actually executed. The backward pass cannot send a gradient through a branch that the forward pass did not take.

The versioned C++ source shows the machinery that replaces the tape’s one Python list. In engine.cpp for PyTorch 2.9.1, a backward invocation is represented by a GraphTask. Work is placed in ready queues associated with devices. The engine counts unresolved dependencies and schedules a node when the contributions expected from its consumers are ready. Input buffers collect those contributions before the node function runs.

That dependency count performs the same logical job as the tape’s reversed topological list. Consider a reused value:

             square
            /      \
x -> branch A      branch B -> add -> loss

The square node cannot run backward as soon as branch A contributes. Branch B still owes another contribution. Running early would send an incomplete gradient to x. The small tape avoids that by giving every consumer a chance to run before its parent. The engine can execute independent ready nodes without imposing one global serial list, but it must preserve the same dependency invariant.

At the boundary, an AccumulateGrad node owns the policy for adding an incoming gradient to a leaf’s existing buffer. The implementation has cases for stealing a new gradient buffer, cloning it when layout or reference-count conditions require that, and adding into an existing buffer. Sparse, dense, and higher-order differentiation conditions make this more complicated than self.grad += contribution. The observable contract in this article is narrower: repeated backward calls accumulate, and clearing the field changes the next call from addition to first assignment.

Device execution adds another constraint absent from the Python tape. engine.cpp documents ready queues per worker device and stream handling for CUDA and other device backends. A gradient tensor produced on one stream cannot be consumed on another until the required ordering is established. Graph order and device completion order are related but not identical. A node being next in the mathematical dependency graph does not mean its device work has already completed.

I did not instrument these C++ queues or run a CUDA backward trace here. The queue and accumulator statements come from the pinned source and official mechanics note. The local probes establish the user-visible consequences: path contributions add, leaves retain gradients, nonleaves do not retain them unless asked, and a second call adds to the first call’s result.

the backward pass asked for values that had been discarded

The matrix rule needs forward X to compute dW and forward W to compute dX. GELU backward needs enough information to evaluate its local slope. Attention backward needs queries, keys, values, and softmax-related state or a way to recompute them. Keeping every useful intermediate makes backward faster but can make activation memory larger than the parameter payload.

Gradient checkpointing chooses selected boundaries to keep. Intermediate values inside a checkpointed region are discarded after forward. During backward, the framework reruns part of that region to recreate the values needed by local derivative rules. The PyTorch checkpoint documentation describes this as trading compute for memory and requires the use_reentrant choice to be explicit in version 2.9. The probe uses use_reentrant=False.

I built twelve residual MLP blocks with width 256 and batch 32, then split them into six checkpoint segments. A saved-tensor hook counted every tensor payload the autograd graph requested to retain. It does not measure allocator reservation or process peak RSS:

                                         plain       checkpointed
saved tensors                               61              21
saved-tensor payload                    13.91 MiB         2.50 MiB
block forward calls                        12              22
median time range over 3 runs        2.56 to 2.97 ms  4.52 to 4.73 ms

The checkpointed graph requests 17.98% as many saved bytes. Block forward calls rise from 12 to 22 because backward recomputes blocks inside segments. Non-reentrant checkpointing can stop recomputation once all needed tensors have been recreated, so the count does not have to reach 24.

The payload counter records every tensor presented to the saved-tensor hook. That is closer to “bytes the backward definitions asked to retain” than to “bytes newly allocated.” Two saved tensor objects may refer to the same storage. An allocator may reserve more memory than the live tensors request. Temporary workspace outside the saved set can raise the actual peak. The hook is still useful for comparing these two executions because the model, inputs, allocator, and counting rule remain fixed.

The shapes explain why activations eventually dominate. A residual tensor with batch B, sequence length T, hidden width D, and element size s occupies:

B * T * D * s bytes

For the pinned GPT-2 configuration, D=768. One fp32 residual state for the five-token prompt occupies:

1 * 5 * 768 * 4 = 15,360 bytes

Thirteen states, one before the first block and one after each of twelve blocks, total 199,680 bytes. This is a derived lower bound for those particular states, not the model’s total saved activation memory. At sequence length 1,024, one such state is 3 MiB and thirteen are 39 MiB. Increasing the batch to eight makes that 312 MiB. Parameter bytes did not change; activation bytes grew with both batch and sequence length.

Attention introduces a different shape. Materializing one score per head for each query and key position uses B*H*T*T elements. With twelve heads, sequence length 1,024, batch one, and fp32 elements:

1 * 12 * 1024 * 1024 * 4 = 50,331,648 bytes = 48 MiB

That is per layer. Twelve dense layers would make the score payload 576 MiB if each layer materialized and retained one full tensor of that shape. The five-token probe’s corresponding tensor is only 1,200 bytes per layer. Sequence length grew by about 205 times, while the score element count grew by about 41,943 times because both position axes grew.

This calculation is not a measurement of the pinned implementation. A fused attention kernel can avoid materializing the complete score matrix, save different statistics, or recompute parts during backward. Masks, dtype, attention layout, and kernel choice change the actual storage. The formula identifies the pressure those implementations try to remove.

The checkpoint probe contains residual MLP blocks, not attention. Its 82 percent hook-payload reduction cannot be applied to a transformer by percentage. Checkpoint boundaries also matter. More segments retain more boundary tensors but shorten recomputation intervals. Fewer segments save fewer boundaries and redo more work. The useful configuration is constrained by the peak live set, not by the total sum of tensor sizes across the whole iteration.

Each process run contains eleven randomized timing rounds. The per-process checkpoint-to-plain median ratio ranges from 1.53 to 1.79. That interval is specific to this CPU model and short benchmark. The output arrays match exactly and the maximum parameter-gradient difference is zero in every verification run. Those controls matter because recomputation can be incorrect if the function depends on changed global state or produces different random behavior. PyTorch preserves random-number-generator state by default for checkpointed regions, which itself has overhead.

The original checkpointing paper, Training Deep Nets with Sublinear Memory Cost, analyzes schedules that retain selected boundaries and recompute between them. My hook result is not a reproduction of its large-network memory or runtime numbers. It is the mechanism on a local graph: fewer saved tensor bytes, more forward calls, identical checked gradients.

the saved value changed while the graph waited

Retaining a tensor is useful only if backward sees the same value that forward used. I made y = sum(x*x), changed x in place, and then called backward:

x = torch.tensor([2.0, 3.0], requires_grad=True)
y = (x * x).sum()

with torch.no_grad():
    x.add_(1.0)

y.backward()

The multiplication’s derivative needs the original x, [2,3]. The storage now contains [3,4]. Silently using the new values would return [6,8] instead of the forward expression’s gradient [4,6]. PyTorch 2.9.1 rejects the call:

one of the variables needed for gradient computation has been modified
by an inplace operation: [torch.DoubleTensor [2]] is at version 1;
expected version 0 instead.

The abbreviated message in the probe continues with a suggestion to enable anomaly detection. The autograd mechanics note explains the check. A tensor has a version counter. Saving a tensor for backward also saves its current version. An in-place operation increments the counter. Retrieving a saved tensor compares the versions and raises when they differ.

torch.no_grad() in the example prevents the update itself from becoming a new differentiable operation. It does not prevent the storage mutation or the version increment. Grad mode answers whether operations are recorded. It does not turn an in-place write into a harmless write.

The reason this check cannot be replaced by a blanket ban is that some backward formulas do not need every forward input. Addition can send the same upstream gradient to both inputs without reading their old numeric values. Multiplication needs the other operand. Exponential backward may save its output rather than its input because the derivative of exp(x) equals the output. The framework’s generated backward definition determines what is saved. Code that mutates an unrelated or unsaved tensor may run, while a seemingly similar mutation of a saved tensor fails.

Views make the storage question less obvious. A transpose or slice can have a different shape and stride while referring to the same underlying bytes. Writing through one view can invalidate a value saved through another. Object identity alone is therefore insufficient; the version follows the shared storage relationship that matters to autograd.

Saved values have another lifetime boundary. After a normal backward call, PyTorch releases graph state that is no longer needed. Calling backward again on the same y produces:

Trying to backward through the graph a second time
(or directly access saved tensors after they have already been freed).

This is distinct from gradient accumulation. The leaf .grad field persists, but the operation graph that could produce another contribution has been released. Passing retain_graph=True to the first backward retains that graph. In the local control, clearing x.grad between calls and running the retained graph a second time reproduces [4,6].

Retaining a graph extends the lifetime of saved tensors and therefore its memory cost. Rebuilding the forward pass usually creates a fresh graph with fresh values, which is why an ordinary training loop does not retain every iteration. A program should retain the graph only when another traversal of that exact execution is necessary.

Higher-order differentiation asks for something stronger. The first backward calculation must itself be recorded as differentiable operations. For f(x)=x³ at x=2, the probe uses torch.autograd.grad with create_graph=True:

first derivative:  3x² = 12
second derivative: 6x  = 12

retain_graph keeps the old derivative graph available. create_graph records a new graph for the derivative computation so another derivative can be taken. They solve different lifetime problems. The thirty-line tape stores local derivative arithmetic as ordinary Python floats inside closures, so its backward pass cannot be differentiated again. Supporting higher-order derivatives would require its backward rules to build Value operations instead of consuming the graph into numbers.

seventeen parameters found a bad corner

I attached the scalar engine to a two-input, four-hidden-unit, one-output MLP with tanh activations. Its seventeen parameters try to fit the four XOR rows. A linear decision boundary cannot separate XOR, so the hidden nonlinear layer is necessary.

I added subtraction and negation to Value, formed a sum of squared errors, and checked one weight before training:

gradcheck W1[0][0]: tape=0.00007575  central-diff=0.00007575  |diff|=9.89e-13

The central difference differs from the tape by 9.89e-13 at h=1e-5. That checks one derivative at one parameter state. It does not prove all seventeen gradients or every training step, but it makes a gross backward-rule error less likely.

Full-batch gradient descent with learning rate 0.5 stalls for seed 0:

step    0  loss 5.260514
step  400  loss 4.000015
step 2000  loss 3.999534
final predictions (want -1,+1,+1,-1):
  0,0 -> -0.9883  (target -1)
  0,1 -> +0.9985  (target +1)
  1,0 -> +0.9999  (target +1)
  1,1 -> +0.9998  (target -1)

Three rows are close to their targets. The (1,1) error is almost 2, whose square accounts for almost all of the loss near 4. The largest absolute parameter gradient after 2,000 steps is 1.757e-3. That is an observed small gradient, not by itself proof that every hidden unit saturated or that no descent path exists.

I repeated the fixed architecture, uniform initialization, learning rate, and 2,000-step budget for ten seeds:

plain SGD lr=0.5, 2000 steps, per seed final loss:
  seed 0: 3.9995  STUCK    seed 1: 0.0001         seed 2: 0.0002
  seed 3: 0.0002          seed 4: 3.9998  STUCK   seed 5: 4.0000  STUCK
  seed 6: 8.0000  STUCK    seed 7: 4.0000  STUCK   seed 8: 3.8899  STUCK
  seed 9: 4.0000  STUCK
stuck 7/10 seeds

Seven finish above the probe’s arbitrary failure threshold of 0.5. This is not a success-rate estimate for XOR networks. Changing width, activation, initialization scale, learning rate, loss, or step count changes the experiment. It establishes sensitivity within one declared protocol.

Changing seed 0 to Adam with global step size 0.05 reaches:

Adam lr=0.05 on seed 0: final loss 0.000089
  0,0 -> -0.9973    0,1 -> +0.9952    1,0 -> +0.9948    1,1 -> -0.9944

The optimizer changes the first update, which changes the next parameter state and therefore every later gradient. This run shows a rescue under one Adam setting. It does not show that Adam is the only optimizer that can solve XOR or that full-batch gradient descent cannot succeed under another schedule.

the first optimizer comparison was rigged by one setting

The update

weights = weights - learning_rate * gradient

is full-batch gradient descent when gradient is computed from every example. It becomes stochastic gradient descent when the gradient comes from a sampled example or minibatch. The existing quadratic probe uses all 200 rows every step, so calling it SGD was incorrect.

Its second input feature is scaled by 100. Mean squared error produces a two-dimensional quadratic whose Hessian, the matrix of second derivatives, has eigenvalues:

lambda_min =     1.840138
lambda_max = 21180.873551
condition number = 11510.48

An eigenvector is a direction in which the Hessian only scales rather than turns a vector. Along a direction with curvature lambda, one gradient-descent error component is multiplied each step by:

1 - learning_rate * lambda

For that component to shrink rather than oscillate outward, its absolute value must remain below one. With positive curvature this requires:

0 < learning_rate < 2/lambda

The most curved direction sets the global ceiling:

2 / 21180.873551 = 9.44248e-5

A larger fixed step is unstable on this exact quadratic. A safe step remains slow along the low-curvature direction. That conflict explains why full-gradient descent makes little progress in 2,000 steps. It does not imply that every neural-network loss has the same spectrum or that gradient descent is obsolete.

The original article compared one momentum setting, one RMSprop setting, and one Adam setting. A bad hyperparameter can make any of them look uniquely broken. I replaced that comparison with declared grids, each starting from zero and receiving 2,000 full-gradient steps:

method      settings tried   finite runs   best final loss
gradient descent       8           7          0.51200307
momentum              28          28          0.00009920
RMSprop               21          21          0.00009920
Adam                  28          28          0.00009920

The best configurations were:

gradient descent: lr=8e-5
momentum:         lr=1e-4, momentum=0.98
RMSprop:          lr=0.01, squared-gradient decay=0.999
Adam:             lr=0.01, beta1=0.9, beta2=0.999

Tuned momentum, RMSprop, and Adam all reach the noise floor of this synthetic fit. The earlier conclusions that momentum necessarily diverged and RMSprop inherently crawled came from one setting each. The grid does not rank the optimizers in general either. It tunes and evaluates on the same deterministic problem, uses a narrow hand-selected search, ignores per-step cost, and reports training loss rather than held-out behavior.

Adam is still worth following because its state appears in training memory budgets. The original Adam paper maintains an exponential average of gradients and an exponential average of squared gradients:

m[t] = beta1*m[t-1] + (1-beta1)*g[t]
v[t] = beta2*v[t-1] + (1-beta2)*g[t]²

Both start at zero, which biases early averages toward zero. Bias correction divides by the missing accumulated mass:

m_hat = m[t] / (1-beta1^t)
v_hat = v[t] / (1-beta2^t)

The update is:

weight -= alpha * m_hat / (sqrt(v_hat) + epsilon)

alpha remains a global step-size parameter. The denominator rescales coordinates according to their squared-gradient history, and m_hat contributes direction and momentum. Calling alpha/(sqrt(v_hat)+epsilon) the complete per-parameter learning rate omits the numerator and can misdescribe the actual step.

For P parameters, Adam commonly retains two state values per parameter in addition to parameters and gradients. If all four collections are fp32, the simple payload is:

parameters:       4P bytes
gradients:        4P bytes
first moment:     4P bytes
second moment:    4P bytes
total:           16P bytes

That arithmetic excludes activations, temporary buffers, allocator overhead, mixed-precision copies, and distributed sharding. Optimizer choice therefore changes both the update path and persistent memory, even when several methods reach the same loss on the two-parameter probe.

The pinned GPT-2 count makes 16P less abstract:

P = 124,439,808 parameter elements

one fp32 collection:
  124,439,808 * 4 bytes
  = 497,759,232 bytes
  = 0.464 GiB

parameters + gradients + two Adam moments:
  4 * 497,759,232 bytes
  = 1,991,036,928 bytes
  = 1.854 GiB

This is derived payload, not measured peak memory. It counts a tied parameter once because model.parameters() reports the unique parameter objects. It does not count optimizer bookkeeping objects, alignment, cached allocator blocks, or a temporary update buffer. Some optimizer implementations allocate moment state lazily on the first step, so memory immediately after model load can be lower than memory after the first update.

Mixed precision does not have one universal bytes-per-parameter number. One arrangement may keep fp16 weights for forward, fp16 or fp32 gradients, an fp32 master weight, and two fp32 moments. With fp16 gradients that simple list is:

fp16 working weights:   2P bytes
fp16 gradients:         2P bytes
fp32 master weights:    4P bytes
fp32 first moment:      4P bytes
fp32 second moment:     4P bytes
total:                 16P bytes

Using fp32 gradient buffers changes the total to 18P. Keeping no separate master copy changes it again. Quantized optimizers, factored state, or state offload add other representations and execution costs. A memory estimate must name which copies exist rather than attach one multiplier to the phrase “mixed precision.”

At larger model sizes, data-parallel replication puts these persistent collections on every rank. Sharding can partition gradients, optimizer state, and parameters so a rank owns only part of a collection, but the missing values then have to be communicated when an operation needs them. The memory problem becomes a scheduling and network problem rather than disappearing. The checkpoint probe addresses activation retention, which scales with batch and sequence shape. Sharding optimizer state addresses persistent parameter-related payload. They act on different terms in the peak.

thirty local slopes multiplied into almost nothing

The scalar engine can isolate depth without an optimizer. I built N repetitions of multiply then activation and measured the input gradient. With tanh and w=1:

tanh chain, w=1.0, dL/dx reaching the input:
  depth   1: dL/dx=7.864e-01
  depth  10: dL/dx=2.228e-01
  depth  50: dL/dx=3.396e-02
  depth 100: dL/dx=1.306e-02

Because w=1, each backward factor here is just 1-tanh², at most one in magnitude. The product reaches 0.013 at depth 100. With w=3, each local factor also includes the multiplication by 3, but tanh saturates close enough to one that its derivative dominates and the product collapses:

tanh chain, w=3.0 (units saturate):
  depth   5: dL/dx=8.204e-07
  depth  10: dL/dx=2.170e-14
  depth  20: dL/dx=1.518e-29
  depth  30: dL/dx=1.062e-44

At depth 6 the binary64 input gradient is 2.50e-8, below fp16’s smallest positive subnormal 5.96e-8, so this value rounds to zero when cast by the local probe.

ReLU provides a different control. Its local derivative is one on a positive branch and zero on a negative branch. With positive input and w=1, every selected branch has derivative one:

relu chain, w=1.0 (active path):
  depth   1: dL/dx=1.000e+00
  depth  50: dL/dx=1.000e+00
  depth 100: dL/dx=1.000e+00

The result remains one because this chain deliberately keeps every unit active and every weight equal to one. It does not test negative branches, changing activation patterns, matrices, or data. Raising every scalar weight to 1.2 makes the product grow:

relu chain, w=1.2:
  depth  10: dL/dx=6.192e+00
  depth  20: dL/dx=3.834e+01
  depth  50: dL/dx=9.100e+03

The exact result is 1.2^50 = 9,100.44. The probe has constructed an exploding product rather than discovered an emergent property of a trained network.

ReLU is not the universal replacement for tanh. The pinned GPT-2 path uses GELU, and many other architectures use other smooth or gated activations. Residual connections create an additional derivative route. If

y = x + F(x)

then locally

dy/dx = 1 + dF/dx

For vectors, the identity matrix is added to the Jacobian of F. This gives the upstream gradient a direct route around F, but it does not guarantee stable norms. The two routes can reinforce or cancel, and normalization and weight structure also matter. The scalar chains isolate products of local factors. They do not explain the training stability of a complete architecture by themselves.

the gradients vanished before the update

The depth experiment ended with a binary64 gradient of 2.50e-8. The number was nonzero, but casting it to fp16 produced zero. That is a second way for a gradient to vanish. The chain rule may deliver a signal and the chosen number format may then erase it.

I did these casts on the CPU. They establish what the formats can represent. They do not measure the speed, instruction selection, or accumulation behavior of an accelerator. The distinction matters because a format can exist in memory while an operator accumulates internally in another format.

Binary16, usually called fp16, has exact landmarks:

smallest positive normal:      2^-14 = 6.104e-5
smallest positive subnormal:   2^-24 = 5.960e-8
largest finite value:                 65504
gap from 1.0 to the next value:       9.766e-4

The normal range is where every represented value has the format’s full significand precision. Subnormal values extend the range towards zero by giving up leading significant bits. Zero is not reached at the smallest subnormal itself. With round to nearest, a positive value below half that subnormal rounds to zero. Values on either side of a midpoint also depend on the tie rule. Saying that everything below 5.960e-8 becomes zero would therefore be wrong.

The local casts make the boundary visible:

grad 1e-04:  fp16=1.000e-04
grad 1e-06:  fp16=1.013e-06
grad 1e-07:  fp16=1.192e-07
grad 1e-08:  fp16=0.000e+00

The value 1e-7 rounds to two subnormal steps. The value 1e-8 is below half of one step and rounds to zero. Long before that final rounding, representation error is already changing the gradient:

relative representation error, |cast(x) - x| / |x|

  x=1e-03: fp16 error=4.044e-04
  x=3e-05: fp16 error=6.288e-04
  x=1e-06: fp16 error=1.328e-02
  x=3e-08: fp16 error=9.868e-01

At 1e-6, the nearest fp16 number differs by 1.3 percent. At 3e-8, it rounds to one minimum subnormal, almost twice the source value. A gradient does not need to become zero before its magnitude is badly distorted.

The chain rule supplies a direct way to move these values. If the loss is replaced by

L_scaled = S * L

then every leaf gradient is multiplied by the same scale S:

dL_scaled/dw = S * dL/dw

The optimizer must divide by S before using the gradient. With S = 65536, the local cast gives:

1e-8 * 65536 -> fp16 6.552e-4
6.552e-4 / 65536     = 9.997e-9

The recovered result is not the original bit pattern. The fp16 rounding error remains. It is far smaller than the error from rounding the unscaled value to zero.

This technique also has an upper boundary. The scale multiplies large gradients along with small ones:

grad 0.1 * 65536 =   6553.6 -> fp16 6552
grad 1.0 * 65536 =  65536.0 -> fp16 inf
grad 2.0 * 65536 = 131072.0 -> fp16 inf

A fixed scale that rescues one tensor may overflow another. The small dynamic scaling probe starts at 65536, checks all scaled gradients for nonfinite values, discards an overflowing attempt, and halves the scale. For a batch containing gradients 2.0 and 1e-8, it settles at 16384:

big:  2.0  * 16384 -> 32768       -> divide -> 2.0000
tiny: 1e-8 * 16384 -> 1.638e-4    -> divide -> 9.997e-9

This is a mechanism demonstration, not a complete scaler. A production implementation also decides when to grow the scale after successful steps, how often to inspect it, and how skipped optimizer and scheduler steps interact. The mixed precision paper describes the fp16 method with an fp32 copy of the weights and loss scaling. Current framework behavior must be checked against the particular framework version, device, and operator.

Loss scaling protects a gradient while it passes through a narrow format. It does not solve the update at the other end. The local trajectory starts a weight at 1.0 and adds 0.0003 two thousand times:

pure fp16 weight after 2000 updates: 1.000000
fp32 master after 2000 updates:      1.600100
arithmetic result before rounding:   1.600000

Near 1.0, adjacent fp16 values are 0.0009765625 apart. The update is less than half that distance, so each individual fp16 addition rounds back to 1.0. The update itself did not underflow. It was lost when added to a much larger value, the same swamping mechanism measured in the earlier floating point experiment.

The fp32 trajectory differs from 1.600000 by about 0.0001 because two thousand fp32 additions round too. More precision reduces the accumulated error; it does not turn floating point into real arithmetic. Casting the fp32 weight to fp16 for a forward operation also rounds the value seen by that operation. The important difference is that the next small update is added to the fp32 state rather than to the rounded fp16 copy.

The spacing grows with magnitude because floating point gives approximately constant relative precision:

weight       next fp16 spacing     does +0.0003 change it?
  1.0             9.77e-4          no
  8.0             7.81e-3          no
 64.0             6.25e-2          no
512.0             5.00e-1          no

This probe intentionally repeats a constant positive update. An optimizer has signed, changing updates and may keep its own fp32 state. The result isolates one arithmetic failure. It does not predict that a whole model would remain frozen for exactly two thousand steps.

Bfloat16 makes a different trade. It retains eight exponent bits, like fp32, but only seven explicit fraction bits. Its smallest positive normal is 2^-126, about 1.175e-38. Its gap above 1.0 is 0.0078125, eight times the fp16 gap:

value       fp16 cast       bfloat16 cast
1e-8        0.000e+00       1.001e-08
1e-10       0.000e+00       1.000e-10

gap above 1.0
fp16        9.766e-4
bfloat16    7.812e-3
fp32        1.192e-7

The two tiny values survive the local bfloat16 casts because its normal range extends much farther towards zero. Its coarser significand means that a direct bfloat16 weight update has an even larger swamping problem near 1.0. Bfloat16 can remove the particular fp16 range pressure that motivated loss scaling while still requiring careful accumulation. Whether it improves an actual training run depends on accelerator support, kernel behavior, and which tensors remain in fp32. None of those performance claims were measured here.

the original disagreement had nowhere left to hide

The opening comparison mixed two different questions. The tape and PyTorch evaluate analytic derivative rules in floating point. The central difference estimates the same derivative by measuring a secant across a finite interval. At h=0.1, that interval was wide enough for curvature to appear in the fourth decimal. Reducing h initially removed truncation error, then subtraction and division amplified floating point error. The observed minimum near h=1e-5 was therefore a property of this expression, dtype, and machine, not a universal step size.

The investigation did reach a real attention block. On the pinned GPT-2 revision, the next token loss sent nonzero gradients through all thirteen recorded residual states and through query, key, and value slices in layers zero, five, and eleven. Moving a layer-zero query weight slice a small distance against the normalized gradient reduced the measured loss from 3.43444896 to 3.42128277; moving with it increased the loss to 3.44740677. Restoring the parameters restored the original loss exactly in that run. This is a local direction check, not evidence that a complete optimizer step improves a training objective.

The checkpoint probe measured saved tensor payload through PyTorch hooks, not peak process memory or GPU memory. The low precision probes measured CPU casts, not tensor core arithmetic. The optimizer grid covered one deterministic quadratic and the XOR sweep covered one tiny network. I did not train the transformer, run mixed precision on a GPU, or measure distributed gradient reduction.

Those limits no longer affect the first number. For L = tanh(a*b + a), the graph records exactly which local derivatives the backward pass needs. Reverse traversal accumulates every path into each leaf. The hand derivative, the tape, PyTorch, a central difference at its observed best interval, and a complex step on this analytic expression now converge on the same prediction. The value at h=0.1 was not a failed gradient. It was a measurement of a curve pretending, over an interval that was still too wide, to be a tangent.