Micrograd: Build and Test a Scalar Autograd Engine

0 · The mathematical core, not framework identity

Karpathy's micrograd is an educational scalar reverse-mode engine. Its key ideas are values, recorded operations, local pullbacks and accumulation. A lecture's view count, invented recording scene or an unmeasured10,000× speedup is not needed to explain them. PyTorch, JAX, TensorFlow, tinygrad and llm.c share chain-rule mathematics but not one identical graph/closure architecture.

Our goal is a complete inspectable engine supporting addition, multiplication, constant powers, tanh, exp, log and ReLU, then a3→4→4→1 network on four local points. The program has explicit power/log domains and single-use backward semantics. It is a teaching implementation, not a replacement for PyTorch's version tracking, tensor broadcasting, device dispatch or higher-order AD.

Prerequisites and outcomes

Bring the chain rule and manual backprop. By the end, explain each stored field, derive each local closure, defend reverse-topological order on a shared graph, fit the four-point network and numerically check softmax cross entropy. Optional pre-reading: the lecture, upstream engine, upstream network. This variant uses tanh outputs; upstream's current network source uses a different activation arrangement.

1 · The Value class — the atom

A node stores data, grad, _prev and _backward; _op labels diagrams. Leaves have no producer and a no-op pullback. The fishing-net metaphor describes multiple paths transmitting signed derivatives, not physical forces. Scalars need no broadcast reduction; a tensor engine must add that rule.

The complete class follows so the method fragments do not accidentally remain unattached top-level functions. Run the fences in order. Local derivatives are captured from forward scalar values; do not mutate parameters until after backward. We intentionally reject repeated backward on the same result instead of reusing stale intermediate accumulators.

import math
import random
 
class Value:
    def __init__(self,data,_children=(),_op=''):
        self.data=float(data)
        self.grad=0.
        self._prev=tuple(_children)
        self._op=_op
        self._backward=lambda:None
        self._used=False
 
    def __repr__(self):
        return f'Value(data={self.data:.6g}, grad={self.grad:.6g})'
 
    @staticmethod
    def wrap(x):
        return x if isinstance(x,Value) else Value(x)
 
    def __add__(self,other):
        other=Value.wrap(other)
        out=Value(self.data+other.data,(self,other),'+')
        def back():
            self.grad+=out.grad
            other.grad+=out.grad
        out._backward=back
        return out
 
    def __mul__(self,other):
        other=Value.wrap(other)
        a,b=self.data,other.data
        out=Value(a*b,(self,other),'*')
        def back():
            self.grad+=b*out.grad
            other.grad+=a*out.grad
        out._backward=back
        return out
 
    def __pow__(self,n):
        if not isinstance(n,(int,float)):
            raise TypeError('constant numeric exponent required')
        a=self.data
        if (a<0 and not float(n).is_integer()) or (a==0 and n<1 and n!=0):
            raise ValueError('power outside differentiable real domain')
        out=Value(a**n,(self,),f'**{n}')
        slope=0. if n==0 else n*a**(n-1)
        def back(): self.grad+=slope*out.grad
        out._backward=back
        return out
 
    def __neg__(self): return self*-1
    def __sub__(self,other): return self+-Value.wrap(other)
    def __rsub__(self,other): return Value.wrap(other)+-self
    def __radd__(self,other): return self+other
    def __rmul__(self,other): return self*other
    def __truediv__(self,other): return self*Value.wrap(other)**-1
    def __rtruediv__(self,other): return Value.wrap(other)*self**-1
 
    def tanh(self):
        t=math.tanh(self.data)
        out=Value(t,(self,),'tanh')
        def back(): self.grad+=(1-t*t)*out.grad
        out._backward=back
        return out
 
    def exp(self):
        e=math.exp(self.data)
        out=Value(e,(self,),'exp')
        def back(): self.grad+=e*out.grad
        out._backward=back
        return out
 
    def log(self):
        a=self.data
        if a<=0: raise ValueError('log requires a positive input')
        out=Value(math.log(a),(self,),'log')
        def back(): self.grad+=out.grad/a
        out._backward=back
        return out
 
    def relu(self):
        a=self.data
        out=Value(max(0.,a),(self,),'relu')
        def back(): self.grad+=(a>0)*out.grad
        out._backward=back
        return out
 
    def topo(self):
        order=[];seen=set()
        def visit(v):
            if v in seen:return
            seen.add(v)
            for parent in v._prev:visit(parent)
            order.append(v)
        visit(self)
        return order
 
    def backward(self):
        order=self.topo()
        if any(v._used for v in order if v._prev):
            raise RuntimeError('rebuild the forward graph before another backward')
        self.grad=1.
        for v in reversed(order):v._backward()
        for v in order:
            if v._prev:v._used=True

2 · Operations — the backward closure contract

2.1 Add

For out=a+b, both derivatives are1. Each input receives out.grad; if a and b are the same object, both contributions must still be added. Deduplicating traversal nodes must not deduplicate the operation's two operand edges.

2.2 Multiply

For out=a*b, a receives b_forward*out.grad and b receives a_forward*out.grad. The closure reads out.grad, not its parent's current gradient. Forward operand values are captured because mutating them before backward would otherwise change the derivative of the computation we recorded.

2.3 Power

For constant n, d(a^n)/da=n*a^(n−1) where the real derivative exists. The n=0 branch has zero slope, including a=0 under this program's constant-function convention. Reject negative bases with noninteger powers and zero raised to negative or singular subunit powers. This implementation does not differentiate with respect to the exponent.

2.4 Mixed arithmetic

__radd__, __rmul__, __rsub__ and __rtruediv__ support 2+x, 2*x, 2-x and 2/x. A numeric target can remain constant: wrapping every target as trainable is unnecessary. Converting a prediction to float(pred.data) before constructing the loss would break its path.

2.5 Tanh

Cache t from forward and use 1−t². It is smooth and bounded; saturation can make the gradient small. It does not imply successful optimization from every initialization.

2.6 Exp and log

Exp's derivative is its output; log's is the reciprocal of its positive input. math.exp can overflow, so shift logits before a softmax/log-sum-exp composition. ReLU uses slope1 for positive inputs and0 otherwise, a chosen subgradient at zero.

3 · Backward — wait for every consumer

A node cannot propagate its full derivative until all consumers have contributed. Postorder DFS builds producer-before-consumer order; reversing it visits consumers first. The traversal costs O(V+E) time and O(V) bookkeeping. Python recursion bounds make this implementation unsuitable for arbitrarily deep scalar chains.

A useful counterexample must have a shared intermediate, not merely a reused leaf whose pullback does nothing. Take a=x*x, b=2*a, c=3*a, L=b+c. At x=3, dL/dx=10x=30. A naive recursive backward that marks a visited after the b path propagates only2 through a, then never propagates c's3. It yields12, or18 if branches reverse. The original x*x+x example does not reliably expose that traversal bug because x is a leaf.

x=Value(3.)
a=x*x;b=2*a;c=3*a;L=b+c
L.backward()
assert L.data==45 and x.grad==30
x2=Value(3.);a2=x2*x2;b2=2*a2;c2=3*a2;L2=b2+c2
seen=set();L2.grad=1
# Deliberately wrong: release a producer before all consumers arrive.
def wrong_walk(v):
    if v in seen:return
    seen.add(v);v._backward()
    for p in v._prev:wrong_walk(p)
wrong_walk(L2)
assert x2.grad!=30
u=Value(3.);(u*u).backward()
assert u.grad==6
# Assignment would retain only one of the two contributions: 3, not 6.
assert 3!=u.grad
try:L.backward()
except RuntimeError:pass
else:raise AssertionError('repeated backward should reject stale graph')

4 · A scalar neuron — check the arithmetic before a network

This is a squared-error scalar neuron, not the nine-parameter BCE network from DL010. Let w=.5,x=1,b=.1,target1; prediction.6, residual−.4 and loss.16. Derivatives are dw=−.8, db=−.8, dx=−.4 and dtarget=.8 if target is tracked.

w,x,b,target=Value(.5),Value(1.),Value(.1),Value(1.)
loss=((w*x+b)-target)**2
loss.backward()
assert math.isclose(loss.data,.16)
assert math.isclose(w.grad,-.8) and math.isclose(b.grad,-.8)
assert math.isclose(x.grad,-.4) and math.isclose(target.grad,.8)
q=Value(2.)
((2-q)+(3/q)).backward()
assert math.isclose(q.grad,-1.75)
q0=Value(0.);(q0**0).backward();assert q0.grad==0

5 · Neuron, Layer and MLP

A Neuron forms a dot product plus bias; a Layer groups neurons; an MLP composes layers. All neurons below use tanh, including the output, matching the ±1 squared-error toy. A logits classifier would make the last layer linear and use the corresponding loss. Shape checks avoid silent truncation by zip.

class Neuron:
    def __init__(self,nin,rng):
        self.w=[Value(rng.uniform(-1,1)) for _ in range(nin)]
        self.b=Value(0)
    def __call__(self,x):
        if len(x)!=len(self.w):raise ValueError('neuron input width')
        return sum((wi*xi for wi,xi in zip(self.w,x)),self.b).tanh()
    def parameters(self):return self.w+[self.b]
 
class Layer:
    def __init__(self,nin,nout,rng):
        self.neurons=[Neuron(nin,rng) for _ in range(nout)]
    def __call__(self,x):return [n(x) for n in self.neurons]
    def parameters(self):return [p for n in self.neurons for p in n.parameters()]
 
class MLP:
    def __init__(self,nin,nouts,seed=0):
        rng=random.Random(seed);sizes=[nin]+list(nouts)
        self.layers=[Layer(a,b,rng) for a,b in zip(sizes,sizes[1:])]
    def __call__(self,x):
        for layer in self.layers:x=layer(x)
        return x[0] if len(x)==1 else x
    def parameters(self):return [p for layer in self.layers for p in layer.parameters()]

5.1 Train it on four points

This network has (3*4+4)+(4*4+4)+(4*1+1)=41 parameters. Rebuild the graph after each update; leaf parameter gradients accumulate until cleared. The data are four three-feature examples, not a two-moons dataset. Sum squared error is used, so switching to a mean changes the scale by four.

model=MLP(3,[4,4,1],seed=0)
Xs=[[2.,3.,-1.],[3.,-1.,.5],[.5,1.,1.],[1.,1.,-1.]]
Ys=[1.,-1.,-1.,1.]
assert len(model.parameters())==41
history=[]
for step in range(200):
    predictions=[model(x) for x in Xs]
    loss=sum((p-y)**2 for p,y in zip(predictions,Ys))
    history.append(loss.data)
    for p in model.parameters():p.grad=0
    loss.backward()
    for p in model.parameters():p.data-=.05*p.grad
final_predictions=[model(x).data for x in Xs]
final_loss=sum((a-b)**2 for a,b in zip(final_predictions,Ys))
assert final_loss<history[0]
assert all(a*b>0 for a,b in zip(final_predictions,Ys))
print('four-point scalar training',history[0],final_loss,final_predictions)

These assertions state the observed seed-specific contract without copying invented loss logs. Fit to four points does not estimate unseen-data performance.

6 · How this differs from PyTorch

1. Granularity: one Python node per scalar operation versus tensor-level kernels and broadcast-aware pullbacks. A tensor matmul replaces a large scalar graph but still needs shape-correct derivatives.
2. Runtime: Python objects/closures versus dispatch, optimized CPU/GPU kernels and engine scheduling. Speed depends on workload and hardware; no measured factor is implied here.
3. Semantics: PyTorch offers version counters, device/dtype management, higher-order differentiation, graph-retention controls and many operations. This engine rejects graph reuse, has no tensor views and has no second-derivative graph.

The mathematical VJP chain rule is shared. Eager PyTorch does not literally instantiate this Value class, and compiled systems need not store Python closures. Retaining scalar loss objects here retains their reachable graph; after PyTorch backward, saved tensors are normally freed even if some graph structure remains referenced. Log numeric loss values to avoid unnecessary ownership.

7 · Two exercises — implemented and checked

ReLU and stable softmax cross entropy

The class already includes ReLU/log. For logits v, choose m=max(data) as a numeric shift and compute log(sum(exp(v-m)))−(v_target−m). Treating m as a constant is valid because the common shift cancels algebraically; it avoids differentiating a max kink unnecessarily. A fused implementation would be preferable for wider tensors.

def softmax_ce(logits,target):
    if not logits or not 0<=target<len(logits):raise ValueError('target')
    m=max(v.data for v in logits)
    return sum((v-m).exp() for v in logits).log()-(logits[target]-m)
 
values=[1000.,1001.,999.]
v=[Value(a) for a in values];ce=softmax_ce(v,1);ce.backward()
ex=[math.exp(a-max(values)) for a in values]
expected=[e/sum(ex)-(i==1) for i,e in enumerate(ex)]
assert all(math.isclose(a.grad,g,abs_tol=1e-12) for a,g in zip(v,expected))
for i in range(3):
    plus=values.copy();minus=values.copy();plus[i]+=1e-4;minus[i]-=1e-4
    n=(softmax_ce([Value(a) for a in plus],1).data-softmax_ce([Value(a) for a in minus],1).data)/2e-4
    assert abs(n-v[i].grad)<1e-7
for a,g in [(-2.,0.),(0.,0.),(2.,1.)]:
    z=Value(a);z.relu().backward();assert z.grad==g

Draw the computation graph without installing Graphviz

Emit DOT text and inspect edges. Stable local IDs come from topological order. Keep duplicate operand edges: x*x has two uses even though traversal visits x once. Rendering DOT is optional and not required for the numerical test.

def graph_dot(loss):
    nodes=loss.topo();ids={v:i for i,v in enumerate(nodes)}
    lines=['digraph G {']
    for v in nodes:
        lines.append(f'n{ids[v]} [label="{v._op or "leaf"}: {v.data:.4g} / grad {v.grad:.4g}"];')
        for p in v._prev:lines.append(f'n{ids[p]} -> n{ids[v]};')
    return '\n'.join(lines+['}'])
assert 'digraph G' in graph_dot(ce)

8 · Reading an autodiff system without forcing one architecture

Ask four questions: which values are differentiable, how are primitive derivative rules represented, how are dependency order and saved values managed, and how are cotangents summed/reduced? Search for the roles, not literal _prev fields in every framework.

8.1 Micrograd, PyTorch, JAX, tinygrad and llm.c

• Micrograd: dynamic scalar graph with local Python pullbacks.
• PyTorch eager: tensor operations create autograd Function nodes; torch.compile can capture supported regions, use AOTAutograd and compile graphs.
• JAX: function transformations use tracing and intermediate representations; grad, jvp, vjp and jit are distinct composable transformations. AD is not merely a computer algebra system producing a handwritten symbolic expression.
• tinygrad: its tagged source describes Tensor operations and a lazy/compiled execution stack. Locate its differentiation and scheduling machinery in that version rather than presuming a legacy Function class.
• llm.c: the GPT-2 C implementation contains explicit forward/backward kernels and parameter/activation buffers. It does not discover an arbitrary user graph with this generic closure engine.

For f:R^n→R^m, a forward-mode JVP propagates one input direction; reverse mode propagates one output cotangent. A full Jacobian can take n JVPs or m VJPs, favoring reverse when m is small relative to n, subject to memory/cost. Micrograd is reverse-mode for a scalar loss. A vector loss can be weighted and summed to seed an arbitrary cotangent; that is different from constructing its entire Jacobian.

8.2 Dynamic, static and captured execution

Dynamic recording follows executed branches, but cannot differentiate an unexecuted branch. Static IRs permit global optimization and explicit control-flow operators; they are not universally worse to debug. Capturing eager code can enable fusion while retaining an eager reference, but guards, graph breaks, dynamic shapes and compilation latency require validation. Eager execution can still use fused library kernels, so “dynamic cannot fuse” is false. Prefer the simplest system that matches the task; do not claim every researcher must use one style.

9 · Failure cases

Assignment instead of accumulation: x*x at3 must get6, not3. Wrong traversal: the diamond in §3 must get30, not12 or18. Stale gradients: rebuild graphs and clear parameter accumulators between optimizer steps, not between intended contributions to the same step. Numeric constants: they are harmless if mixed operators wrap them; extracting prediction data is the graph break. Mutation: cached scalar derivatives in this variant refer to the forward snapshot, but production tensor mutation needs stricter version/alias controls. Retained logs: store loss.data, not every Value root, when only numbers are needed.

Changing control-flow length changes the executed tape. A symbolic system can also represent branches or powers, so this is not something “impossible for symbolic derivatives.” The distinction is what this engine records and evaluates, not what mathematics can express.

10 · Retention scaffold — worked answers

1. Fields: data is the forward scalar; grad is accumulated sensitivity; parents define dependencies; the closure distributes output sensitivity using local derivatives. _op is optional labeling.
2. Multiply: add b_forward*out.grad to a and a_forward*out.grad to b, even when a is b.
3. Order: wait for every consumer; the shared-a diamond proves that early DFS release loses a contribution below a.
4. Why +=: x*x at3 has two equal paths of3, totaling6. A single assignment cannot represent both.
5. Framework difference: scalar versus tensor granularity, Python versus optimized dispatch, and deliberately reduced semantics; reverse-mode chain-rule mathematics is shared, not every internal class.
6. Neuron: L=.16,dw=db=−.8,dx=−.4 for the stated squared-error fixture.
7. Extensions: §7 defines shifted log-sum-exp CE, derives P-one_hot, numerically checks all three coordinates, and checks ReLU on negative/zero/positive inputs.
8. AD modes: n input-basis JVPs versus m output-basis VJPs for a full Jacobian; reverse fits a many-parameter scalar training loss. Neither computes arbitrary second derivatives for free.

Rebuild the class tomorrow, then recheck the shared graph before adding operations. DL013 replaces scalar work with arrays and measures the result rather than promising100× speed.

Sources and execution boundary

Upstream engine, network, PyTorch2.9 mechanics, JAX v0.6.2 advanced autodiff, tinygrad v0.11.0, llm.c, Baydin et al.. Exact dated source snapshots bind mutable upstream URLs to hashes. This lesson's scalar CPU program was tested; no JAX/tinygrad/compiler/large-model benchmark is claimed.