In partnership with

SnackOnAI Engineering | Senior AI Systems Researcher | Technical Deep Dive | August 9, 2026

This issue dissects HIPS/autograd at version 1.9.1: the tracer, the Box and Node types, reverse and forward mode, the VJP rule library, and the failure modes you hit in production. Every number below came from running it on this machine, not from the README.

Covered: why operator overloading buys you while loops and recursion for free, exactly what the tape costs, where forward mode beats reverse mode, and three spellings of the same function that return three different gradients.

Excluded: GPU execution and JIT compilation, because autograd has neither by design. Excluded: a JAX tutorial, though JAX is the direct descendant and the lessons transfer to PyTorch eager mode unchanged.

What It Actually Does

autograd differentiates ordinary Python and NumPy. Not a DSL, not a graph you build and then run. You write the function, call grad, and get a function.

>>> import autograd.numpy as np
>>> from autograd import grad
>>> def tanh(x):
...     return (1.0 - np.exp(-2 * x)) / (1.0 + np.exp(-2 * x))
>>> grad(tanh)(1.0)
np.float64(0.419974341614026)

The trick is that autograd.numpy is not a reimplementation of NumPy. It is NumPy with every function wrapped in a decorator that watches for a special argument type. When it sees one, it records the call. When it does not, it calls straight through to real NumPy with zero overhead.

That is the entire idea. The rest is bookkeeping.

Written by Dougal Maclaurin, David Duvenaud, Matt Johnson, and Jamie Townsend at Harvard's HIPS group. Three of those four went on to build JAX. At 7.5k stars and 937 forks, it is now maintained by Agriya Khetarpal, Fabian Joswig, and Jamie Townsend, and it is still shipping: 1.9.1 supports NumPy 2.

The Architecture, Unpacked

Caption: Focus on the box marked "straight line only." Everything autograd can and cannot do follows from the fact that the tape records values flowing through primitives, never the Python that decided which primitives to call.

The Code, Annotated

The whole library, in one decorator

From autograd/tracer.py, lightly trimmed:

def primitive(f_raw):
    @wraps(f_raw)
    def f_wrapped(*args, **kwargs):
        boxed_args, trace, node_constructor, _ = find_top_boxed_args(args)
        if boxed_args:
            # Strip the boxes off, so f_raw sees plain ndarrays and runs at
            # full NumPy speed. autograd never reimplements the math.
            argvals = subvals(args, [(argnum, box._value) for argnum, box in boxed_args])
            parents = tuple(box._node for _, box in boxed_args)

            ans  = f_wrapped(*argvals, **kwargs)          # real work happens here
            node = node_constructor(ans, f_wrapped, argvals, kwargs, argnums, parents)
            return new_box(ans, trace, node)              # rebox, carrying the new node
        else:
            return f_raw(*args, **kwargs)   # ← THIS is the trick: untraced calls are FREE
    return f_wrapped

Caption: The else branch is why you can leave autograd.numpy imported in production code that never calls grad. No box in the arguments means one dict-free type check and a direct call to NumPy.

Two details reveal how hard that type check was optimized. Box uses __slots__ to avoid a per-instance dict. And the membership test carries a comment from the authors:

box_types = Box.types
isbox = lambda x: type(x) in box_types  # almost 3X faster than isinstance(x, Box)

Measured on this machine over two million iterations: type(x) in box_types took 25.4 ns, isinstance(x, Box) took 68.9 ns, a 2.71x gap. The comment is accurate. It is also a tell: this line runs once per argument per primitive call, so a 40 ns saving is worth writing down.

The tape is a list of closures, not a list of operations

class VJPNode(Node):
    __slots__ = ["parents", "vjp"]
    def __init__(self, value, fun, args, kwargs, parent_argnums, parents):
        self.parents = parents
        vjpmaker = primitive_vjps[fun]
        # ← THIS is the design decision. The VJP closure is built NOW, during
        #   the forward pass, capturing ans/args/kwargs. The backward pass is
        #   then just "call the closures in reverse order."
        self.vjp = vjpmaker(parent_argnums, value, args, kwargs)

def backward_pass(g, end_node):
    outgrads = {end_node: (g, False)}
    for node in toposort(end_node):
        outgrad = outgrads.pop(node)          # pop, so dead nodes free immediately
        ingrads = node.vjp(outgrad[0])
        for parent, ingrad in zip(node.parents, ingrads):
            outgrads[parent] = add_outgrads(outgrads.get(parent), ingrad)
    return outgrad[0]

Caption: Building the VJP eagerly is why higher-order derivatives are free. The backward pass is itself made of traced primitive calls, so tracing it again yields the second derivative with no extra machinery.

It In Action

Watching a second derivative get traced

autograd ships examples/print_trace.py, which swaps in a Node subclass that prints instead of differentiating. Same tracer, no autodiff. Running it on grad(grad(fun)) where fun(x) = sin(x + x) / 1:

A = 1.0
B = add(A,A)      = 2.0
C = sin(B)        = 0.9092974268256817
D = add(C,C)      = 1.8185948536513634
E = divide(D,2)   = 0.9092974268256817
F = cos(B)        = -0.4161468365471424      ← the BACKWARD pass starts here
G = multiply(1.0,F) = -0.4161468365471424    ← and it is being TRACED as forward ops
H = add(vspace,G,G) = -0.8322936730942848

Caption: Lines A through E are the function. Lines F through H are its derivative, recorded on the tape as if they were ordinary forward computation. That is the entire mechanism behind grad(grad(grad(f))).

Cost of each additional order, on tanh over a hundred points:

order

time (ms)

vs f(x)

value at x=1

0

0.0067

1.0x

0.7615941560

1

0.1139

17.1x

0.4199743416

2

0.3577

53.7x

-0.6397000084

3

1.0062

151.0x

0.6216266808

4

3.1792

477.1x

0.6650910448

5

9.3371

1401.3x

-5.5568935585

Each order costs roughly 3x the previous one, because each order traces the backward pass of the one below it. Not exponential in a scary sense, just a clean multiplicative constant you can plan around.

Logistic regression, end to end

Input: 2,000 samples, 20 features, synthetic labels from a known weight vector. Full-batch gradient descent using value_and_grad.

def loss(w):
    z = np.dot(X, w)
    return -np.mean(y * z - np.log(1 + np.exp(z))) + 1e-3 * np.dot(w, w)

w  = onp.zeros(D)
vg = value_and_grad(loss)      # one trace, both outputs, no duplicate forward pass
for i in range(300):
    L, g = vg(w)
    w = w - 1.0 * g

Output:

initial loss 0.693147
iter    1  loss 0.693147  |grad| 3.799197e-01
iter   10  loss 0.355080  |grad| 8.570764e-02
iter   50  loss 0.292544  |grad| 1.505422e-02
iter  300  loss 0.287923  |grad| 9.241711e-05
300 iterations in 118.1 ms  (0.394 ms per value_and_grad on 2000x20)
train accuracy 87.2%   cosine(w, w_true) = 0.9959
max |autograd - central difference| = 1.155e-10
check_grads(order=2) passed

The gradient agrees with central differences to 1.155e-10, which is the finite-difference method's error, not autograd's. Autograd's answer is exact to floating point.

Why exactness is the smaller half of the win

Reverse mode gets the gradient of a scalar loss with respect to all parameters in a constant multiple of one forward pass, regardless of how many parameters there are. Finite differences need two evaluations per parameter. Measured on the same logistic loss at four sizes:

parameters

autograd grad (ms)

equivalent finite differences (ms)

speedup

error

10

0.299

0.4

1x

2.20e-11

100

0.259

6.0

23x

6.76e-11

1,000

0.596

370.8

622x

4.96e-11

5,000

1.943

7,782.0

4,006x

3.85e-11

Caption: At ten parameters autograd is a wash. At five thousand it is four thousand times faster and more accurate. This single table is the reason gradient-based ML is possible at all, and it is the practical content of the Baydin et al. survey.

Why This Design Works, And What It Trades Away

Why it works. The tape records primitives, not Python. autograd therefore never needs to understand your if, your while, your recursion, or your closures. It just watches which primitives actually executed.

Newton's method with a data-dependent stopping condition, differentiated:

def newton_sqrt(x):
    g = x / 2.0
    for _ in range(60):
        if abs(getval(g * g - x)) < 1e-14:   # getval peeks through the box
            break
        g = (g + x / g) / 2.0
    return g

d/dx sqrt(  2.0) = 0.353553390593   analytic = 0.353553390593
d/dx sqrt(  9.0) = 0.166666666667   analytic = 0.166666666667
d/dx sqrt(100.0) = 0.050000000000   analytic = 0.050000000000

Exact to twelve digits, through an iteration count that varies with the input. Recursion works too: grad of a recursively defined x**8 gives 8x^7 correctly. This is what makes differentiating a fluid simulation possible.

What it trades away.

The overhead is per primitive call, so it never amortizes. Measured cost of grad(f) relative to f:

workload

primitives

f (ms)

grad f (ms)

ratio

matmul 100x100

3

0.0084

0.1087

12.95x

matmul 400x400

3

0.0516

0.2503

4.85x

matmul 1000x1000

3

0.3482

0.8083

2.32x

Python loop, ten iterations

30

0.0468

0.6988

14.93x

Python loop, hundred iterations

300

0.4379

6.1968

14.15x

Python loop, thousand iterations

3,000

4.4588

58.8052

13.19x

Caption: Top three rows, same three primitives, arrays growing. The ratio collapses toward 2.3x as real FLOPs swamp the bookkeeping. Bottom three rows, primitive count growing. The ratio refuses to move off roughly 14x. That is the whole performance story in one table.

There is no cache. Every call to grad(f) re-traces from scratch. A twenty-layer network, same input, same shapes, call one and call thirty cost the same 1.167 ms against a 0.199 ms forward pass. No compilation, no shape specialization, no fusion. This is precisely the hole jax.jit was built to fill.

The tape is your memory bill. At 150,000 traced primitives the tape cost 30.3 MB of resident memory and 6.83 seconds to build and unwind. That is roughly 45 microseconds and 200 bytes per primitive on a trivial operation.

Mutation is out. y[0] = 5.0 on a boxed array raises TypeError: 'ArrayBox' object does not support item assignment. So does float(box). Both are correct: a mutable tape entry cannot be replayed, and a float has nowhere to store a node.

Technical Moats

The tracer is not the moat. It is 195 lines and you could rewrite it in a weekend.

The moat is numpy_vjps.py: 975 lines of vector-Jacobian product rules, plus 303 more for the forward-mode versions and 351 for linear algebra. Every one of them has to be correct under broadcasting, correct for complex inputs, correct at boundary values, and correct when composed with itself for higher-order derivatives.

unbroadcast_f alone is the kind of thing you only get right by getting it wrong. When you write x + b with x of shape (1000, 20) and b of shape (20,), the gradient flowing back to b must be summed over the broadcast axis. Miss it and you get a shape error, or worse, a silently wrong gradient that still trains.

The second moat is check_grads, which verifies any function against finite differences to arbitrary order, including the mixed forward-reverse paths. check_grads(loss, modes=['rev'], order=2) passing is a meaningfully strong statement about the correctness of every rule on the path.

The third is the vspace abstraction, which lets gradients flow through dicts, lists, and tuples of arrays without any special casing in the tracer. That is why you can define a neural network's parameters as a nested list and grad returns the same nested shape back.

Insights

Insight One: "Autograd is slow" is the wrong sentence

The correct sentence is that autograd charges a fixed toll per primitive call and nothing per FLOP. Reread the table above. On a 1000x1000 matmul the gradient costs 2.32x the forward pass, which is close to the theoretical floor for reverse mode. On a Python loop doing tiny operations it costs 13x to 15x and stays there no matter how long the loop gets.

The implication is not "pick a different AD library." It is "call fewer, bigger primitives." Vectorize the loop and the overhead disappears into the FLOPs.

That lesson is not autograd-specific. PyTorch eager mode has the identical architecture, a tape of recorded operations built by operator overloading, and the identical failure mode: a training loop full of small tensor ops spends most of its time in the interpreter and dispatcher. Every engineer who has watched a PyTorch model go three times faster under torch.compile without changing a line of math has met this exact tax. autograd is just the version small enough to read.

Insight Two: the "no longer actively developed" line everyone quotes is stale, and the new work has teeth

Search for autograd and you will find the note that the main developers moved to JAX. It is repeated in issue threads and blog posts as though the project stopped in 2020.

Version 1.9.1 ships NumPy 2 support, requires-python >= 3.10, and an __array_ufunc__ implementation on ArrayBox under three maintainers who were not the original authors.

That last item quietly changed the contract. The README's first instruction has always been to import autograd.numpy instead of numpy. Test it now and plain NumPy mostly works: numpy.sum(numpy.tanh(x)) on a boxed array produces a correct gradient, identical primitives, identical timing, because ufunc dispatch routes it back into autograd's wrappers.

Mostly. Which is worse than never.

Takeaway

Three ways to write the same function. Three different gradients. Same point.

Take f(x) = sum(clip(x, 0, 2)) evaluated at x = [0, 1, 2], where two of the three inputs sit exactly on a bound.

autograd.numpy.clip(x, 0, 2)                  →  [0.0, 1.0, 0.0]
numpy.clip(x, 0, 2)          (via dispatch)   →  [1.0, 1.0, 1.0]
anp.minimum(anp.maximum(x, 0), 2)             →  [0.5, 1.0, 0.5]

All three compute the same values. All three are defensible subgradients at a kink. None of them agree.

The cause is two lines of rule source. clip's VJP zeros the gradient wherever the output touches a bound, which catches inputs that were already at the bound and never actually clipped:

defvjp(anp.clip, lambda ans, x, a_min, a_max:
       lambda g: g * anp.logical_and(ans != a_min, ans != a_max))

maximum and minimum instead split the gradient evenly at exact ties:

def balanced_eq(x, z, y):
    return (x == z) / (1.0 + (x == y))   # ← 0.5 when x and y tie

And routing through plain numpy.clip decomposes into elementwise primitives that take the pass-through convention. Three code paths, three conventions, zero documentation of the difference.

This is not a bug report. Subgradient choice at a kink is genuinely arbitrary and every framework picks one. The point is that it is arbitrary per rule, invisible at the call site, and reachable by refactoring that looks purely cosmetic. If your parameters ever initialize at exactly zero and pass through a clip(w, 0, None), autograd hands you a zero gradient and your weights never move. The spelling determined the training dynamics.

Check check_grads on your own kinks before you trust a refactor. Finite differences will not save you here, because they disagree at kinks too.

TL;DR For Engineers

  • Cost of grad(f) is 2.32x on a 1000x1000 matmul and 13x to 15x on a loop of tiny ops, regardless of loop length. Overhead is per primitive, never per FLOP. Vectorize or pay forever.

  • Reverse mode beats finite differences by 1x at ten parameters, 622x at a thousand, 4,006x at five thousand, and is exact rather than approximate.

  • Forward mode is not a curiosity. For R^1 → R^200, forward computed the Jacobian in 0.063 ms against reverse mode's 4.804 ms, a 76x win. The crossover is inputs versus outputs, not preference.

  • Every grad call re-traces. No cache, no fusion, no shape specialization. That absence is exactly the gap jax.jit and torch.compile fill.

  • Higher-order derivatives cost about 3x per order and come free architecturally, because the backward pass is itself traced.

Explain It Like I'm New

Almost all machine learning is one loop: guess some numbers, measure how wrong you are, nudge the numbers in the direction that reduces the error, repeat. The nudge direction is called the gradient, and finding it is the hard part.

You could find it by trial and error. Wiggle one number, see if the error goes down, put it back, wiggle the next. That works, and for a handful of numbers it is fine. For a model with a billion numbers it would take longer than the universe has existed.

Automatic differentiation is the alternative, and it is neither trial and error nor the symbolic algebra you did in calculus class. The idea is that any program, however complicated, is ultimately a long chain of simple operations: add, multiply, take a logarithm. Calculus already tells us the derivative of each of those. The chain rule tells us how to stitch them together.

So the computer watches your program run, writes down every simple operation in order, then walks that list backwards applying the chain rule. The cost of the whole backward walk is a small constant multiple of running your program once, no matter how many numbers you are adjusting. That is the entire reason modern AI is affordable.

Autograd's contribution was making this feel like nothing. You write ordinary Python. You call one function. You get the derivative. No special syntax, no graph to build, no restrictions on loops or branches.

The reason this matters is that the same idea has since escaped machine learning. Physics simulators, financial models, robotics controllers, and drug design pipelines are all being rewritten so their gradients can be taken, which turns "simulate this" into "optimize this."

See It In Action

  • Automatic Differentiation, Matt Johnson (Deep Learning Summer School Montreal, MILA). One of autograd's authors builds a working autodiff system live. The clearest available explanation of why the tape holds closures rather than operations.

  • examples/print_trace.py (HIPS/autograd). Not a video, the single most useful file in the repo. Fifty lines that reuse the tracer with a printing Node instead of a differentiating one, proving the tracer knows nothing about calculus.

  • examples/fluidsim/fluidsim.py (HIPS/autograd). Backpropagation through a full fluid simulation to optimize an initial condition. The best argument that "differentiate arbitrary Python" is a real capability and not a slogan.

  • Dougal Maclaurin's PhD thesis, chapter four (dougalmaclaurin.com, Harvard). The design rationale from the person who wrote the first version, including how the primitive VJPs are defined and why complex numbers needed special handling.

Community Conversation

  • The stale maintenance note (python-control issue). A representative thread where maintainers of another library conclude autograd is abandoned, based on the README line about the authors moving to JAX. Useful as evidence of how long a single sentence can misdirect an ecosystem.

  • Dougal Maclaurin on autograd's lineage (dougalmaclaurin.com). Notes that autograd was ported to Lua and Julia and helped inspire PyTorch, which borrowed the name for its autodiff module. The clearest first-person account of how one Harvard side project shaped three frameworks.

  • Where AD should be going (van Merriënboer et al., arXiv:1810.11530). Argues that tape-based operator overloading, autograd's approach, blocks ahead-of-time optimization, and proposes a graph IR supporting closures and recursion so differentiation can happen by source transformation instead. Reads as a direct critique of the architecture dissected above.

  • Baydin and Pearlmutter's original pitch (arXiv:1404.7456). Written when ML had largely not noticed AD existed, and argues reverse mode both predates and generalizes backpropagation. Worth reading for how recently this was a fringe position.

Read The Small One To Understand The Big Ones

autograd is the right tool for a narrower set of jobs than it was in 2015. If you need GPUs, use JAX or PyTorch. If you need throughput on small operations, you need a compiler, and autograd deliberately does not have one.

It remains the right tool for gradients of scientific Python that already runs on CPU, for research code where retracing every call is a feature because the graph genuinely changes, and for any situation where you need to read the autodiff implementation to trust it.

It is also, and this is the underrated part, the fastest way to actually understand the tools you do use in production. The 13x tax on small operations, the absence of a trace cache, the subgradient conventions that differ per rule: every one of those exists in PyTorch eager mode and in JAX before jit. In autograd you can see the line of code responsible.

Two hundred lines. One afternoon. Then go read your profiler again.

References

autograd differentiates arbitrary Python by wrapping NumPy primitives in a decorator that records a tape of closures, which is why loops, branches, and recursion work for free and why the overhead is fixed per operation rather than per FLOP. Measured here, that means 2.32x on a large matmul and a stubborn 13x to 15x on many small ones, a distinction that determines whether the library fits your workload. The same architecture underlies PyTorch eager mode, making autograd the smallest readable version of a tax most engineers pay without seeing it.

Sponsored Ad If you enjoy practical AI insights, check out SnackOnAI and support the newsletter by subscribing, sharing, and exploring our sponsored ad, it helps us keep building and delivering value 🚀

Own Search With Podcasts

Your competitors are fighting over the same keywords. The smartest brands are building the authority that search engines, AI platforms, and customers trust everywhere.

Every relevant podcast appearance can produce branded mentions, backlinks, transcripts, citations, clips, expert content, and third-party proof that keeps compounding across search and AI discovery.

PodPitch searches millions of podcasts, finds the shows that matter to your market, develops the angle, sends personalized pitches, and follows up automatically until your experts are booked.

Growth teams are already using podcast appearances to build distributed authority that cannot be manufactured by publishing another generic SEO article.

Only 20 SEO, AEO, and GEO demo spots are available this month. Once they’re claimed, the offer disappears.

Start building searchable authority now, before your competitors own the conversations shaping your market.

Appear on 3,853,234 Podcasts

Recommended for you

View all
caret-right
© 2026 Snack On AI.
beehiivPowered by beehiiv