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 . 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