Generating an image with a diffusion or flow-matching model typically involves integrating a velocity field for hundreds of sequential steps. That computation suits a GPU poorly and an analog neuromorphic substrate well, where one fast, low-power physical relaxation could replace the whole sequence. Two things are needed: a generative algorithm whose computation is a relaxation, and a way to train it without backpropagation, whose global gradient transport is hard to realise in analog. Energy Matching supplies the first, sampling by gradient-descent relaxation on a time-independent potential - but obtaining the gradient that descent follows is itself a backward pass, once per integration step. We prove that a predictive coding network with a single linear injection term added to its energy relaxes to a state whose visible layer carries that gradient inherently, with an error that vanishes for small injection strengths. An Energy Matching model trained by backpropagation can therefore be loaded into this network unchanged and sampled by relaxation alone, matching the original model's FID to within 0.04 at CIFAR-10 and ImageNet-32 scale. The network can also be trained without backpropagation: Equilibrium Propagation reproduces backpropagation's parameter updates to a cosine similarity of 0.9998, and training end to end with those updates reaches an FID of 3.64 on unconditional CIFAR-10 against 3.54 for backpropagation. To our knowledge this is the first CIFAR-10 generative model trained with Equilibrium Propagation, and the first at competitive quality by any backpropagation-free rule. This quality does not require group normalisation or attention: an architecture of only convolutions and pointwise nonlinearities, both native to an analog crossbar, reaches an FID of 3.51 when trained by backpropagation.
Figures & tables
Figure 1: The same network, read two ways. (a) Energy Matching generates by x˙=−∇xV , and each evaluation of ∇xV is a backward pass through the network ( Jk=∂fk/∂hk−1 ). (b) In the PotentialPCN the same blocks fk make predictions, but the hidden states (shaded) are not computed from them: instead they are free variables that relax by descending EPCN to stationary values h∗ . Each error unit ek=fk(hk−1)−hk compares a prediction against a state, and the energy is the sum of squared errors plus a single linear term γαhN , which at hN∗ equals γV up to O(γ2) . The same term fixes eN∗=γα , and the recursion ek∗=Jk+1⊤ek+1∗ carries it down through the same transposed Jacobians as (a), so the gradient of EPCN at the visible node is γ∇xV+O(γ2) . Dividing by γ , the visible dynamics x˙=−γ−1∇xEPCN reproduce those of (a) to O(γ) : the backward pass has been performed by the relaxation.
Figure 2: Relative deviation ϵ=∥vPCN−vFFN∥/∥vFFN∥ of the relaxed PCN velocity from the exact FFN autograd velocity, on the CIFAR-10 replication checkpoint of Table 1 , evaluated at 512 flow-matching interpolants xt constructed from training images exactly as during training. The deviation is O(γ) (dashed) until it meets the arithmetic floor of the number format; float64 meets no floor in the measured range, and at the operating γ=10−4 (dotted) the deviation lies below the TF32 floor of our runs.
FID ↓
Dataset
Checkpoint
Feedforward
PCN relaxation
CIFAR-10
Phase 1, τs=1.0
6.42
6.42
Phase 2, τs=3.25
3.54
3.55
ImageNet-32
Phase 1, τs=1.0
8.33
8.33
Phase 2, τs=2.5
6.75
6.78
Table 1: Backprop-trained checkpoints sampled with feedforward autograd velocities ( vFFN ) versus PCN relaxation velocities ( vPCN ), at γ=10−4 (CIFAR-10) or 10−6 (ImageNet-32). Phases 1 and 2 are the Energy Matching training stages of Section 2 ; sampling seeds are matched within each row.
Method
Gradient
Params
FID ↓
ImprovedCD ( Du et al., 2021 )
backprop
—
25.1
SDiT (spiking) ( Yang et al., 2026 )
backprop
77M
11.34
CLEL-large ( Lee et al., 2022 )
backprop
31M
8.61
OT-CFM ∗ ( Tong et al., 2024 )
backprop
37M
4.04
DDPM++ ( Kim et al., 2022 )
backprop
62M
3.45
Energy Matching ( Balcerak et al., 2025 )
backprop
50M
3.34
Table 2: Unconditional CIFAR-10 FID, with our models after both training phases at each checkpoint’s optimal τs . The upper block is taken from the cited papers; ∗ marks reproductions by Balcerak et al. (2025) . ‡ Phase-2 negatives drawn with a feedforward twin (see text). Protocol and per-arm differences: App. J , App. I .
Appendix figures & tables13 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 3: Decomposition and budget independence of the velocity deviation, measured as in Figure 2 . (a) The relative velocity deviation ϵ=∥vPCN−vFFN∥/∥vFFN∥ resolved into components parallel (magnitude error) and perpendicular (direction error) to vFFN . Both scale as O(γ) , with direction the larger, ϵ⊥2 carrying roughly two-thirds of ϵ2 (the curves combine in quadrature). (b) At dt=1 the deviation is independent of the number of relaxation sweeps Kh at every γ tested: a single sweep already reaches the equilibrium.
Figure 4: The same sixteen noise realisations sampled by feedforward autograd velocities (top) and by PCN relaxation velocities (bottom), on our CIFAR-10 phase-2 replication checkpoint at γ=10−4 and τs=3.25 .
Figure 5: Relative velocity deviation ϵ on the authors’ phase-1 ImageNet-32 checkpoint, measured as in Figure 2 . The O(γ) law and both arithmetic floors reproduce their CIFAR-10 values, and the clamped potential, with the clamp incorporated into the top prediction function, follows the same law and shares the same floor. At the operating γ=10−6 (dotted) the deviation lies far below the floor of the sampling-run arithmetic.
Figure 6: cos(ΔθEP,ΔθBP) against the free-phase relaxation budget Tfree , at six values of γ , in TF32 (left) and float64 (right). CIFAR-10 replication checkpoint, 512 flow-matching interpolants, Tnudge=6 , Kh=1 , β=30 . Every γ that reaches full fidelity does so at the same budget; the plateau heights, by contrast, are specific to this β (Table 3 ).
β
10−5
3×10−5
10−4
10−3
10−2
10−1
0.3
0.0803
0.2353
0.7701
0.9933
0.9896
0.3391
1
0.3437
0.7149
0.9536
0.9988
0.9952
0.8840
3
0.6657
0.9459
0.9924
0.9996
0.9960
0.3935
10
0.9489
0.9920
0.9991
0.9996
0.9949
0.6788
30
0.9937
0.9989
0.9998
0.9990
0.9836
0.3902
100
0.9990
0.9998
0.9997
0.9969
0.9722
0.2807
Appendix
Table 3: Parameter-update fidelity cos(ΔθEP,ΔθBP) over nudge strength β (rows) and injection strength γ (columns), varied about the settings of our EP phase-1 training run (Table 6 ), in TF32 (the arithmetic of our training runs), on the CIFAR-10 replication checkpoint at 512 flow-matching interpolants. Shaded: cos<0.99 . Every failure in the grid belongs to the EP estimator, not the correspondence: the implicit solve, which uses no nudge, is exact to four decimal places up to γ=10−2 .
Figure 7: cos(ΔθEP,ΔθBP) against nudge strength at γ=10−4 , 512 interpolants, Tfree=14 , Tnudge=6 , Kh=1 . The quadratic nudge’s relaxation diverges above β=1 ; its values beyond that point are not true gradient estimates.
value
nCG (IFT)
Tnudge (EP)
Kh (EP)
1
1.0000
0.9994
0.9998
2
1.0000
0.9998
0.9998
3
1.0000
0.9998
0.9998
5
1.0000
—
0.9998
10
1.0000
0.9998
—
20
1.0000
—
—
Appendix
Table 4: Budget axes that are flat from a single sweep, at γ=10−4 on the CIFAR-10 replication checkpoint of Table 1 , 512 flow-matching interpolants, otherwise the settings of our EP training run ( Kh=1 , Tfree=14 , Tnudge=6 , β=30 ) with the named axis varied. Left: conjugate-gradient steps for the implicit solve, cos(ΔθIFT,ΔθBP) ; centre and right: nudged-phase sweeps and hidden-equilibration steps for EP, cos(ΔθEP,ΔθBP) . Dashes mark values not in that axis’s swept grid.
Phase-1 FID ↓
Gradient
50 k steps
100 k steps
145 k steps
Backpropagation
11.22
6.13
6.42
Implicit solve
12.35
6.45
6.70
EP
12.60
6.23
6.87
Appendix
Table 5: Phase-1 FID against training step, from 50k samples with EMA weights at τs=1.0 . All three arms improve substantially to 100 k and are slightly worse by 145 k under the constant learning rate, and the EP arm’s 100 k checkpoint is better than any arm’s 145 k, so its contrastive phase does not start from a weaker model. Repeated evaluations of one checkpoint span up to 0.15 , and we treat differences below about 0.2 as unresolved.
shared
Optimiser
Adam, β=(0.9,0.95) , gradient clipping at norm 1
Learning rate
1.2×10−3 , linear warmup over 10 k steps
Effective batch
1024 (phase 1), 512 (phase 2)
EMA decay
0.9999 (phase 1), 0.99 (phase 2)
Steps
145 k (phase 1; EP’s contrastive phase starts from its 100 k checkpoint) +2 k (phase 2)
Flow objective
u=x1−x0 , exact OT per batch, σ=0
Appendix
Table 6: Training settings. The upper block is shared by all arms; the lower gives the relaxation settings of the two PCN arms, with EP’s two phases separately (sweep budgets defined in App. B ).
Figure 8: CIFAR-10 samples from the backpropagation arm, FID 3.54 at τs=3.25 .
Figure 9: CIFAR-10 samples from the implicit-solve arm, FID 3.64 at τs=3.25 . Parameter gradients come from implicit differentiation of the PCN’s equilibrium rather than from backpropagation through a forward pass.
Figure 10: CIFAR-10 samples from the EP arm, FID 3.64 at τs=3.25 — a model trained without backpropagation in either phase.
Figure 11: CIFAR-10 samples from the norm- and attention-free architecture of Section 5.1 , FID 3.51 at τs=5.0 . Every operation in this network is a fixed weighted sum or a pointwise nonlinearity.