Initial commit: ept — backprop-free equilibrium transformer (EP)

Code (ep_run/), organized docs (docs/{method,campaign,hardware,outreach,paper}),
analysis scripts (scripts/), ONBOARDING.md entry point. Large data/checkpoints
git-ignored (share separately).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_014FAPDWQ49M5Ye3NpTndTpn

1 files changed, 3311 insertions, 0 deletions

diff --git a/refs/scurria_2602.03670v2.txt b/refs/scurria_2602.03670v2.txt
new file mode 100644
index 0000000..d9eba6a
--- /dev/null
+++ b/refs/scurria_2602.03670v2.txt
@@ -0,0 +1,3311 @@
+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
+
+
\ No newline at end of file