path: root/refs/scurria_2602.03670v2.txt

Equilibrium Propagation for Non-Conservative Systems

Antonino Emanuele Scurria 1 Dimitri Vanden Abeele 1 Bortolo Matteo Mognetti 2 Serge Massar 1

arXiv:2602.03670v2 [cs.LG] 1 Jun 2026

Abstract

from inference, the transmission of nonlocal error signals,
and synchronous layer-wise computations with explicit gradient storage. These constraints have no clear analog in
physical systems, making backpropagation challenging to
implement in neuromorphic or analog hardware. Consequently, understanding how credit assignment can instead
emerge from intrinsic system dynamics, through local interactions and continuous relaxation, is a central question in
neuroscience and machine learning.

Equilibrium Propagation (EP) is a physicsinspired learning algorithm that uses stationary
states of a dynamical system both for inference
and learning. In its original formulation it is
limited to conservative systems, i.e. to dynamics which derive from an energy function. Given
their applications, it is important to extend EP
to non-conservative systems, i.e. systems with
non-reciprocal interactions. Previous attempts to
generalize EP to such systems failed to compute
the exact gradient of the cost function. Here we
propose a framework that extends EP to arbitrary
non-conservative systems, including feedforward
networks. We keep the key property of equilibrium propagation, namely the use of stationary
states both for inference and learning. However,
we modify the dynamics in the learning phase by
a term proportional to the non-reciprocal part of
the interaction so as to obtain the exact gradient
of the cost function. This algorithm can also be
derived using a variational formulation that generates the learning dynamics through an energy
function defined over an augmented state space.
Numerical experiments show that this algorithm
achieves better performance and learns faster than
previous proposals.

Equilibrium Propagation (EP) (Scellier &Bengio, 2017)
represents one of the most promising advances in this direction. It formulates supervised learning as a contrast between
two stationary states of a dynamical system: a ‘free’ phase
where the system evolves autonomously, and a ‘nudged’
phase where outputs are weakly pushed toward their targets.
The local change in neural states between these phases recovers the exact gradient of the cost function with respect to
parameters. This enables spatially local learning exploiting
the continuous relaxation of the system without a distinct
backward circuit or explicit weight transport.
Since its introduction, several works have sought to improve the practicality and biological realism of EP. Algorithmic adaptations include enforcing temporal locality to
avoid state storage (Ernoult et al., 2020; Falk et al., 2025),
deriving agnostic updates for black-box energies (Scellier et al., 2022), and substituting nudging with clamping
(Stern et al., 2021). Theoretically, the framework has been
extended to stochastic systems (Scellier &Bengio, 2017;
Massar &Mognetti, 2025) and Lagrangian dynamics for
time-varying inputs (Massar, 2025; Pourcel et al., 2025;
Berneman &Hexner, 2025). In parallel, simulations have
explored suitable substrates, ranging from spiking (Martin et al., 2021; O’Connor et al., 2019) and resistive networks (Kendall et al., 2020) to coupled oscillators (Wang
et al., 2024; Rageau &Grollier, 2025), as well as quantum
systems (Wanjura &Marquardt, 2025; Massar &Mognetti,
2025; Scellier, 2024). Experimental realizations have been
demonstrated in memristor crossbars (Yi et al., 2023), selfadjusting electrical circuits (Dillavou et al., 2022; 2024),
elastic networks (Altman et al., 2024), and classical Ising
models trained on quantum annealers (Laydevant et al.,
2024).

1. Introduction
Standard neural network optimization relies on error backpropagation, an algorithm whose computational mechanism
is difficult to reconcile with biological (Crick, 1989) and
physical implementations (Indiveri &Liu, 2015). Specifically, backpropagation requires a backward pass distinct
1
Laboratoire d’Information Quantique (LIQ) CP224, Université
libre de Bruxelles (ULB), Av. F. D. Roosevelt 50, 1050 Bruxelles,
Belgium 2 Interdisciplinary Center for Nonlinear Phenomena and
Complex Systems CP231, Université libre de Bruxelles (ULB), Av.
F. D. Roosevelt 50, 1050 Bruxelles, Belgium. Correspondence to:
Antonino Emanuele Scurria <antonino.scurria@ulb.be>.

Proceedings of the 43 rd International Conference on Machine
Learning, Seoul, South Korea. PMLR 306, 2026. Copyright 2026
by the author(s).

Despite these recent developments and the theoretical el-

1

Equilibrium Propagation for Non-Conservative Systems

egance of EP, its standard formulation remains restricted
to conservative systems. In these systems, dynamics are
derived from an energy function, which inherently enforces
symmetry (e.g., symmetric synaptic connections Jij = Jji )
through the action-reaction principle. This constraint precludes the use of EP in a broad class of models characterized
by non-conservative forces. This includes the feedforward
architectures dominant in modern AI, biological circuits,
as well as physical systems that reach stationary states far
from thermodynamic equilibrium, such as nonlinear optical
systems driven by external lasers (Cin et al., 2025), optoelectronic systems (Kalinin et al., 2025), exciton-polariton
condensates (Sajnok &Matuszewski, 2025), active metamaterials (Brandenbourger et al., 2019) and active colloids
(Bishop et al., 2023; Osat &Golestanian, 2023) (see (Bowick
et al., 2022) for a review).

a framework where the original dynamics serve for inference, while a new augmented dynamic is used to compute
gradients of the cost Eq. (2). In this augmented phase, the
output neurons are nudged towards their targets (as in standard EP), while a local corrective term – proportional to the
antisymmetric part of the Jacobian at the free equilibrium
∂
JF (x0 , θ, u) = ∂x
F (x0 , θ, u) – is added to the forces. The
exact gradients of the cost with respect to parameters are
then obtained by contrasting stationary states of the augmented system.
Second, we introduce Dyadic EP, a ‘variational’ approach
to learning in non-conservative systems. This method involves doubling the number of variables in the system’s
state space and subsequently introducing a new energy function in this extended space. This approach takes advantage
of the extended space to execute the positive and negative
nudging phases in parallel, recovering the same computational cost as AsymEP. Derived from first principles, this
approach is inspired by established methods for mapping
dissipative dynamical systems onto conservative ones by
doubling the degrees of freedom (Bateman, 1931; Galley,
2013; Aykroyd et al., 2025). A more comprehensive study
of the theoretical framework and its application to feedforward networks can be found in (Scurria, 2026). Our method
is related to the Dual Propagation algorithm (Høier et al.,
2023; Høier &Zach, 2023; 2024) and constitutes an independent, first-principles generalization of Dyadic Learning
(Nest &Høier; Høier et al., 2024)—previously limited to
Hopfield networks—to arbitrary force fields.

Formally, we consider a dynamical system governed by
a non-reciprocal force field F (x, θ, u), which relaxes to a
stationary configuration x0 satisfying:
F (x0 , θ, u) = 0,

(1)

where x represents the state variables, θ the learnable parameters and u the static input. Our goal, given a target y(u), is
to compute the gradient of the cost function C(x0 , y) at this
equilibrium,
dC 0
(x , y),
(2)
dθ
and update θ to minimize the cost.
Previous attempts to extend EP to non-conservative dynamics include the Vector Field (VF) algorithm (Scellier et al.,
2018). However, as noted by the authors, this method provides an unbiased gradient of the cost Eq. (2) only in the
conservative case. To mitigate this, (Laborieux &Zenke,
2024) proposed adding a penalty to keep the Jacobian close
to symmetry, essentially forcing the system to be as conservative as possible. Alternative methods related to VF,
which similarly do not compute the exact gradient, were
proposed in (Farinha et al., 2020; Costa &Santos, 2025) and
for specific systems in simulation (Cin et al., 2025; Sajnok
&Matuszewski, 2025).

Third, we validate our framework on MNIST (LeCun, 1998),
Fashion-MNIST, and CIFAR-10. In continuous Hopfield
networks initialized with symmetric connection matrices,
AsymEP achieves better accuracy and learns faster than
EP and VF. Additionally, when we constrain the network
to have a strong degree of structural asymmetry, in which
case EP is inapplicable, AsymEP outperforms VF. Finally,
when we restrict connections to a feedforward structure, our
algorithm effectively trains all parameters; in contrast, VF
is limited to training the last layer, acting essentially as an
Extreme Learning Machine (Huang et al., 2006; Wang et al.,
2022) with poor performance.

Conversely, generalizations of backpropagation can handle
non-reciprocal forces and compute the exact gradient of
the cost Eq. (2) but inherit the same challenges in physical
implementations. For instance, Backpropagation Through
Time (Werbos, 1990) unfolds the network in time to apply standard backpropagation, Recurrent Backpropagation
(Almeida, 1990; Pineda, 1987) avoids this memory requirement but still requires a specific circuit to propagate errors,
and the continuous Adjoint Method (Chen et al., 2018) additionally requires integrating the dynamics backward in time
which is not physically possible for a dissipative system.

In summary, this theoretical work proposes two generalizations of EP beyond conservative systems to arbitrary differentiable dynamics that compute in their stationary states.

2. Equilibrium Propagation Overview
2.1. Conservative Systems
We first review standard Equilibrium Propagation (EP)
(Scellier &Bengio, 2017). We consider a network described
by an energy function E(x, θ, u), such that the force field is

In this paper, we first propose Asymmetric EP (AsymEP),
2

Equilibrium Propagation for Non-Conservative Systems

derived from the potential E:
FE (x, θ, u) = −

∂
E(x, θ, u).
∂x

stationary point, i.e., that Eq. (7) holds. Second, EP implic∂
itly assumes that the Jacobian JE (x0 , u) = ∂x
FE (x0 , u) is
invertible. In this work, we assume this condition holds and
will not state it explicitly hereafter. Third, for simplicity,
we omit the dependency on the input u and target y in the
following equations.

(3)

0
The objective is to compute the total gradient dC
dθ (x , y) of a
(quadratic) cost function C(x, y) evaluated at the minimum
energy configuration of the system. This free equilibrium
denoted x0 (which depend implicitly in θ and u), satisfies
the stationarity condition:

∂
− E(x0 , θ, u) = 0.
∂x

2.2. Vector Field
The Vector Field (VF) algorithm, introduced in (Scellier
et al., 2018), is an early attempt to adapt EP to nonreciprocal forces. This method relies on the observation
that, for conservative systems, linearizing the right-hand
side of Eq. (9) around the equilibrium point x0 yields

(4)

To compute gradients, we introduce the augmented energy
functional:
ET (x, θ, β, u, y) = E(x, θ, u) + βC(x, y),

1
β→0 2β

(5)

lim

where β is a scalar nudging parameter. The stationary configuration of this augmented system is obtained by integrating
the dynamics
∂ET (x, θ, β, u)
dx
=−
,
dt
∂x

where FE = −∂x E(x, θ) is the conservative force. It is
therefore tempting to use the right-hand side of Eq. (10) for
parameter updates of non-conservative systems, for which
no energy function E exists.

(6)

until the energy minimum is reached. This new fixed point
xβ , called nudged equilibrium, satisfies:
∂E(xβ , θ, u)
∂C(xβ , y)
+β
= 0.
∂x
∂x

The VF algorithm adopts precisely this approach. It uses
the nudged counterpart of Eq. (7),

(7)

F (xβ , θ) − β

The training procedure, as improved in (Laborieux et al.,
2021), uses two nudged phases with factors ±β (with
β ̸= 0). Starting from x0 , the system relaxes to two
nearby perturbed equilibria, x+β and x−β . The displacement x+β − x−β is then used to compute the parameter
update in the learning rule:


1 ∂E(xβ , θ, u) ∂E(x−β , θ, u)
−
, (8)
∆θ = −ϵ
2β
∂θ
∂θ

∂C β
(x ) = 0,
∂x

(11)

in conjunction with the learning rule Eq. (10):

∆θ = ϵ

∂F 0
(x , θ)
∂θ

⊤  β

x − x−β
.
2β

(12)

However, as noted in (Scellier et al., 2018), Eq. (12) does
0
not align with the true gradient dC
dθ (x ) and is exact only if
the force is conservative. To see this, let JF (x, θ) denote
the Jacobian of the vector field F (x, θ) (in components
i (x,θ)
(JF (x, θ))ij = ∂F∂x
). Differentiating the equilibrium
j
0
condition F (x , θ) = 0 with respect to θ gives

where ϵ > 0 is the learning rate. The theoretical foundation
of EP is the result that, in the limβ→0 of Eq. (8), we get:
dC(x0 , y)
d ∂E(xβ , θ, u)
=
,
dθ
dβ
∂θ


∂E(xβ , θ) ∂E(x−β , θ)
−
∂θ
∂θ

⊤  β
 (10)
x − x−β
∂FE 0
(x , θ)
= lim −
,
β→0
∂θ
2β



(9)

JF (x0 , θ)

see Appendix D.1. The error of the above method is O(β 2 ).
This error can be further reduced using holomorphic equilibrium propagation (Laborieux &Zenke, 2022).

dx0
∂F 0
+
(x , θ) = 0.
dθ
∂θ

(13)

Consequently, the exact gradient of the cost is

Thus, EP recovers the exact gradient of the cost function
using only local computations. In this manner, learning
implements gradient descent without an explicit backward
pass, and credit assignment is realized through the system’s
intrinsic relaxation dynamics.

⊤

dx0 ∂C 0
dC 0
(x ) =
(x )
dθ
dθ ∂x

⊤ 


−1 ∂C 0
∂F 0
⊤ 0
=−
(x , θ)
JF (x , θ)
(x ) .
∂θ
∂x
|
{z
}|
{z
}

Three remarks can be made at this point. First, EP does not
require the system to be at an energy minimum, but only at a

pre-synaptic

post-synaptic

(14)
3

Equilibrium Propagation for Non-Conservative Systems

Algorithm 1 Asymmetric EP (AsymEP)

The terms ’pre-synaptic’ and ’post-synaptic’ in Eq. (14)
are used by analogy with neuronal transmission: the presynaptic factor captures the local influence of θ on the force
F , while the post-synaptic factor is the sensitivity of the
cost to state perturbations.

1: Inputs: Force field F (x, θ), cost function C(x), nudg-

ing parameter β, learning rate ϵ.
2: repeat
3:
1. Free Phase: Evolve to stationary state
4:
Evolve the system dynamics
5:

If instead we differentiate the nudged equilibrium condition
in Eq. (11) with respect to β and evaluate at β = 0, we
obtain
β

JF (x0 , θ)

dx
dβ

−
β=0

∂C 0
(x ) = 0,
∂x

dx
= F (x, θ),
dt

(15)

6:
7:
8:
9:

which gives

−1 ∂C 0
dxβ
= JF (x0 , θ)
(x , y).
dβ β=0
∂x

(17)

until convergence to the stationary state x0 .
2. Jacobian Decomposition
Compute the Jacobian at equilibrium:
∂F 0
(x , θ),
(18)
∂x
and decompose it in its antisymmetric part:
JF (x0 , θ) =

(16)
10:
11:

The right-hand side of Eq. (16) represents the effective postsynaptic term used by the VF algorithm (Eq. 12). Comparing this with the exact post-synaptic term derived in Eq. (14),
we see that they coincide only if JF = JF⊤ , i.e., only if the
system is conservative.

AJ (x0 , θ) = 12 (JF (x0 , θ) − JF (x0 , θ)⊤ ). (19)
12:
13:
14:

Now, let SJ (x0 , θ) and AJ (x0 , θ) denote the symmetric
and antisymmetric parts of the Jacobian at the free (unnudged) equilibrium, respectively. Then, we show in Appendix A that the gradient error increases with the spectral

−1
radius of SJ (x0 , θ)
AJ (x0 , θ). Consequently, large
antisymmetric contributions degrade the gradient estimation, confirming empirical observations in the Appendix of
(Ernoult et al., 2020). In fact, in the pathological limit where
the Jacobian would be purely antisymmetric SJ (x0 , θ) = 0,
the update of VF gives the negative of the true gradient,
maximizing the cost rather than minimizing it.

15:
16:
17:
18:

3. Nudged Phase: Augmented Dynamics
Integrate the dynamics twice starting from x0
dx
∂C
= F (x, θ) − β
(x) − 2AJ (x0 , θ) (x − x0 ),
dt
∂x
(20)
until convergence to two new stationary states
x±β
A .
4. Parameter Update
Update the parameters according to:

∆θ = ϵ

3. Asymmetric EP

⊤
∂F 0
(x , θ)
∂θ

xβA − x−β
A
2β

!
.

(21)

19: until convergence of θ
20: Output: Optimized parameters θ.

Here, we introduce Asymmetric EP (AsymEP), see Algorithm 1, which removes the gradient estimate error inherent
to VF by adding a local correction term to the augmented
inference dynamics. The new nudged equilibrium xβA satisfies:
∂C β
F (xβA , θ) − β
(x ) − 2AJ (x0 , θ) (xβA − x0 ) = 0, (22)
∂x A

where JFA (x, θ) is the Jacobian of the modified dynamical
system Eq. (20). At the equilibrium point x0 , JFA is equal
to the transpose of the original Jacobian:
JFA (x0 , θ)

As in VF, we then obtain two perturbed states x±β
A for opposite nudging ±β and apply the contrastive learning rule
of Eq. (12).

=

JF (x0 , θ) − 2AJ (x0 , θ)

=

SJ (x0 , θ) − AJ (x0 , θ)

=

JF⊤ (x0 , θ).

(24)

We now show that AsymEP gives rise to the correct learning
rule, i.e. that right-hand side of Eq. (21) is proportional to
0
the gradient of the cost function dC
dθ (x ) at the equilibrium
0
point x (Eq. 14). To this end, note that the same reasoning
leading to Eq. (16) leads to

where we have used the decomposition Eq. (44) of the original Jacobian J into its symmetric and antisymmetric components. Therefore, the left hand side of Eq. (23) is equal to
the true post-synaptic term


−1 ∂C 0
dxβA
= JFA (x0 , θ)
(x ).
dβ β=0
∂x


−1 ∂C 0
dxβA
= JF⊤ (x0 , θ)
(x ),
dβ β=0
∂x

(23)
4

(25)

Equilibrium Propagation for Non-Conservative Systems

until a stationary point (z β , z ′β ) is reached. Upon convergence, we follow the standard EP paradigm in using the
difference z β − z ′β to compute the post-synaptic term. Un′
′
der the change of variables m = z+z
2 and d = z − z , we
prove in Appendix D that m follows the original dynamics
F (ensuring valid inference), while d relaxes to a "physical"
error signal proportional to the cost gradient.

which, using Eq. (14), proves the result. Additionally, although implied by the equality with the true gradient, we
explicitly demonstrate the equivalence of the gradient estimates obtained by AsymEP and Backpropagation Through
Time in Appendix B following (Ernoult et al., 2019).
Note that the corrective term −2AJ (x0 , θ)(x − x0 ) in
Eq. (20) is spatially local: AJ vanishes for unconnected
neurons, and (x − x0 ) is available at the synapse given the
memory mechanism already required by Eq. (12). This
correction can create backward connections (Section 5.3).
However, in physical realizations, both feedforward and
feedback connections must be physically present, though
feedback may be deactivated during inference.

It is important to notice that while Dyadic EP introduces a
distinct formulation, it remains consistent with the general
theoretical setting of EP and matches the computational
cost of AsymEP. Note also that we start the evolution of
the free phase (β = 0) with the identical initial condition
for z and z ′ , (i.e., d = 0). This guarantees that integrating Eq. (32) leads to a symmetric stationary point where
z 0 = z ′0 . Finally, we underline that the modified variational update rule in Eq. (34) is equivalent to the standard
symmetric EP update rule in Eq. (8) (see Appendix D).

4. Dyadic EP
We now introduce Dyadic EP (Algorithm 2), a variational
algorithm that computes the exact cost gradient in the limit
of infinitesimal nudging. It maps the original n-variable
dynamics F (x, θ) onto a 2n-variable system (z, z ′ ) defined
by an energy H(z, z ′ , θ) and cost D(z, z ′ ). We show in
Appendix E that AsymEP can be seen as the first-order
projection of Dyadic EP onto the original n-dimensional
state space.

Now, to make this concrete, consider a continuous Hopfield
network (see also Eq. (35)) with an asymmetric connection
matrix J. After some calculations (see Appendix F), the
augmented energy of the system can be re-expressed as:
1
1
HT = − ρ(z)⊤ Sρ(z) + ρ(z ′ )⊤ Sρ(z ′ ) − ρ(z)⊤ Aρ(z ′ )
2
2
1
β
2
′ 2
+ (∥z∥ − ∥z ∥ ) + (C(z, y) + C(z ′ , y)) ,
2
2
(29)
where S and A are the symmetric and antisymmetric parts
of J, respectively and ρ is an element-wise non-linearity.
An interesting analogy can be drawn with standard learning
rules in discrete Hopfield networks (Hopfield, 1982). For
a sequence of binary memories {ξ 1 , . . . , ξ m } where ξ µ ∈
corresponds to the standard autoassociative
{−1, 1}n , S P
Hebbian rule µ ξ µ (ξ µ )⊤ , creating stable attractors at each
pattern. In contrast, A corresponds to the heteroassociative
rule (e.g., a cycle between ξ µ and ξ ν given by ξ ν (ξ µ )⊤ −
ξ µ (ξ ν )⊤ ), encoding transitions between patterns.

The new system is defined by the energy H and cost function
D, given in terms of F and C by:


z + z′
H(z, z ′ , θ) = −(z − z ′ )⊤ F
,θ ,
2


′
z+z
,
(26)
D(z, z ′ ) = C
2
where z, z ′ ∈ Rn . In order to learn, we introduce the augmented energy
HT (z, z ′ , θ, β) = H(z, z ′ , θ) + βD(z, z ′ ).

(27)

The equilibrium configuration corresponds to a saddle point
of HT , where z minimizes and z ′ maximizes the energy.
This poses no issue for EP, which requires only that the
joint state (z, z ′ ) reaches a stationary state. Although this
min-maximization can be interpreted as z evolving forward
and z ′ backward in time, in practice they evolve forward
simultaneously, as we integrate the coupled equations:


dz
∂HT
z + z′
=−
=F
,θ
dt
2
∂z
⊤


z − z′
∂F
β ∂C z + z ′
+
−
,
2
∂z z+z′
2 ∂z
2
2


dz ′
∂HT
z + z′
=+
=F
,θ
′
dt
2
∂z



′ ⊤
z−z
∂F
β ∂C z + z ′
−
+
,
2
∂z ′ z+z′
2 ∂z ′
2
2
(28)

For this specific energy, the update rule given by Eq. (34)
can be re-expressed as:


⊤
1
∆J ∝ −
ρ(z ′β ) − ρ(z β ) ρ(z ′β ) + ρ(z β ) . (30)
2β
In the limit β → 0, this gives:
!
β
d
∆J ∝
⊙ ρ′ (m)ρ(m)⊤ .
β

(31)

matching the learning rule in (Pineda, 1987), with
β
limβ→0 dβ being the error signal.

5. Numerical Experiments
In this section, we numerically validate AsymEP (Algorithm 1). The neuronal dynamics follows the one introduced
5

Equilibrium Propagation for Non-Conservative Systems

where ∥ · ∥F denotes the Frobenius norm. Note that this
metric does not capture the asymmetry of the Jacobian,
which depends on the state x.

Algorithm 2 Dyadic EP
1: Inputs: Force field F (x, θ), cost function C(x, y),

nudging parameter β, learning rate ϵ
2: repeat
3:
1. Free Phase: Evolve to stationary state
4:
Evolve the system dynamics, starting from identi-

For numerical experiments, we restricted the network to a
layered architecture with a single hidden layer to facilitate
comparison with prior work. Accordingly, J in contains
only input-to-hidden connections, while J dyn is block offdiagonal, encoding bidirectional interactions between the
hidden and output layers. Both J in and J dyn are trained.

cal initial conditions z(0) = z ′ (0) = z0 ,
5:

6:
7:
8:

dz
∂H
=−
,
dt
∂z

11:
12:
13:

(32)

We first use MNIST (LeCun, 1998) (60k train, 10k test)
followed by Fashion-MNIST to validate AsymEP, and then
we further validate AsymEP and Dyadic EP by comparing
them to Backpropagation on a convolutional feedforward,
with CIFAR-10. Inputs are normalized using min-max to
[−1, 1] and targets are one-hot encoded in {−1, 1}. All
hyperparameters are detailed in Appendix G, along with
additional details and numerical results.

until stationary states z 0 , z ′0 are reached.
2. Nudged Equilibrium
Evolve the system dynamics, starting from the
solution of the free phase z 0 = z ′0 :

9:

10:

dz ′
∂H
=+ ′,
dt
∂z

dz
∂HT
=−
,
dt
∂z

dz ′
∂HT
,
=+
dt
∂z ′

(33)

until two nudged stationary states z β , z ′β are
reached.
3. Parameter Update
Update the parameters according to:


1 ∂H(z β , z ′β , θ)
∆θ = −ϵ
(34)
β
∂θ

5.1. Symmetric Initialization
We start by comparing AsymEP with standard EP and
VF. All algorithms are initialized with an identical symmetric matrix J dyn . EP maintains this symmetry throughout training, while VF and AsymEP induce asymmetry in
J dyn . Since EP and VF already achieve strong performance
on MNIST, the purpose of this experiment is to validate
AsymEP and compare it against EP and VF rather than
outperform the state of the art.

14: until convergence of θ
15: Output: Optimized parameters θ.

in (Scellier &Bengio, 2017), and is generalized to allow
for non-reciprocal forces as in (Scellier et al., 2018). For
clarity, we express the forces in a form that explicitly separates the contributions of the external input and the recurrent
interactions:


F (x) = ρ′ (x) ⊙ J in u + J dyn ρ(x) − x,
(35)

Figure 1 compares the three algorithms as a function of
hidden-layer dimension after 1 and 20 epochs. AsymEP
consistently outperforms the baselines, suggesting it learns
faster and better.
Figure 2 studies the evolution of the asymmetry ratio rstr .
The results are reported for 50 hidden neurons. As expected,
EP preserves the initial weight symmetry. In contrast, VF
and AsymEP induce non-trivial evolution of rstr following
two distinct patterns, resulting in three distinct network
configurations. A complementary figure is available in Appendix G.1.

where u ∈ RNin denotes the input and x ∈ RNdyn the neuronal state, comprising both hidden and output units. The
matrices J in ∈ RNdyn ×Nin and J dyn ∈ RNdyn ×Ndyn define the
input and recurrent connectivity, respectively. The activation
function ρ(·) is taken to be the hyperbolic tangent, applied
element-wise.

5.2. Fixed Asymmetry Ratio

If J dyn is symmetric, we can define the energy:
E(x) =

While the previous section focused on networks compatible
with all three algorithms (EP, VF, AsymEP), we now turn
to architectures with strong structural asymmetry. In this
regime, EP is inapplicable by construction, and, as we show,
VF performs poorly, contrary to AsymEP which remains
effective.

1
1
∥x∥2 − ρ(x)⊤ J dyn ρ(x) − ρ(x)⊤ J in u, (36)
2
2

which is identical to that of (Scellier &Bengio, 2017), provided that the input neurons are activated as ρ(u).
Equation (35) naturally motivates a quantitative measure of
structural asymmetry rstr , defined as:

To this end, we consider a class of networks where the
asymmetry ratio rstr defined in Eq. (37) is kept fixed. Let S̃
and à be arbitrary symmetric and antisymmetric matrices
in RNdyn ×Ndyn respectively. We enforce a fixed rstr via the

⊤

rstr =

∥(J dyn − J dyn )/2∥F
,
∥J dyn ∥F

(37)
6

Equilibrium Propagation for Non-Conservative Systems

where γ ∈ R is a learnable global scale.
Using VF and AsymEP, we train a layered network with one
hidden layer of 50 neurons (in which case S̃ and à are block
off-diagonal) for different values of rstr to investigate the
impact of structural asymmetry. We compare two training
regimes: training only the input weights J in (and the scale
γ), versus training all parameters including J dyn . The first
regime trains only the external forces from the input ρ′ (x) ⊙
J in u (which correspond to a symmetric contribution in the
Jacobian) applied to our non-conservative system, while
the second additionally trains J dyn and therefore the nonsymmetric part of the Jacobian directly.

(a) Results after one epoch.

Figure 3 summarizes the results. We find that AsymEP
maintains robust performance across all asymmetry levels
(e.g., achieving an accuracy of 93.8 ± 0.4% at rstr = 0 and
94.9 ± 0.2% at rstr = 0.875 when training all parameters)
and can even learn when the recurrent connection matrix
J dyn is completely antisymmetric (rstr = 1). Additionally,
training all parameters shows significant improvement over
training only J in .
In contrast, VF performs well at low asymmetry ratios
but degrades as asymmetry increases, eventually dropping
to chance levels (e.g., accuracies of 5 ± 3% and 8 ± 4%
at rstr = 1 for input-only and all-parameter training, respectively). When only J in is trained, VF accuracy collapses around rstr ≈ 0.5, whereas training all parameters
delays this collapse until rstr ≈ 0.8. Our analysis in Appendix G.2.1 reveals that VF adjusts the dynamics such that
the asymmetry of the Jacobian’s off-diagonal terms remains
strictly lower than the structural asymmetry ratio. The training appears to adjust the neuronal state such that neurons
connected by strongly asymmetric weights have low activation. As shown in Appendix G.2.1, AsymEP learns faster
than VF across all levels of asymmetry.

(b) Results after 20 epochs.
Figure 1. Comparison of algorithm performance on MNIST using
a layered architecture with one hidden layer and symmetric initialization. Squares denote AsymEP, circles EP, and triangles VF.
Test accuracy (averaged over 10 runs) is shown after one epoch
(Fig. 1a) and 20 epochs (Fig. 1b).

Finally, Appendix G.3 opens with a brief theoretical discussion of the stability of these non-conservative dynamics,
followed by simulations on all-to-all topologies with constrained rstr and input projections J in . Even in this worstcase setting, AsymEP reduces oscillations and improves
stability.
5.3. Feedforward Architectures
Figure 2. Evolution of the asymmetry ratio rstr (defined in Eq. (37))
during training on MNIST for AsymEP, EP and VF, initialized
from a symmetric configuration. The models use 50 hidden neurons.

We now consider a purely feedforward architecture. Here
VF trains only the last layer: with no backward connections,
the output nudging signal cannot reach earlier layers, so for
every layer but the last the nudged stationary states coincide
with the free states, giving zero weight updates. As only
the output layer is trained, the system essentially becomes
an Extreme Learning Machine (Huang et al., 2006; Wang
et al., 2022). In contrast, AsymEP introduces a correction
that generates effective backward connections, allowing the

following parameterization of the recurrent parameters:
"q
#
S̃
Ã
dyn
2
J =γ
1 − rstr
+ rstr
,
(38)
∥S̃∥F
∥Ã∥F
7

Equilibrium Propagation for Non-Conservative Systems

tivity structures inspired by (Millidge et al., 2023), while
keeping the number of trainable parameters fixed.
Experiments are conducted on Fashion-MNIST using a twohidden-layer network with hidden dimensions 500 and 200.
Network states are denoted (x0 , x1 , x2 , x3 ), where x0 is
the input and x3 = xL the output. Forward and backward
connections are denoted by Wk and Bk , respectively, with
W1 = J in .
We consider three classes of dynamics. First, the Continuous
Hopfield (CH) dynamics introduced previously:


dxk
= −xk +ρ′ (xk )⊙ Wk ρ(xk−1 )+(1−δk,L )Bk ρ(xk+1 ) .
dt
(40)
Second, Predictive Coding (PC) dynamics, defined through
the prediction errors ek = xk − Wk ρ(xk−1 ), whose fixed
point ek = 0 corresponds to a standard feedforward network:

Figure 3. Impact of the structural asymmetry ratio rstr on accuracy
(top) and standard deviation over 10 runs (bottom) on MNIST.
We compare VF (orange) and AsymEP (blue) under two training
regimes: training only J in (dashed) or all parameters (solid).

dxk
= −ek + (1 − δk,L ) (ρ′ (xk ) ⊙ (Bk ek+1 )) .
dt

nudging signal to influence all layers. We make this explicit
for a network with one hidden layer.

(41)

Third, a standard dynamics chosen for direct comparison
with backpropagation:

Let the state x be partitioned in hidden h and output o
layers. The recurrent connection matrix is then J dyn =

0
0
. The forces of the system are:
Wh→o 0
 β



Fh = ρ′ (h) ⊙ J in u + λ(Wh→o )⊤ (o − o0 ) − h





 β
0
Fo = ρ′ (o) ⊙ Wh→o ρ(h) − λWh→o (h − h )
(39)



∂C


−o
+ λβ
∂o

dxk
= −xk + Wk ρ(xk−1 ) + (1 − δk,L )Bk ρ(xk+1 ). (42)
dt
For each dynamics, we examine three connectivity scenarios.
⊤
• In the asymmetric case (Bk ̸= Wk+1
), the backward
weights Bk are randomly initialized and kept fixed
while only the forward weights are trained, ensuring a
fair comparison (i.e., identical number of parameters);
in PC, the learning rule for Bk is zero when only inputs
are clamped.

where λ is 0 during the free inference and 1 during the
nudged phase (Eq. 20). The force on the hidden layer Fhβ
now depends on the output layer through the term ρ′ (h) ⊙
⊤
(Wh→o ) (o − o0 ), enabling the nudge (the term β ∂C
∂o ) to
influence the hidden layer. This implicitly assumes that the
hardware implementation supports the physical activation
of these backward connections.

⊤
• In the symmetric / conservative case (Bk = Wk+1
), the
CH and PC dynamics derive from an energy functional,
while the standard dynamics remains non-conservative
due to its non-symmetric Jacobian.

• In the feedforward case (Bk = 0), the PC and standard dynamics coincide; for the standard dynamics, the
AsymEP learning rule mirrors backpropagation, with
1
∆xβk = 2β
(xβ − x−β ) acting as the propagated error
signal.

We validate this using a single hidden layer of only 20 neurons on MNIST. After training, VF saturates with 64.3 ±
2.0% accuracy, whereas AsymEP reaches 92.7 ± 0.5% accuracy. We expect this discrepancy to increase with network
depth, since this increases the number of layers unable to
learn under VF. A figure with the accuracy during training
can be found in Appendix G.4.2.

Table 1 shows that AsymEP consistently outperforms VF
in both asymmetric and feedforward settings, in final accuracy, learning speed, and stability. After a single epoch
it already provides on average a 15% accuracy gain with
an order-of-magnitude reduction in variance. Remarkably,
AsymEP with asymmetric connectivity also surpasses EP
on symmetric networks despite training only the forward

5.4. Advantages of Non-Conservative Dynamics
AsymEP is not tied to a specific neural dynamics. To further
assess the benefits of training non-conservative dynamics
using AsymEP, we compare several dynamics and connec8

Equilibrium Propagation for Non-Conservative Systems

6. Discussion and Conclusion

weights, suggesting that relaxing symmetry constraints may
improve expressivity. Supplementary results are provided
in Appendix G.5.

In this work, we extended Equilibrium Propagation (EP)
to non-conservative systems that reach stationary states by
deriving two mathematically equivalent algorithms that recover the exact gradient of the cost function in the limit of
infinitesimal nudging.

Table 1. Test accuracy on Fashion-MNIST (%) at Epoch 50 (mean
± std 10 runs). BP on a standard feedforward architecture using
MSE and SGD achieve 87.37 ± 0.29%.
EP

AsymEP

Asym
86.78 ± 0.14
Feedfor
86.05 ± 0.12
Sym
84.30 ± 0.13
Asym
86.20 ± 0.17
PC
Sym
84.78 ± 0.14
Asym
82.91 ± 0.48
Standard
Feedfor
86.25 ± 0.16
CH

The first approach, Asymmetric EP, preserves the original
inference dynamics. It introduces a corrective force during
the nudged phase that remains spatially local, as the antisymmetric Jacobian is null for unconnected neurons and the
perturbation from equilibrium is available at the synapse
level. Unlike standard methods like Recurrent Backpropagation (Almeida, 1990; Pineda, 1987), this avoids explicit
digital weight transposition. However, a physical mechanism to obtain the local corrective force at the synapse
level remains a subject for future work. We also note that
AsymEP shares the temporal non-locality of standard EP.

VF
85.20 ± 0.12
77.76 ± 0.37
80.71 ± 6.17
75.52 ± 1.69
78.58 ± 0.28

Finally, to investigate how AsymEP scales with depth, we
trained deeper fully connected networks with two and three
hidden layers of 500 neurons on Fashion-MNIST, reaching
86.41 ± 0.22% and 87.8 ± 0.15% test accuracy respectively.

The second approach, Dyadic EP, doubles the state space
to map non-reciprocal dynamics onto an energy landscape—conceptually reminiscent of multi-compartment cortical neurons, where apical dendrites integrate feedback
(analogous to z − z ′ ) separately from basal feedforward
input (analogous to z + z ′ ) (Guerguiev et al., 2017). Additionally, this expanded space also enables the positive and
negative nudging phases to run in parallel. This offers a
pathway to implement a version of EP that is local in time,
but would require a doubling of the degrees of freedom
on the physical hardware. More fundamentally, the energy
defined on the extended state shows that the tools and theoretical guarantees obtained for EP should also apply to
the case of non-reciprocal forces, and that the variational
principle behind EP is universal in the sense that it can be
applied to all networks which operate in a stationary state.

5.5. Feedforward Training on CIFAR-10: BP vs. Dyadic
EP vs. AsymEP
To test whether our framework scales beyond shallow networks, we conclude with a deep, purely feedforward CNN
architecture trained on CIFAR-10. We compare backpropagation (BP), VF, AsymEP and Dyadic EP in a controlled
setting where the gradient estimator is the only difference
between runs: all methods share the same configuration,
with the BP gradient replaced by the contrast of stationary
states for the EP-based methods (see App. G.6 for details).
Each configuration is trained for 40 epochs over 5 seeds.

Furthermore, Dyadic EP is not restricted to the EP community and could suggest a more physically plausible alternative to the stationary-state Adjoint Method (for fixed
inputs) (Chen et al., 2018): by solving the forward and adjoint equations simultaneously via relaxation, it circumvents
a separate backward-in-time pass.

Table 2 reports the final test accuracy. Both of our algorithms scale to this regime, closely tracking the BP baseline
throughout training and matching its final accuracy: a paired
t-test finds no significant difference between Dyadic EP and
BP (p = 0.75), and only a sub-percent gap for AsymEP.
In contrast, VF makes slight initial progress (peaking near
30%) before collapsing to chance level (10%). Additional
details can be found in Appendix G.6

Finally, our experiments on MNIST, Fashion-MNIST, and
CIFAR-10 confirm that AsymEP and Dyadic EP consistently outperform EP and VF, and notably enables effective
training of feedforward networks.
Our work thus opens new avenues for learning in neuromorphic hardware, dissipative physical systems, and neural
architectures where asymmetry is intrinsic rather than incidental.

Table 2. Test accuracy on CIFAR-10 (%) at epoch 40 (mean ± std
over 5 seeds).
Method

Test Acc. (%)

Backpropagation
Dyadic EP
AsymEP
VF

90.66 ± 0.25
90.69 ± 0.14
89.74 ± 0.14
10.00 ± 0.00

9

Equilibrium Propagation for Non-Conservative Systems

Impact Statement

References

This paper presents results that advance the field of machine
learning. There are many potential societal consequences
of our work, none of which we feel must be specifically
highlighted here.

Almeida, L. B. A learning rule for asynchronous perceptrons with feedback in a combinatorial environment. In
Artificial neural networks: concept learning, pp. 102–111.
1990.

Acknowledgments

Altman, L. E., Stern, M., Liu, A. J., and Durian, D. J. Experimental demonstration of coupled learning in elastic
networks. Physical Review Applied, 22(2):024053, 2024.

AES is fully funded by the Horizon Europe Marie
Skłodowska-Curie Doctoral Network ’Postdigital Plus’
(Grant 101169118). DVA acknowledges the support of
the French Community of Belgium through a FRIA fellowship. SM acknowledges financial support by the Fonds de la
Recherche Scientifique–FNRS, Belgium under EOS Project
No. 40007536. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif
(CÉCI), funded by the Fonds de la Recherche Scientifique
de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and
by the Walloon Region.

Aykroyd, C., Bourgoin, A., and Poncin-Lafitte, C. L. Hamiltonian treatment of non-conservative systems. arXiv
preprint arXiv:2507.18658, 2025.
Bateman, H. On dissipative systems and related variational
principles. Physical Review, 38(4):815, 1931.
Berneman, M. and Hexner, D. Equilibrium propagation for
dissipative dynamics. Advanced Intelligent Systems, pp.
e202501310, 2025.
Bishop, K. J., Biswal, S. L., and Bharti, B. Active colloids
as models, materials, and machines. Annual Review of
Chemical and Biomolecular Engineering, 14(1):1–30,
2023.

“ἁρμονίη ἀφανὴς φανερῆς κρείττων”

Bowick, M. J., Fakhri, N., Marchetti, M. C., and Ramaswamy, S. Symmetry, thermodynamics, and topology
in active matter. Physical Review X, 12(1):010501, 2022.
Brandenbourger, M., Locsin, X., Lerner, E., and Coulais, C.
Non-reciprocal robotic metamaterials. Nature communications, 10(1):4608, 2019.
Cesa-Bianchi, N. and Lugosi, G. Prediction, learning, and
games. Cambridge university press, 2006.
Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud,
D. K. Neural ordinary differential equations. Advances
in neural information processing systems, 31, 2018.
Cin, N. D., Marquardt, F., and Wanjura, C. C. Training
nonlinear optical neural networks with scattering backpropagation. arXiv preprint arXiv:2508.11750, 2025.
Costa, P. and Santos, P. A. Directed equilibrium propagation
revisited. Mathematics, 13(11), 2025. ISSN 2227-7390.
Crick, F. The recent excitement about neural networks.
Nature, 337, 1989.
Dillavou, S., Stern, M., Liu, A. J., and Durian, D. J. Demonstration of decentralized physics-driven learning. Physical Review Applied, 18(1):014040, 2022.
Dillavou, S., Beyer, B. D., Stern, M., Liu, A. J., Miskin,
M. Z., and Durian, D. J. Machine learning without a processor: Emergent learning in a nonlinear analog network.
Proceedings of the National Academy of Sciences, 121
(28):e2319718121, 2024.
10

Equilibrium Propagation for Non-Conservative Systems

Ernoult, M., Grollier, J., Querlioz, D., Bengio, Y., and Scellier, B. Updates of equilibrium prop match gradients
of backprop through time in an rnn with static input.
Advances in neural information processing systems, 32,
2019.

Indiveri, G. and Liu, S.-C. Memory and information processing in neuromorphic systems. Proceedings of the
IEEE, 103(8):1379–1397, 2015.

Ernoult, M., Grollier, J., Querlioz, D., Bengio, Y., and Scellier, B. Equilibrium propagation with continual weight
updates. arXiv preprint arXiv:2005.04168, 2020.

Kalinin, K. P., Gladrow, J., Chu, J., Clegg, J. H., Cletheroe,
D., Kelly, D. J., Rahmani, B., Brennan, G., Canakci, B.,
Falck, F., et al. Analog optical computer for ai inference
and combinatorial optimization. Nature, 645(8080):354–
361, 2025.

Falk, M. J., Strupp, A. T., Scellier, B., and Murugan,
A. Temporal contrastive learning through implicit nonequilibrium memory. Nature Communications, (16),
2025.

Kendall, J., Pantone, R., Manickavasagam, K., Bengio,
Y., and Scellier, B. Training end-to-end analog neural
networks with equilibrium propagation. arXiv preprint
arXiv:2006.01981, 2020.

Farinha, M. T., Pequito, S., Santos, P. A., and Figueiredo,
M. A. T. Equilibrium propagation for complete directed
neural networks. In Proceedings of the 28th European
Symposium on Artificial Neural Networks, Computational
Intelligence and Machine Learning (ESANN 2020), 2020.

Laborieux, A. and Zenke, F. Holomorphic equilibrium
propagation computes exact gradients through finite size
oscillations. Advances in Neural Information Processing
Systems, 35:12950–12963, 2022.
Laborieux, A. and Zenke, F. Improving equilibrium propagation without weight symmetry through jacobian homeostasis. In Proceedings of the International Conference on Learning Representations (ICLR) 2024, Virtual
(ICLR), May 2024.

Galley, C. R. Classical mechanics of nonconservative systems. Physical review letters, 110(17):174301, 2013.
Guerguiev, J., Lillicrap, T. P., and Richards, B. A. Towards
deep learning with segregated dendrites. elife, 6:e22901,
2017.

Laborieux, A., Ernoult, M., Scellier, B., Bengio, Y., Grollier,
J., and Querlioz, D. Scaling equilibrium propagation to
deep convnets by drastically reducing its gradient estimator bias. Frontiers in neuroscience, 15:633674, 2021.

Høier, R. and Zach, C. A lagrangian perspective on dual
propagation. In Proceedings of the First Workshop on Machine Learning with New Compute Paradigms at NeurIPS
2023, New Orleans, LA, USA, Dec 2023.

Laydevant, J., Marković, D., and Grollier, J. Training an
ising machine with equilibrium propagation. Nature Communications, 15(1):3671, 2024.

Høier, R. and Zach, C. Two tales of single-phase contrastive
hebbian learning. In Salakhutdinov, R., Kolter, Z., Heller,
K., Weller, A., Oliver, N., Scarlett, J., and Berkenkamp, F.
(eds.), Proceedings of the 41st International Conference
on Machine Learning, volume 235 of Proceedings of
Machine Learning Research, pp. 18470–18488. PMLR,
21–27 Jul 2024.

LeCun, Y. The mnist database of handwritten digits.
http://yann. lecun. com/exdb/mnist/, 1998.
Martin, E., Ernoult, M., Laydevant, J., Li, S., Querlioz, D.,
Petrisor, T., and Grollier, J. Eqspike: spike-driven equilibrium propagation for neuromorphic implementations.
Iscience, 24(3), 2021.

Høier, R., Staudt, D., and Zach, C. Dual propagation: accelerating contrastive hebbian learning with dyadic neurons.
In Proceedings of the 40th International Conference on
Machine Learning, ICML’23. JMLR.org, 2023.

Massar, S. Equilibrium propagation for learning in lagrangian dynamical systems. Physical Review E, 112
(3):035304, 2025.

Høier, R., Kalinin, K., Ernoult, M., and Zach, C. Dyadic
learning in recurrent and feedforward models. In NeurIPS
2024 Workshop Machine Learning with new Compute
Paradigms, 2024.

Massar, S. and Mognetti, B. M. Equilibrium propagation:
the quantum and the thermal cases. Quantum Studies:
Mathematics and Foundations, 12(1):6, 2025.

Hopfield, J. J. Neural networks and physical systems with
emergent collective computational abilities. Proceedings
of the national academy of sciences, 79(8):2554–2558,
1982.

Millidge, B., Song, Y., Salvatori, T., Lukasiewicz, T., and
Bogacz, R. Backpropagation at the infinitesimal inference limit of energy-based models: Unifying predictive
coding, equilibrium propagation, and contrastive hebbian
learning. In International Conference on Learning Representations (ICLR), 2023.

Huang, G.-B., Zhu, Q.-Y., and Siew, C.-K. Extreme learning
machine: theory and applications. Neurocomputing, 70
(1-3):489–501, 2006.
11

Equilibrium Propagation for Non-Conservative Systems

Nest, T. and Høier, R. Dyadic learning in asymmetric
convnets. In New Frontiers in Associative MemoriesWorkshop at ICLR 2026.

Wang, Q., Wanjura, C. C., and Marquardt, F. Training
coupled phase oscillators as a neuromorphic platform
using equilibrium propagation. Neuromorphic Computing
and Engineering, 4(3):034014, 2024.

Osat, S. and Golestanian, R. Non-reciprocal multifarious
self-organization. Nature Nanotechnology, 18(1):79–85,
2023.

Wanjura, C. C. and Marquardt, F. Quantum equilibrium
propagation for efficient training of quantum systems
based on onsager reciprocity. Nature Communications,
16(1):6595, 2025.

O’Connor, P., Gavves, E., and Welling, M. Training a spiking neural network with equilibrium propagation. In The
22nd international conference on artificial intelligence
and statistics, pp. 1516–1523. PMLR, 2019.

Werbos, P. J. Backpropagation through time: what it does
and how to do it. Proceedings of the IEEE, 78(10):1550–
1560, 1990.

Pineda, F. Generalization of back propagation to recurrent
and higher order neural networks. In Neural information
processing systems, 1987.

Yi, S.-i., Kendall, J. D., Williams, R. S., and Kumar, S.
Activity-difference training of deep neural networks using
memristor crossbars. Nature Electronics, 6(1):45–51,
2023.

Pourcel, G., Basu, D., Ernoult, M., and Gilra, A. Lagrangianbased equilibrium propagation: generalisation to arbitrary boundary conditions & equivalence with hamiltonian echo learning. arXiv preprint arXiv:2506.06248,
2025.
Rageau, T. and Grollier, J. Training and synchronizing
oscillator networks with equilibrium propagation. Neuromorphic Computing and Engineering, 2025.
Sajnok, K. and Matuszewski, M. Near-equilibrium propagation training in nonlinear wave systems. arXiv preprint
arXiv:2510.16084, 2025.
Scellier, B. Quantum equilibrium propagation: Gradientdescent training of quantum systems. arXiv preprint
arXiv:2406.00879, 2024.
Scellier, B. and Bengio, Y. Equilibrium propagation: Bridging the gap between energy-based models and backpropagation. Frontiers in computational neuroscience, 11:24,
2017.
Scellier, B., Goyal, A., Binas, J., Mesnard, T., and Bengio,
Y. Generalization of equilibrium propagation to vector
field dynamics. arXiv preprint arXiv:1808.04873, 2018.
Scellier, B., Mishra, S., Bengio, Y., and Ollivier, Y. Agnostic
physics-driven deep learning. arXiv:2205.15021v1, 2022.
Scurria, A. E. A physical theory of backpropagation: Exact
gradients from the least-action principle. 2026.
Stern, M., Hexner, D., Rocks, J. W., and Liu, A. J. Supervised learning in physical networks: From machine
learning to learning machines. Physical Review X, 11(2):
021045, 2021.
Wang, J., Lu, S., Wang, S.-H., and Zhang, Y.-D. A review
on extreme learning machine. Multimedia Tools and
Applications, 81(29):41611–41660, 2022.
12

Equilibrium Propagation for Non-Conservative Systems

A. Gradient Estimation Error in VF

where s denotes the dynamical state of the system. This
symmetry is the linchpin of the equivalence proof, as the
gradient expressions derived for BPTT and standard EP
differ precisely by a transpose operation applied to ∂F
∂s .

In this appendix, we quantify the gradient estimation error
introduced by VF in the limit where the Jacobian asymmetry
is small.

This observation aligns with our analysis in the main text:
VF fails in non-conservative systems due to the missing
transpose in the post-synaptic term (see Eq. (16)). Following
the derivation in Ernoult et al. (2019) (viz., Appendix A, Eqs.
(31–33)), the recursive relations for the gradients in BPTT
are given by:

Comparing the post-synaptic update terms in Eqs. (12) and
(14) gives the following error in the gradient of the cost:
⊤
∂F 0
(x , θ)
∂θ


−1

−1  ∂C 0
(x , y), (43)
× JF (x0 , θ)
− JF⊤ (x0 , θ)
∂x


Error = −

∇BPTT
(0) =
s

To quantify this error, we decompose the Jacobian JF (x, θ)
into its symmetric part SJ (x, θ) and antisymmetric part


SJ (x, θ) = 12 JF (x, θ) + JF⊤ (x, θ) ,
(44)


AJ (x, θ) = 12 JF (x, θ) − JF⊤ (x, θ) .

(JF )

=

(JF⊤ )−1 =

∞
X

∇BPTT
(t) =
s
∇BPTT
(t) =
θ

(−1)

n=0
∞
X

(SJ−1 AJ )n

SJ−1 ,




∂F
(x, s⋆ , θ)
∂s

⊤

∂F
(x, s⋆ , θ)
∂θ

⊤

∇BPTT
(t − 1),
s

(50)

∇BPTT
(t − 1),
s

(51)

where θ represents the optimization parameters, ℓ is the
cost function, s⋆ is the free equilibrium state (satisfying
F (s⋆ ) = 0), y is the target, and x is the input. The index t
denotes the unrolled time steps, initialized at s(0) = s⋆ .

!
n

(49)

and for all t = 1, . . . , K,

Assuming the asymmetry AJ (x, θ) is small, we can make
a series expansion in SJ−1 AJ (omitting the dependencies
for clarity). Applying the Neumann expansion for small
∥SJ−1 AJ ∥ gives
−1

∂ℓ
(s⋆ , y),
∂s

(45)

In contrast, the gradients computed by VF follow the recursion (viz., Ernoult et al. (2019), Appendix A, Eqs. (24–26)):

!
(SJ−1 AJ )n

SJ−1 .

(46)

∆EP
s (0) = −

n=0

∂ℓ
(s⋆ , y),
∂s

(52)

and for all t ≥ 0,

Subtracting the two series and assuming convergence, we
finally obtain
!
∞
X

2n+1
−1
−1
⊤ −1
(JF ) − (JF ) = −2
SJ AJ
SJ−1 .

∂F
(x, s⋆ , θ) ∆EP
s (t),
∂s
⊤

∂F
(x,
s
,
θ)
∆EP
∆EP
(t
+
1)
=
⋆
s (t).
θ
∂θ

∆EP
s (t + 1) =

n=0

(47)

(53)
(54)

Comparing these two sets of equations confirms that the only
difference are Eqs. (50) and (53), specifically the transpose
of the Jacobian ∂F
∂s (ignoring the global sign difference in
Eqs. (49) and (52)).

B. Equivalence between AsymEP and BPTT
In this appendix, we sketch the equivalence between the
gradient estimate computed by AsymEP and Backpropagation Through Time (BPTT) (Werbos, 1990) for a Recurrent
Neural Network with fixed inputs. Our derivation relies on
the proof provided by Ernoult et al. (2019), which established that standard (conservative) EP computes gradients
identical to those of BPTT. To facilitate direct comparison,
we adopt their notation for this section.

In AsymEP, we modify the dynamics by adding a correction
term dependent on the antisymmetric part of the Jacobian.
Denoting the force of this augmented system by F A , the
Jacobian at the free equilibrium satisfies:

The proof provided by Ernoult et al. (2019) relies on the
assumption that the vector field F (i.e., transition function)
is derived from a scalar potential function, which implies
that

⊤
∂F
∂F
=
,
(48)
∂s
∂s

By substituting this corrected Jacobian into the recursive
relations, AsymEP recovers the exact transpose required
by BPTT. Consequently, our method extends the equivalence between EP and BPTT to the general case of nonconservative force.

∂F A
(x, s⋆ , θ) =
∂s

13



∂F
(x, s⋆ , θ)
∂s

⊤
.

(55)

Equilibrium Propagation for Non-Conservative Systems

C. Out-of-Equilibrium Mechanics

C.3. Symmetry Breaking as Credit Assignment

Here we sketch the physical picture behind the doubledenergy construction of Eq. (26). The full derivation from
Hamilton’s least-action principle, together with its connection to the Bateman–Galley formalism for non-conservative
classical mechanics (Bateman, 1931; Galley, 2013; Aykroyd
et al., 2025), can be found in (Scurria, 2026).

On the diagonal manifold z = z ′ the doubled system enjoys
a gauge symmetry: the auxiliary variable z ′ is redundant and
the difference d is identically zero. Credit assignment is implemented by deliberately breaking this symmetry through
the task cost. Adding βD(z, z ′ ) = β C(m) to H exerts
opposite forces on z and z ′ and drives them apart, so that
the difference d ceases to be redundant and begins to carry
information about the loss landscape.

C.1. The Helmholtz Obstruction
The natural physical route to a variational principle for a
dynamical system ẋ = F (x, θ) is to seek a scalar potential
E such that F = −∂x E. The classical Helmholtz integrability condition states that such an E exists if and only if the
Jacobian JF is symmetric everywhere. Whenever the interactions are non-reciprocal — as in feedforward networks,
active matter, or driven optical systems — JF acquires
a non-zero antisymmetric part and the Helmholtz condition fails identically. No scalar potential on the original
n-dimensional state space can then generate the dynamics,
and the “energy minimisation” route at the heart of standard
EP is blocked at the structural level. The obstruction is not
a matter of computational convenience: it reflects the fact
that the rotational component of F carries information that
no scalar function of x alone can record.

D. Proofs for Dyadic EP
We now demonstrate that Dyadic EP correctly trains the
parameters θ of the original force field F (x, θ), giving the
0
)
in the limit of infinitesimal nudging.
exact gradient dC(x̄
dθ
D.1. Proof of EP
First, recall that standard EP does not strictly require the
system to settle at an energy minimum; it requires only that
the system reaches a stationary state (a fixed point of the
dynamics). Indeed, using the notation of Section 2.1, EP
relies on the key identity:
d2
d2
ET (xβ , θ) =
ET (xβ , θ).
dθdβ
dβdθ

C.2. Variational Reconstruction on a Doubled Space

Expanding the total derivative with respect to β gives:

Applying the Bateman–Galley formalism circumvents this
obstruction by enlarging the configuration space. The single
state x ∈ Rn is replaced by a conjugate pair (z, z ′ ) ∈ R2n ,
and the rotational component of F — which has no scalar
generator on the original n-dimensional space — is absorbed into a bilinear coupling between z and z ′ on the
doubled space, where it does admit a variational description. The physical motion is recovered on the diagonal
submanifold z = z ′ (the so called ’physical limit’), while
the off-diagonal direction d = z − z ′ supplies the additional
degree of freedom needed to encode non-reciprocity.

d
ET (xβ , θ) =
dβ




z + z′
,θ ,
2



∂ET (xβ , θ)
∂x

⊤

dxβ
∂ET (xβ , θ)
+
dβ
∂β

= C(xβ ).

(58)

Where the first term vanishes because the system is at a
∂
ET (xβ , θ) = 0; this holds even if
stationary state, i.e., ∂x
the system is not at a minimum of ET . Similarly, for the
derivative with respect to θ:
d
∂ET (xβ , θ)
ET (xβ , θ) =
,
dθ
∂θ

Specializing this reconstruction to the overdamped (firstorder) regime relevant to relaxational neural dynamics yields
the bilinear energy
H(z, z ′ , θ) = −(z − z ′ )⊤ F

(57)

(59)

where we additionally assume that the cost function does
not depend explicitly on the parameters θ. Substituting these
results into Eq. (57) in the limit of infinitesimal nudging
(β → 0) recovers the fundamental relation given by Eq. (9).

(56)
D.2. Proof of Dyadic EP

which is precisely Eq. (26). The symmetric midpoint m =
(z + z ′ )/2 plays the role of the physical coordinate of the
doubled system, while d is the auxiliary direction along
which non-reciprocity is stored. On the submanifold z = z ′
the coupling proportional to (z − z ′ ) vanishes identically
and both states evolve under the original field F , so the
doubling leaves the on-shell physics unchanged. We refer
the reader to (Scurria, 2026) for the full construction.

We analyze now the stationary states of Dyadic EP by introducing the change of variables:
m=

z + z′
,
2

d = z − z′.

(60)

In these coordinates, the augmented energy HT becomes
HT (m, d, θ, β) = −d⊤ F (m, θ) + βC(m)
14

(61)

Equilibrium Propagation for Non-Conservative Systems

and the dynamics in Eq. (28) can be rewritten as:
∂HT
dm
=−
= F (m, θ),
dt
∂d
dd
∂HT
∂
=−
= dT JF (m, θ) − β
C(m).
dt
∂m
∂m

In Dyadic EP, we instead employ the single-phase update:


1 ∂H(z β , z ′β , θ)
∆θ ∝ −
(70)
β
∂θ

(62)
(63)

This choice avoids the overhead of evolving two coupled
equations in the extended space, which would be computationally equivalent to evolving four equations in the original
space (two for +β and two for −β). Using Eq. (70), we
evolve only one coupled equation for +β in the extended
space; this corresponds to two equations in the original
space, thereby achieving the same computational complexity as AsymEP. Furthermore, this single-phase formulation
suggests a pathway toward making the update local in time,
provided appropriate hardware is used to implement the
augmented phase.

β

The stationary states (mβ , d ) are the solutions to:
F (mβ , θ) = 0,
∂
βT
d JF (mβ , θ) − β
C(mβ ) = 0.
∂m

(64)
(65)

This leads to the following observations:
1) The stationary state of m is independent of β and coincides with the stationary state of the original system:
z β + z ′β
= mβ = m0 = x0 .
2

Mathematically, these two approaches yield the same gradient estimate because the equations for dβ are linear. Explicitly we have :


∂H(z β , z ′β , θ)
∂F z β + z ′β
= −(z β − z ′β )⊤
,θ
∂θ
∂θ
2

⊤

−1
∂F 0 

= −β
z ,θ
JF⊤ (z 0 , θ)
∂θ


∂C 0
×
(z ) ,
(71)
∂x

(66)

2) The Jacobian of the extended system defined in Eq. (26)
is invertible, provided JF is invertible. This is most evident
from Eq. (63).
3) The stationary state value of d is given by:



−1 ∂C 0
β
d = β JF⊤ (m0 , θ)
(x )
∂x

(67)

where we have used Eqs. (66) and (67). Inspection of
Eq. (71) confirms that, up to corrections of order β 2 , we
obtain exactly the same gradient as in AsymEP.

0

In particular, when β = 0, we have d = 0, which implies
that the free stationary states coincide: z 0 = z ′0 .

E. AsymEP versus Dyadic EP

4) The cost at the stationary state of the extended system
is equal to the cost at the stationary state of the original
system:
D(m0 ) = C(x0 ).
(68)

In this appendix, we demonstrate that Asymmetric Equilibrium Propagation (AsymEP) emerges naturally as the firstorder projection of the 2N -dimensional Dyadic Equilibrium
Propagation onto a single N -dimensional state space. We
then formalize the physical trade-offs between the two architectures.

Consequently, the gradients of the cost with respect to the
parameters are identical.
Since both the original and extended systems, given respectively in Eq. (28) and Eq. (1-2), share the same cost at their
respective stationary states, and because the Jacobians of
both models are invertible, applying EP update rule to the
extended system give the correct gradient estimate for the
parameters θ of the original system.

E.1. AsymEP as the Linear Projection of Dyadic EP
As established in Appendix D.2, transforming the 2N dimensional extended space (z, z ′ ) into the mean state
′
m = z+z
and the difference state d = z − z ′ exactly
2
decouples the stationary dynamics. Because the stationary
state of m is the free state of the original system (mβ = x0 ),
the cost function drives the difference variable to a stationary
β
state d satisfying:

The final step of the proof is to establish the equivalence
between the standard parameter update rule in Eq. (8) and
the modified rule used by Dyadic EP in Eq. (34). Indeed, if
we were to apply the standard update rule in the extended
space, the update would be:
1
∆θ ∝ −
2β



β

′β

−β

∂H(z , z , θ) ∂H(z , z
−
∂θ
∂θ

′−β

, θ)

β

JF⊤ (x0 , θ)d = β


.

∂C 0
(x )
∂x

(72)

To recover this exact error signal in an N -dimensional space,
we postulate a modified dynamical system FA (x) compris-

(69)
15

Equilibrium Propagation for Non-Conservative Systems

F. Derivation of the Hopfield-like Energy

ing the standard EP dynamics and a spatial correction Γ(x):
FA (x) = F (x) − β

∂C
(x) + Γ(x)
∂x

In this section, we derive the explicit energy functional for
the Continuous Asymmetric Hopfield dynamics defined in
Eq. (35). The force field is given by:

(73)

Let ∆x = xβA − x0 denote the displacement from the
free equilibrium. Expanding the stationarity condition
FA (xβA ) = 0 to first order around x0 yields:
JF (x0 , θ)∆x − β

∂C 0
(x ) + Γ(xβA ) ≈ 0
∂x

F (x) = ρ′ (x) ⊙ (Jρ(x)) − x.

We omit external inputs J in for brevity, as they appear symmetrically in the Jacobian. The variational Hamiltonian is
defined as:




z + z′
z + z′
H(z, z ′ ) = −(z − z ′ )⊤ F
+ βC
.
2
2
(79)

(74)

To ensure the first-order displacement matches the Dyadic
β
EP error signal (i.e., ∆x ≈ d ), we substitute Eq. (72) into
the expansion:


Γ(xβA ) = JF⊤ (x0 , θ) − JF (x0 , θ) ∆x

(75)

= −2AJ (x0 , θ)(xβA − x0 )

(76)

To analyze this expression, we introduce the midpoint m =
z+z ′
and the difference d = z − z ′ . Since the separation
2
between z and z ′ is induced solely by the nudging parameter
β, the difference scales as ∥d∥ ∼ O(β). We therefore
neglect terms of order O(∥d∥3 ) (i.e., or equivalently O(β 3 ))
as they do not contribute to the gradient of the cost.

This uniquely recovers the AsymEP augmented dynamics.
Finally, to eliminate the O(β 2 ) error, AsymEP evaluates the
centered difference of two opposite nudges:
0
x±β
A =x ±β

dxA
+ O(β 2 )
dβ β=0

(78)

The activation at the midpoint can be approximated as:
ρ(m) =

(77)

ρ(z) + ρ(z ′ )
+ O(∥d∥2 ).
2

(80)

Similarly, the difference in activations is:

Subtracting these states cancels the O(β 2 ) error, yielding
β
+β
−β
1
3
2 (xA − xA ) = d + O(β ), successfully recovering the
exact post-synaptic update term.

ρ(z) − ρ(z ′ ) = ρ′ (m) ⊙ d + O(∥d∥3 ).

(81)

Inverting this relation, we express the state difference as:
E.2. Physical Trade-offs and the Extended Space

z − z ′ = (ρ(z) − ρ(z ′ )) ⊙ ρ′ (m) + O(∥d∥3 ).

We can view AsymEP and Dyadic EP as a space-time tradeoff of the same underlying physical optimization problem.

(82)

We substitute these expansions into the interaction term
of the Hamiltonian, Hint = −(z − z ′ )⊤ (ρ′ (m) ⊙ Jρ(m)).
Applying the identity a⊤ (b ⊙ c) = (a ⊙ b)⊤ c, we obtain:

AsymEP preserves the original N -dimensional state space
of the network at the cost of temporal non-locality. The system must evolve sequentially, requiring physical memory
not only to store the free equilibrium x0 for the asymmetric correction, but also to store the successive stationary
states required to evaluate the contrastive gradient update.
AsymEP thus serves as the direct, spatially minimal extension of EP.

⊤

Hint = − ((z − z ′ ) ⊙ ρ′ (m)) Jρ(m)


ρ(z) + ρ(z ′ )
≈ −(ρ(z) − ρ(z ′ ))⊤ J
.
2

(83)

Expanding the product gives:

Dyadic EP provide a learning signal that is local in both
space (where z − z ′ encodes the gradient) and time (allowing the nudged phases to execute in parallel) at the cost
of doubling the state space. In particular, capturing nonconservative forces in this extended space requires a specific bilinear coupling, rather than a trivial superposition
of uncoupled subsystems. It can be seen as a blueprint for
future neuromorphic hardware.

Hint = −

1h
ρ(z)⊤ Jρ(z) + ρ(z)⊤ Jρ(z ′ )
2
i

− ρ(z ′ )⊤ Jρ(z) − ρ(z ′ )⊤ Jρ(z ′ ) .

(84)

We decompose the connectivity matrix J into its symmetric
part S and antisymmetric part A. The first and last terms
simplify to ρ(z)⊤ Sρ(z). The cross terms satisfy:

Ultimately, the reduction of Dyadic EP to AsymEP via the
variables m and d proves the universality of EP’s variational
principle.

ρ(z)⊤ Jρ(z ′ ) − ρ(z ′ )⊤ Jρ(z) = ρ(z)⊤ (J − J ⊤ )ρ(z ′ )
= ρ(z)⊤ (2A)ρ(z ′ ).
16

(85)

Equilibrium Propagation for Non-Conservative Systems

Thus, the interaction term reduces to:
1
1
Hint = − ρ(z)⊤ Sρ(z) + ρ(z ′ )⊤ Sρ(z ′ )
2
2
− ρ(z)⊤ Aρ(z ′ ) + O(∥d∥3 ).

The input parameters are then updated using the standard
learning rule (21). In particular, the presynaptic term associated with the input weights is given by,
∂Fi
= δik ρ′ (xi )ul .
in
∂Jkl

(86)

Finally, for the nudging term, we expand the cost function
around the midpoint:
C(m) =

1
(C(z) + C(z ′ )) + O(∥d∥2 ).
2

The presynaptic terms associated with the dynamical paramdyn
eters Jij
depend on the experiment.
G.1. Symmetric Initialization

(87)

G.1.1. L EARNING RULES

When multiplying by β, the remainder term becomes β ·
O(∥d∥2 ). Since ∥d∥ ∼ O(β), this remainder is of order
O(β 3 ) and can be consistently discarded alongside the thirdorder terms from the interaction expansion.

For clarity, we write the learning rules for VF and AsymEP.
For the input weights, using (93), we have:
i
1 h +β
in
′ 0
∆Jik
∝
(xi − x−β
(94)
i )ρ (xi )uk ,
2β

Combining all these components, the final Hamiltonian is:

while for the recurrent weight, we get:
i
1 h +β
dyn
0
′ 0
(xi − x−β
)ρ(x
)
.
∆Jij
)ρ
(x
∝
i
j
i
2β

1
1
H(z, z ′ ) = − ρ(z)⊤ Sρ(z) + ρ(z ′ )⊤ Sρ(z ′ )
2
2
1
⊤
′
− ρ(z) Aρ(z ) + (∥z∥2 − ∥z ′ ∥2 )
2
β
′
+ (C(z) + C(z )).
(88)
2

in
∆Jik
∝

(89)
(90)

To complement Fig. 2, we report the evolution of the accuracy of the three methods in Fig. 4. We consider a layered
network with 50 hidden neurons. While this capacity is
insufficient for state-of-the-art performance, it amplifies the
difference in accuracy between models to aid visualization.
Models are trained for 20 epochs starting from a symmetric
configuration, the natural setting for both VF and EP. With
this initialization, AsymEP consistently outperforms the
other methods and learns faster by exploiting the additional
degrees of freedom of the asymmetric network.

As in the main text, the neuronal dynamics are governed by
the vector field:


X dyn
Fi = ρ′ (xi ) 
Jij ρ(xj ) + bi (u) − xi ,
(91)
j

G.2. Fixed Asymmetry Ratio

where the input-dependent bias bi (u) is precomputed for
each MNIST input u as:
Jilin ul .

(96)

G.1.2. S UPPLEMENTARY N UMERICAL R ESULTS

G. Experimental Details

bi (u) =

 i
1 h +β
ρ(xi ) − ρ(x−β
)
uk ,
i
2β

and for the recurrent weights:
i
1 h +β
dyn
−β
−β
∆Jij
∝
ρ(xi )ρ(x+β
j ) − ρ(xi )ρ(xj ) . (97)
2β

This system recovers the original continuous Hopfield dynamics when z = z ′ (assuming β = 0).

X

(95)

For EP, we have:

The saddle-point dynamics, given by Eq. 32, generated by
this Hamiltonian are:
β ∂C
dz
= ρ′ (z) ⊙ (Sρ(z) + Aρ(z ′ )) − z −
,
dt
2 ∂z
′
dz
β ∂C
= ρ′ (z ′ ) ⊙ (Sρ(z ′ ) + Aρ(z)) − z ′ +
.
dt
2 ∂z ′

(93)

This section details the implementation for the fixed asymmetry ratio experiments presented in Section 5.2, followed
by complementary numerical results regarding learning
speed and induced Jacobian asymmetry.

(92)

l∈in

This term projects the input space into the recurrent subspace. The bias yields a diagonal contribution to the Jacobian JF = ∂F
∂x , and therefore does not contribute to the
antisymmetric correction used in the augmented dynamics
Eq. (20) of AsymEP.

G.2.1. L EARNING RULES
Parametrization and notation. To enforce a fixed asymmetry ratio, we explicitly parameterize the independent elements of Eq. (38). We introduce two parameter vectors θS
17

Equilibrium Propagation for Non-Conservative Systems

Parameter

Sym. Init. / Feedforward
sec. 5.1 & 5.3

Fixed rstr
sec. 5.2

Fixed rstr & rin
app. G.3

0.05
0.01
0.5
0.5
20
10
40 / 20
64
n.a.
784 - n.a. -10
tanh
s ∼ U (−1, 1)
θ ∼ N (0, N1 )
10

0.05
0.01
0.3
0.5
30
10
30
√64
60
784-50-10
tanh
s ∼ U(−1, 1)
θ ∼ N (0, N1 )
10

0.0125
0.0025
0.3
0.5
40
10
40
√64
60
all-to-all, 500 hid
tanh
s ∼ U(−1, 1)
θ ∼ N (0, N1 )
10

Learning Rate (Input-Hidden)
Learning Rate (Hidden-Output)
Time Step (Dynamics Integration)
Nudging Parameter (β)
Free-phase Steps (nfree )
Nudged-phase Steps (nnudge )
Number of Epochs
Batch Size
Scaling Parameter γ
Structure
Activation function ρ
Initial Recurrent State s
Initial Parameters θ
Number of Runs (training + inference)

Table 3. Trained Model Hyperparameters on MNIST. N is the total number of neurons, U(−1, 1) is a uniform distribution, and N (µ, σ 2 )
is a Gaussian distribution. For the rstr parametrization, we choose more cautious hyperparameters for training and inference compared to
the symmetric initialization, due to increasingly non-conservative and potentially oscillatory dynamics.

elements of S̃, the full matrices are constructed as:
S
S̃ij = δij ξi + (1 − δij )θk(max(i,j),min(i,j))
,
A
Ãij = ϵij θk(max(i,j),min(i,j))
,

(99)
(100)

where ϵij is the Levi-Civita symbol. The dynamical parameters are then given by:
dyn
Jij
= γ(cS S̃ij + cA Ãij ),

with normalization coefficients
p
2
1 − rstr
cS =
,
FS

cA =

(101)

rstr
,
FA

defined in terms of the Frobenius norms:
v
uN
M
uX
X

2
F =t
ξ2 + 2
θS ,

Figure 4. Evolution of the mean accuracy and standard deviation
(over 10 runs) during training on MNIST for AsymEP, EP, and VF.
Models use 50 hidden neurons.

S

i

i=1

k

(102)

(103)

k=1

v
u M
u X

2
θkA .
FA = t2

(104)

k=1

and θA of size M = Ndyn (Ndyn − 1)/2, which encode the
off-diagonal elements of the symmetric and antisymmetric
components S̃ and Ã, respectively. The correspondence
between matrix and vector indices is given by:

k(i, j) =

(i − 1)(i − 2)
+ j,
2

Presynaptic computation. The dependence of the normalization coefficients on the parameters introduces additional regularization terms in the learning rule compared
to the parameterization of (Scellier &Bengio, 2017). The
gradients of the normalization coefficients are:

(1 ≤ j < i ≤ Ndyn )
(98)

∂cS
θkS
=
−2c
S
2,
∂θkS
(FS )

where the condition j < i selects the strictly lower triangular
elements. Introducing an additional vector ξ for the diagonal

∂cA
θA
= −2cA k 2 .
A
∂θk
(FA )
18

∂cS
ξm
= −cS
2,
∂ξm
(FS )

(105)
(106)

Equilibrium Propagation for Non-Conservative Systems

Parameter

Comparison Dyn.
sec. 5.4

2 hidden layers
sec. 5.4

3 hidden layers
sec. 5.4

Learning Rate (Input-Hidden)
Learning Rate (Hidden-Hidden)
Learning Rate (Hidden-Output)
Time Step (Dynamics Integration)
Nudging Parameter (β)
Free-phase Steps (nfree )
Nudged-phase Steps (nnudge )
Number of Epochs
Batch Size
Layer Structure
Activation function ρ
Initial Recurrent State s
Initial Parameters θ
Number of Runs (training + inference)

0.0016
0.0016
0.0016
0.4
0.3
40
20
50
64
784-500-200-10
tanh
s ∼ U(−1, 1)
θ ∼ N (0, N1 )
10

0.0013
0.0013
0.0013
0.3
0.5
40
20
40
64
784-500-500-10
tanh
s ∼ U(−1, 1)
θ ∼ N (0, N1 )
10

0.6
0.6
0.6
0.0075
0.20
60
30
40
64
784-500-500-500-10
tanh
s ∼ U (−1, 1)
θ ∼ N (0, N1 )
10

Table 4. Trained Model Hyperparameters on Fashion-MNIST. N is the total number of neurons, U(−1, 1) is a uniform distribution, and
N (µ, σ 2 ) is a Gaussian distribution. For the rstr parametrization, we choose more cautious hyperparameters for training and inference
compared to the symmetric initialization, due to increasingly non-conservative and potentially oscillatory dynamics.

Combining these with the derivatives of the matrices S̃ and
Ã, we have:
∂ S̃ij
= δip δjq + δiq δjp ,
∂θkS

∂ S̃ij
= δij δkj
∂ξk

∂ Ãij
= δip δjq − δiq δjp ,
∂θkA

(where p > q):

N
θkA X
∂Fi
′
=
γc
ρ
(x
)
−2
Ãij ρ(xj )
A
i
(FA )2 j=1
∂θkA

(107)


+ δip ρ(xq ) − δiq ρ(xp ) .

(108)

(111)

where k corresponds to the index pair (p, q) with p > q, as
defined in Eq. (98). The full presynaptic terms are then:

Initialization. To ensure the stability of the system, we
initialize our parameters suchhthat the
i variance of dynamdyn
ical parameters scales as Var Jij ∝ 1/Ndyn . This is a
conservative choice for the layered architectures used in our
dyn
experiments, where many entries of Jij
are zero.

• For the diagonal parameters ξm :

N
∂Fi
ξm X
S̃ij ρ(xj )
= γcS ρ′ (xi ) −
∂ξm
(FS )2 j=1

+ δim ρ(xm ) .

In practice, we initialize the parameter vectors θS , θA , and
ξ with identical variances σ 2 . For large Ndyn , the expected
Frobenius norms approximate to E[FS,A ] ≈ Ndyn σ. Consequently, the normalization coefficients become:

(109)

p

• For the off-diagonal symmetric parameters θkS (where
p > q):

cS ≈


N
∂Fi
θkS X
′
=
γc
ρ
(x
)
−2
S̃ij ρ(xj )
S
i
(FS )2 j=1
∂θkS

2
1 − rstr
,
Ndyn σ

cA ≈

rstr
.
Ndyn σ

(112)

Since the symmetric and antisymmetric components are statistically independent, the variance of the weights is derived
as follows:


+ δip ρ(xq ) + δiq ρ(xp ) .

• Diagonal elements (i = j):

(110)
h
i
2
1 − rstr
Var Jiidyn = γ 2 c2S σ 2 ≈ γ 2
.
2
Ndyn

• For the off-diagonal antisymmetric parameters θkA
19

(113)

Equilibrium Propagation for Non-Conservative Systems

• Off-diagonal elements (i ̸= j):
i
h


γ2
dyn
= γ 2 c2S + c2A σ 2 ≈ 2 ,
Var Jij
Ndyn

a zero-cost baseline (perfect prediction) during learning.
Specifically, for each method and value of rstr , we calculate the cumulative loss by summing the batch-averaged
costs of the first 5 epochs (out of 30, to avoid saturation
effects), and reporting the mean and standard deviation over
10 independent training runs. Mathematically, for each run:

(114)

i
h
dyn
∝ 1/Ndyn , we set:
To satisfy Var Jij
γ=

p
Ndyn

(115)
Cumul. Loss =

Note that by random matrix theory, diagonal elements do
not affect stability in the large Ndyn limit.

(116)

Ãij = ϵij θk(max(i,j),min(i,j)) .

(117)

NX
batches


X


epoch=1 k=1

(x0 ,u)∈Bk


C(x0 , u) 
,
|Bk |

(120)
where Bk represents the k-th batch, and |Bk | denotes the
number of examples in the batch. The parameters are updated after each batch step; consequently, the free equilibrium x0 is inferred using the updated parameters and the
current example u.

Potential Simplification. Although the parameterization
above is fully general, a simpler construction is possible
by removing self-connections (ξ = 0) and enforcing identical parameterization for the symmetric and antisymmetric
components, i.e., θS = θA = θ. The matrix elements then
become:
S̃ij = (1 − δij )θk(max(i,j),min(i,j)) ,

5
X

In this case, the Frobenius norms are equal (FS = FA ), and
we can omit the explicit normalization:
q
dyn
2 S̃ + r à .
Jij
= 1 − rstr
(118)
ij
str ij
For a parameter θk corresponding to indices (p, q) with
p > q, the presynaptic term is given by:
q

∂Fi
2 +r
1 − rstr
= ρ′ (xi )
str δip ρ(xq )
∂θk
q

 (119)
2 −r
1 − rstr
+
δ
ρ(x
)
str
iq
p .
While this parameterization works in simulations and keeps
the number of parameters constant for all rstr , it constrains
the asymmetry to be “homogeneous”, by which we mean
that the asymmetry ratio is identical for every pair of neurons; hence, the network cannot learn to be symmetric in one
region and antisymmetric in another. Therefore, we choose
to explore the more general case of (38) in our experiments.

Figure 5. Cumulative loss as defined by (120) over the first 5
epochs of training, for different asymmetry ratios rstr . We compare
VF (orange) and AsymEP (blue), under two training regimes:
training only J in (dashed) or all parameters (solid).

In Fig 5, we observe that learning slows down for both algorithms when rstr ≳ 0.6. This behavior likely results from
the increased difficulty of reaching a stationary state as the
dynamics become strongly asymmetric. With a fixed number of inference steps, incomplete convergence degrades the
accuracy of the gradient estimates, thereby slowing down
the learning. Fig 5 shows that while VF can eventually
achieve competitive accuracy, it is consistently slower than
AsymEP as soon as asymmetry is introduced.

G.2.2. S UPPLEMENTARY N UMERICAL R ESULTS
To complement the results of Fig 3, we analyze the training
efficiency as a function of the asymmetry ratio rstr and investigate the robustness of VF by monitoring the Jacobian
asymmetry.
Training efficiency. We first study the training efficiency
of the two algorithms as a function of the asymmetry ratio rstr . Inspired by the related concept in (Cesa-Bianchi
&Lugosi, 2006), we define the cumulative loss as the accumulated difference between the free equilibrium cost and

Jacobian asymmetry. We next examine how the structural asymmetry rstr is reflected in the Jacobian of the dy20

Equilibrium Propagation for Non-Conservative Systems

namics (35), given by:
∂Fi
dyn ′
= (1 − δij )ρ′ (xi )Jij
ρ (xj )
∂sj
h
i
+ δij ρ′ (xi )(Jiidyn ρ′ (xi )) + ρ′′ (xi )bi − 1 .
(121)
In our layered architecture, the self-connections are zero
(Jiidyn = 0). For the following analysis, we neglect all diagonal terms in the Jacobian (including external inputs and
potential), since they do not contribute to the antisymmetric
correction (20) and thus to the discrepancy between the performance of VF and AsymEP. Consequently, we define the
following asymmetry ratio based solely on the off-diagonal
Jacobian JF,off :
(122)

Figure 6. Asymmetry ratio of the Jacobian rjac defined in equation
(122) after training for different asymmetry ratios rstr . We compare
VF (orange) and AsymEP (blue), under two training regimes:
training only J in (dashed) or all parameters (solid).

The results are presented in Fig 6. For each trained model
and ratio rstr , we compute rjac averaged over the stationary
states of the first batch (64 images) across 10 independent
runs. We observe that when structural asymmetry is strong
and all parameters are trained, VF partially compensates for
the asymmetry by adjusting the neuronal states. This can be
understood by rewriting the ratio as:


dyn
dyn ⊤
ρ′ (xi ) Jij
− (Jji
) ρ′ (xj )
F
rjac =
.
(123)
dyn ′
ρ′ (xi )Jij
ρ (xj )

Consequently, local stability requires the largest real eigenvalue of the effective weight matrix to be strictly less than 1.
Assuming weights are initialized independently with variance σ 2 , Girko’s circular law dictates that the eigenvalues
of√an asymmetric matrix uniformly populate a disk of radius
σ n in the complex plane. In contrast, imposing symmetry
forces the eigenvalues
√ onto the real line, broadening the
spectral radius to 2σ n according to Wigner’s semicircle
law. As a result, asymmetric networks can stably accommodate larger variance in the weight initializations than their
symmetric counterparts.

Compared to the structural asymmetry ratio in Eq. (37),
a value of rjac < rstr indicates that the neuronal states effectively dampen the structural asymmetry, rendering the
dynamics more symmetric. This symmetrization of the Jacobian appears without imposing an additional symmetrization penalty and could be enhanced using the method of
(Laborieux &Zenke, 2022). This mechanism likely explains
the superior performance of ‘All (VF)’ compared to ‘J in
(VF)’ in Fig 3, as the former is able to use the additional
degrees of freedom to reduce the effective asymmetry at
high rstr .

Asymmetry nevertheless introduces imaginary eigenvalues
and, consequently, damped oscillations. To study this effect
experimentally in a controlled setting, we constrain the input
projections J in . In the experiments of the main text, fixing
the structural asymmetry ratio rstr still allowed AsymEP
to reduce oscillations by aligning and increasing the input
projections J in , thereby adding stabilizing diagonal contributions to the Jacobian. To isolate the network’s ability to
suppress oscillations independently of the magnitude of the
input drive, we further constrain the relative scale of J in and
J dyn by imposing

rjac =

⊤
∥JF,off − JF,off
∥F
,
∥JF,off ∥F

F

rin =

G.3. Stability analysis with Fixed Asymmetry Ratio &
Constrained Inputs Projection

∥J in ∥F
∥J in ∥F
=
,
∥J dyn ∥F
γ

(124)

where ∥J dyn ∥F = γ following Eq. (101). Defining unscaled
input projections J˜in , we set

A complete stability analysis of the non-conservative dynamics trainable with AsymEP is beyond the scope of this
work. Nevertheless, for the class of continuous Hopfield
networks considered here, standard arguments from random
matrix theory suggest that asymmetry inherently improves
asymptotic stability.

J in = rin γ

J˜in
∥J˜in ∥F

(125)

G.3.1. L EARNING RULES

In the dynamics defined by Eq. (91), the linear leak term
−xi shifts the spectrum of the system’s Jacobian by −1.

Reusing the notations of the previous section, we write
Jilin = γcin J˜ilin with the normalization cin = rin /Fin , where
21

Equilibrium Propagation for Non-Conservative Systems

Fin = ∥J˜in ∥F . Applying the chain rule yields:
"
#
in X
∂Fi
J˜kl
′
in
˜
= γcin ρ (xi ) δik ul − 2
J um .
in
Fin m im
∂ J˜kl

(126)

And for γ we have:
1
∂Fi
= (Fi + xi ).
∂γ
γ

(127)

G.3.2. S UPPLEMENTARY N UMERICAL R ESULTS
Figure 7. Comparison of AsymEP and VF on a feedforward network. Test accuracy on MNIST is shown as a function of training
epochs for a single-hidden-layer network with 20 neurons. Curves
report the mean and standard deviation over 10 runs. Best accuracies are 92.7% ± 0.5% (AsymEP) and 64.3% ± 2.0% (VF).

Table 5 reports a worst-case control experiment where the
structural asymmetry is fixed at rstr = 0.7 while the input
scale ratio rin is varied. The experiment uses an all-to-all
architecture on MNIST (excluding direct input-to-output
connections). The output variance during extended inference (steps 30-50) confirms that the system successfully
learns to suppress oscillations even when rin is severely restricted. Any small residual oscillations can be mitigated by
time-averaging over the inference steps.

G.5. Advantages of Non-Conservatives Dynamics
In Section 5.4, we compare three (non-)conservative dynamics under varying constraints. To further evaluate learning
speed, Table 6 reports network performance after a single epoch. These results confirm our earlier observation:
AsymEP learns faster than VF.

Finally, rin can be interpreted as a measure of the external
signal magnitude relative to the internal recurrent dynamics.
These results indicate that the system remains capable of
learning and stabilizing even under a low external input
drive. Even when the input projection ∥J in ∥F is 100 times
smaller than the recurrent connections ∥J dyn ∥F , the network
still achieves 36.34 ± 6.25% accuracy, which is well above
chance level (∼ 10%).

G.6. Feedforward CIFAR-10 Experiments
This appendix details the architecture and hyperparameters
of the deep feedforward experiments comparing backpropagation (BP), VF, AsymEP and Dyadic EP on CIFAR-10
(see subsection 5.5)

G.4. Feedforward Network
Architecture. We use a nine-layer convolutional network
(denoted CNN9). The first eight layers are convolutional
with 3 × 3 kernels and zero-padding; spatial downsampling
is performed by strided convolutions (stride 2 on layers 2, 4,
6, 8 and stride 1 otherwise), so no pooling is used. The channel widths follow the sequence 3 → 64 → 64 → 128 →
128 → 256 → 256 → 512 → 512, reducing the spatial
resolution from 32 × 32 to 2 × 2. The last layer is a fully
connected readout mapping the 512 × 2 × 2 feature map
to the 10 class logits. All hidden units use a ReLU nonlinearity.
p Weights are initialized with the Kaiming scheme
(σ = 2/fan-in) and biases at zero.

G.4.1. L EARNING RULES
For clarity, we write the learning rules for VF and AsymEP
in a feedforward architecture with one hidden layer using
the notation of Section 5.3. For the input weights connecting
to the hidden layer, we get the usual formula:
in
∆Jik
∝

i
1 h +β
−β
0
(hi − hi )ρ′ (hi )uk ,
2β

(128)

while for the feedforward weights connecting the hidden to
the output layer, we get:
∆(Wh→o )ji ∝

i
1 h +β
0
′ 0
(oj − o−β
j )ρ (oj )ρ(hi ) .
2β

Training setup. All methods are trained for 40 epochs
with batch size 64 and repeated over 5 seeds. Inputs are
normalized per channel and augmented with random 32 ×
32 crops (padding 4), random horizontal flips and Cutout
(one 8 × 8 patch). Parameters are updated with SGD with
momentum 0.9, weight decay 5 × 10−4 and gradient-norm
clipping at 1, under a cosine learning-rate schedule decaying
from 3.5 × 10−2 to 2 × 10−4 . Test accuracy is reported
on an exponential moving average of the weights (decay

(129)

Note that EP is not applicable in this case.
G.4.2. S UPPLEMENTARY N UMERICAL R ESULTS
In addition to the final accuracy reported in Sec. 5.3, we
show in Fig. 7 the evolution of the accuracy over 20 epochs
for AsymEP and VF.
22

Equilibrium Propagation for Non-Conservative Systems
Table 5. Output variance and final test accuracy on MNIST (%) across different values of rin with rstr = 0.7. (mean ± std over 10 runs)
(500 hiddens, all-to-all, no input-output).

rin
0.01
0.10
0.50
1.00

Output variance
Untrained
Epoch 80
(3.38 ± 0.90) × 10−4
(2.33 ± 0.48) × 10−4
(1.34 ± 0.32) × 10−5
(6.27 ± 1.24) × 10−7

(5.22 ± 2.34) × 10−5
(1.39 ± 0.17) × 10−4
(1.06 ± 0.25) × 10−6
(1.75 ± 0.50) × 10−8

Table 6. Test accuracy on Fashion-MNIST (%) at Epoch 1 (mean
± std 10 runs). The table compares three classes of network
dynamics: Continuous Hopfield (CH), Predictive Coding (PC),
and Standard dynamics. Each is evaluated under three connec⊤
tivity structures: Asymmetric (Asym, Bk ̸= Wk+1
), Symmet⊤
ric/conservative (Sym, Bk = Wk+1 ), and Feedforward (Feedfor,
Bk = 0).
EP

AsymEP

Asym
74.91 ± 0.45
Feedfor
74.36 ± 0.29
Sym
74.57 ± 0.43
Asym
77.83 ± 0.47
PC
Sym
76.23 ± 0.39
Asym
76.87 ± 0.51
Standard
Feedfor
77.92 ± 0.51
CH

VF
48.98 ± 4.09
48.84 ± 3.46
57.75 ± 3.37
61.50 ± 4.06
63.98 ± 0.73

0.9995). The targets are smoothed (ε = 0.1), which for
the EP methods amounts to nudging toward the smoothed
one-hot vector.
Relaxation hyperparameters. The four methods differ
only in the gradient estimator: BP uses automatic differentiation, while the EP-based methods contrast stationary
states of the corresponding relaxation dynamics. VF uses
an integration step η = 1.0, nudging β = 0.1, and at most
K = 1000 relaxation steps with an early-stopping threshold
of 9 × 10−6 on the mean state update. Dyadic EP uses
the same settings except for a nudging strength β = 0.1.
AsymEP uses a smaller step η = 0.5, nudging β = 0.1,
and up to K = 250 relaxation steps with a threshold of
1 × 10−4 .

23

Test Acc. (%)
Epoch 80
36.34 ± 6.25
90.54 ± 0.19
94.96 ± 0.10
96.30 ± 0.09