Equilibrium Propagation and Hamiltonian Inference in the Diffusive Fitzhugh-Nagumo Model

# Equilibrium Propagation and Hamiltonian Inference in the Diffusive Fitzhugh-Nagumo Model
Source: [https://arxiv.org/html/2605.21568](https://arxiv.org/html/2605.21568)
###### Abstract

In this work, we extend the Equilibrium Propagation framework to skew\-gradient systems and show an equivalence between deep Energy\-Based Models and Hamiltonian neural networks\. We focus on networks of diffusively coupled Fitzhugh\-Nagumo neurons as a prototypical example\. We show that since stationary solutions of the Fitzhugh\-Nagumo model are described by self\-adjoint operators, the methods of equilibrium propagation for performing credit assignment can be applied\. Furthermore, for Fitzhugh\-Nagumo networks with the topology of a deep residual network, we show that the steady state solutions admit a \(spatial\) Hamiltonian, and thus the methods of Hamiltonian Echo Backpropagation can be applied\. We end by deriving an explicit layer\-wise Hamiltonian recurrence relation governing inference for stationary solutions of both deep Fitzhugh\-Nagumo networks and deep Energy\-Based Models\.

## IIntroduction

The question of how biological brains coordinate their synaptic updates across the brain in order to perform learning is a major open question in neuroscience, and is also relevant to machine learning\. In machine learning, one typically uses stochastic gradient descent \(SGD\) to perform coordinated updates to all synapses with the goal of minimizing a loss or error function\. SGD has the benefit of being both conceptually simple and highly effective in training deep or hierarchical neural networks, and is even able to train spiking and neuromorphic architectures\(Leeet al\.,[2016](https://arxiv.org/html/2605.21568#bib.bib30); Wunderlich and Pehle,[2021](https://arxiv.org/html/2605.21568#bib.bib29)\)\.

However, in neuroscience, it is well\-known that the backpropagation algorithm, which is used to estimate the gradients that SGD requires to coordinate its weight updates, is not biologically plausible\(Grossberg,[1987](https://arxiv.org/html/2605.21568#bib.bib31)\)\. The primary reason for this is that backpropagation is a nonlocal algorithm: it requires an explicit "backward" graph along which to propagate gradient information, which must be matched to the forward graph at all times\. This has led some neuroscientists to conclude that SGD itself is not biologically plausible as a learning algorithm\. However, various alternatives to backpropagation have been proposed that are able to provide gradient estimates for biologically plausible synaptic updates\. A summary of recent work can be found inLillicrapet al\.\([2020](https://arxiv.org/html/2605.21568#bib.bib23)\)\.

### I\-ARelated Work

One promising alternative to backpropagation for training biologically plausible neural networks is Equilibrium Propagation \(EqProp\)\(Scellier and Bengio,[2017](https://arxiv.org/html/2605.21568#bib.bib1)\)\. EqProp applies to the class of networks known as Energy\-Based Models \(EBMs\), such as Hopfield networks and their modern variants\(Krotov and Hopfield,[2016](https://arxiv.org/html/2605.21568#bib.bib38)\)\.

Despite its restriction to EBMs, EqProp possesses several properties which make it attractive as a candidate for neural learning\. Perhaps most importantly, it is fully local in its weight updates, i\.e\. the update for a given synapse in the network depends only on the local activity difference of the neurons to which it is connected\. EqProp also possesses lower variance in the gradient estimate compared to node and weight perturbation and REINFORCE\-like algorithms, which makes it scalable to neural networks of practical relevance\(Høieret al\.,[2026](https://arxiv.org/html/2605.21568#bib.bib32); Scellieret al\.,[2023](https://arxiv.org/html/2605.21568#bib.bib34); Laborieuxet al\.,[2021](https://arxiv.org/html/2605.21568#bib.bib33)\)\. At the same time, it yields low\-bias gradient estimates compared to algorithms such as Feedback Alignment and Target Propagation\(Laborieuxet al\.,[2021](https://arxiv.org/html/2605.21568#bib.bib33)\)\. Lastly, the same network performs both inference and gradient estimation, so there is no need for a separate backward graph for propagating gradients\.

Another approach, with a long history of use in neuroscience, is Predictive Coding\(Friston,[2003](https://arxiv.org/html/2605.21568#bib.bib5); Bogacz,[2017](https://arxiv.org/html/2605.21568#bib.bib6); Millidgeet al\.,[2021](https://arxiv.org/html/2605.21568#bib.bib37)\)\. EqProp and Predictive Coding are closely related, as both can be formulated in terms of a bi\-level optimization problem: an inner optimization problem in which an energy function is minimized during inference, and an outer optimization problem where the weights are updated to perform learning\. As a result, Predictive Coding and EqProp share many theoretical and empirical properties, such as equivalence to backpropagation under certain conditions\(Millidgeet al\.,[2022](https://arxiv.org/html/2605.21568#bib.bib4)\), and the fact that both algorithms optimize a well\-defined objective function during learning\. It is worth noting that Predictive Coding, unlike EqProp, is capable of training feedforward networks\.

However, compared to Predictive Coding and many other alternatives to backpropagation, EqProp has the advantage of not requiring a separate, learned backward graph solely for propagating gradient information\. The mechanism by which EqProp achieves this is by exploiting the underlyingself\-adjointnessof Energy\-Based Models\. This self\-adjointness means these models, differentially, are equal to their own backward graph: their underlying activation Jacobians are symmetric\. This is a general property of so\-calledgradient systems, which follow dynamics determined by the gradient of a scalar function with respect to their state variables\. Note that this gradient following is different from gradient descent learning\. In gradient descent learning, it is thesynapse– i\.e\., weight – dynamics which follow the gradient of a loss function\. In contrast, in EBMs it is the neurons – i\.e\., activations – which follow the gradient of an energy function to perform inference\.

### I\-BGradient Systems and Equilibrium Propagation

Gradient systems, where the dynamics follow the gradient of an energy function, have the following interesting property: the gradient of anarbitraryscalar function defined on a subset of the system’s state variables can be propagated backwards, physically through the network, to any other subset of the state variables without requiring a separate backward graph\(Scellier and Bengio,[2017](https://arxiv.org/html/2605.21568#bib.bib1); Bengio and Fischer,[2015](https://arxiv.org/html/2605.21568#bib.bib2)\)\. Notably, this property enables one to perform local credit assignment in this class of networks, where the function defined on the state variables can be the loss or reward function of a neural network model\.111In other words, it shows that the loss gradient for any node in the network can be accurately calculated using only information that is locally available to that node, without requiring a separate backward graph containing copies of the transposed weightsWkTW\_\{k\}^\{T\}and activation function derivativesf′​\(h\)f^\{\\prime\}\(h\)\.

This property was originally proved by Scellier and Bengio\(Scellier and Bengio,[2017](https://arxiv.org/html/2605.21568#bib.bib1)\)in energy\-based models of neural networks, and was quickly shown to extend to many physical systems of interest such as optical networks\(Hugheset al\.,[2018](https://arxiv.org/html/2605.21568#bib.bib11)\), resistive electrical circuits\(Kendallet al\.,[2020](https://arxiv.org/html/2605.21568#bib.bib8)\), Ising machines\(Laydevantet al\.,[2024](https://arxiv.org/html/2605.21568#bib.bib9)\), mechanical mass\-spring oscillators\(Altmanet al\.,[2024](https://arxiv.org/html/2605.21568#bib.bib10)\), quantum circuits\(Wanjura and Marquardt,[2025](https://arxiv.org/html/2605.21568#bib.bib12); Scellier,[2024](https://arxiv.org/html/2605.21568#bib.bib13)\), and other systems of a physical origin\(Stern and Murugan,[2023](https://arxiv.org/html/2605.21568#bib.bib14)\)\.

More recently, Hamiltonian systems were also shown to support an analogous self\-adjoint algorithm for propagating gradient information, called Hamiltonian Echo Backpropagation \(HEB\)\(Lopez\-Pastor and Marquardt,[2023](https://arxiv.org/html/2605.21568#bib.bib15)\), which exploits the time\-reversibility of a Hamiltonian system to perform credit assignment across time\.

These efforts represent significant progress in solving the credit assignment problem in biophysical models of neural networks, and have led to new approaches for the design of neuro\-inspired computing architectures\(Yiet al\.,[2023](https://arxiv.org/html/2605.21568#bib.bib17); Martinet al\.,[2021](https://arxiv.org/html/2605.21568#bib.bib18)\)\.

From these properties, we believe methods which use the self\-adjoint structure of EBMs and Hamiltonian models to perform gradient estimation are promising candidates for explaining how gradient\-based learning could originate in simple biological networks\. However,there remains a gap between the classes of systems covered by EqProp and HEB, and real, biophysical models of neural circuits\. Notably,there is, to date, no existing method which is able to simultaneously incorporate nonlinear dissipation, gain, and time\-varying components\(all present in real biophysical neuronal networks\) into a single unified framework capable of performing credit assignment using self\-adjoint methods\.

### I\-CThe Fitzhugh\-Nagumo Model

Unlike gradient systems and conventional Hamiltonian systems, real biophysical neurons tend to simultaneously possess nonlinear dissipation, gain, and energy\-storage components such as capacitances\. This makes their treatment using techniques like Equilibrium Propagation and Hamiltonian Echo Backpropagation difficult as, in general, their dynamics are not self\-adjoint, and therefore cannot be given directly in terms of the gradient of a scalar energy function\.

As a canonical example of this, the Fitzhugh\-Nagumo model\(FitzHugh,[1961](https://arxiv.org/html/2605.21568#bib.bib19); Nagumoet al\.,[1962](https://arxiv.org/html/2605.21568#bib.bib20)\)is the simplest biophysical neuron model capable of generating action potentials\. It possesses these three categories of nonlinear dissipation, gain, and energy\-storage components, and as a result is capable of displaying a wide array of complex dynamics, including Turing patterns, spirals, traveling waves, solitons \(traveling spikes or standing pulse solutions\), and critical dynamics\(Cebrián\-Lacasaet al\.,[2024](https://arxiv.org/html/2605.21568#bib.bib21)\)\. The Fitzhugh\-Nagumo equations are given as:

d​ud​t\\displaystyle\\frac\{du\}\{dt\}=u−u3−v\\displaystyle=u\-u^\{3\}\-vd​vd​t\\displaystyle\\frac\{dv\}\{dt\}=ε​\(u−α​v−β\)\\displaystyle=\\varepsilon\(u\-\\alpha v\-\\beta\)
As one can see, the system of equations defined by the Fitzhugh\-Nagumo model is not self\-adjoint: its Jacobian is not symmetric with respect to the state variables\. However, it is part of a large class of \(active\) reaction\-diffusion systems which possess what has been termed "skew\-gradient" structure\(Yanagida,[2002b](https://arxiv.org/html/2605.21568#bib.bib7),[a](https://arxiv.org/html/2605.21568#bib.bib3)\)\. Skew\-gradient systems can be described as partitioning the system’s state variables into two parts, one describing an "activator" species, and the other describing an "inhibitor" species, with one species following the positive gradient and the other following the negative gradient of an energy function\. For example, in the Fitzhugh\-Nagumo model, the activator corresponds to the membrane potential of the neuron and the inhibitor corresponds to the recovery variable of the neuron\. Instead of being gradient systems, where all state variables follow the gradient of an energy function, they are mini\-maximizers with respect to the energy function\(Yanagida,[2002a](https://arxiv.org/html/2605.21568#bib.bib3)\)\. Some variables \(the activators\) seek to maximize the energy, while others \(the inhibitors\) seek to minimize the energy\.

In this work, we extend Equilibrium Propagation and Hamiltonian Echo Backpropagation to the stationary \(steady\-state\) solutions of such skew\-gradient systems, and analyze and train a skew\-gradient based deep neural network based on the Fitzhugh\-Nagumo model using Equilibrium Propagation\. By exploiting the underlying spatial Hamiltonian structure of the stationary states, we then derive a forward layer\-wise recursion which is able to perform inference in a single forward pass, given appropriate initial conditions, without relying on iterative temporal convergence to a fixed point, as is typical in the inference process of energy\-based models\. Finally, we apply these ideas to the more conventional energy\-based models ofScellier and Bengio \([2017](https://arxiv.org/html/2605.21568#bib.bib1)\); Bengio and Fischer \([2015](https://arxiv.org/html/2605.21568#bib.bib2)\)\.

## IIMain Result

A network of Fitzhugh\-Nagumo neurons with diffusive \(resistive\) coupling can be written in the following form, as described in Yanagida\(Yanagida,[2002b](https://arxiv.org/html/2605.21568#bib.bib7),[a](https://arxiv.org/html/2605.21568#bib.bib3)\)\.

T​ut=D​Δ​u\+f​\(u\)T\\mathrm\{u\}\_\{t\}=D\\Delta\\mathrm\{u\}\+\\mathrm\{f\(u\)\}\(1\)whereTTandDDare positive, diagonal matrices,Δ\\Deltais a Laplacian operator, which describes the neuron\-to\-neuron coupling, andf​\(u\)f\(u\)is a nonlinear function which can be expressed as the gradient of a scalar functionFF, usually called thefree energy, multiplied by a matrixQQ:

f​\(u\)=Q​∂F∂u,\\mathrm\{f\(u\)\}=Q\\frac\{\\partial F\}\{\\partial\\mathrm\{u\}\},\(2\)withQQsatisfyingQ2=𝕀2Q^\{2\}=\\mathbb\{I\}\_\{2\}\. For example, in our case, we chooseu=\[uv\]\\mathrm\{u\}=\\begin\{bmatrix\}u\\\\ v\\end\{bmatrix\}, andQ=\[100−1\]Q=\\begin\{bmatrix\}1&0\\\\ 0&\-1\\end\{bmatrix\}\. For notational convenience, we will refer to the uprightu\\mathrm\{u\}as the combined set of activator and inhibitor concentrations, and the scriptuuandvvfor the activator and inhibitor concentrations, respectively\.

Thus, the Fitzhugh\-Nagumo model with a single activator speciesuuand a single inhibitor speciesvvis given as follows:

\[τ100τ2\]​\[utvt\]=\[d100d2\]​\[Δ​uΔ​v\]\+\[100−1\]​\[∇uF∇vF\]\\begin\{bmatrix\}\\tau\_\{1\}&0\\\\ 0&\\tau\_\{2\}\\end\{bmatrix\}\\begin\{bmatrix\}u\_\{t\}\\\\ v\_\{t\}\\end\{bmatrix\}=\\begin\{bmatrix\}d\_\{1\}&0\\\\ 0&d\_\{2\}\\end\{bmatrix\}\\begin\{bmatrix\}\\Delta u\\\\ \\Delta v\\end\{bmatrix\}\+\\begin\{bmatrix\}1&0\\\\ 0&\-1\\end\{bmatrix\}\\begin\{bmatrix\}\\nabla\_\{u\}F\\\\ \\nabla\_\{v\}F\\end\{bmatrix\}
At a conceptual level, we will take the perspective of Nagumo and physically interpret the Fitzhugh\-Nagumo model as an electrical circuit with a tunnel diode nonlinearity and resistive coupling\. In the electrical circuit interpretation, we can think of the activator and inhibitor concentrationsuuandvvas the node voltages \(or loop currents\) of a circuit, as shown in Fig\. 1B\.

We will use the discrete Laplacian\(Muoloet al\.,[2024](https://arxiv.org/html/2605.21568#bib.bib22)\)in place of the usual continuous Laplacian, since we are primarily interested in discrete networks\. We can construct the Laplacian of a generic weighted, undirected \(diffusive\) graph using the incidence matrix of the graph, and explicitly write out the full set of discrete Fitzhugh\-Nagumo equations in terms of the graph Laplacian of the network\.

\[τ1​utτ2​vt\]=\[d1​𝕀n00d2​𝕀n\]​\[L1​uL2​v\]\+\[𝕀n00−𝕀n\]​\[∇uF∇vF\]\\begin\{bmatrix\}\\tau\_\{1\}u\_\{t\}\\\\ \\tau\_\{2\}v\_\{t\}\\end\{bmatrix\}=\\begin\{bmatrix\}d\_\{1\}\\mathbb\{I\}\_\{n\}&0\\\\ 0&d\_\{2\}\\mathbb\{I\}\_\{n\}\\end\{bmatrix\}\\begin\{bmatrix\}L\_\{1\}u\\\\ L\_\{2\}v\\end\{bmatrix\}\+\\begin\{bmatrix\}\\mathbb\{I\}\_\{n\}&0\\\\ 0&\-\\mathbb\{I\}\_\{n\}\\end\{bmatrix\}\\begin\{bmatrix\}\\nabla\_\{u\}F\\\\ \\nabla\_\{v\}F\\end\{bmatrix\}
Here,𝕀n\\mathbb\{I\}\_\{n\}is then×nn\\times nidentity matrix, anduuandvvarenn\-dimensional vectors for a network withnnnodes\. The LaplaciansL1L\_\{1\}andL2L\_\{2\}govern the network topology for the activator and inhibitor, respectively, and can be written in terms of the corresponding branch\-node incidence matrixBiB\_\{i\}and a diagonal matrix of positive conductances \(weights\) of the edgesYiY\_\{i\}, wherei=1i=1corresponds to the activator andi=2i=2the inhibitor\.

![Refer to caption](https://arxiv.org/html/2605.21568v1/Fitzhugh_Nagumo_Figure_1.png)Figure 1:A\) Time dynamics \(time increasing downwards\) of the Fitzhugh\-Nagumo model defined on a 1\-dimensional spatial interval\. The parameter settings of the model are tuned to display convergence to a stationary 1\-dimensional Turing pattern, i\.e\. spatial oscillations\. B\) Discrete 1D spatial model of coupled Fitzhugh\-Nagumo neurons \(activator only shown\), corresponding to a 1\-dimensional path graph\. Activator variableuuat positioniiin the path is denoteduiu^\{i\}\. Inhibitorviv^\{i\}is not explicitly shown, and takes the same form, but with no nonlinearity\. The nonlinear current through the tunnel diode is denoted byfif^\{i\}\. The momentum variablepip^\{i\}is the voltage drop between two adjacent nodesui\+1−uiu^\{i\+1\}\-u^\{i\}along the path\. C\) Coupled path graphs with coupling connections denoted bygj​kig^\{i\}\_\{jk\}, connecting nodejjof layeriito nodekkof layeri\+1i\+1, forming a deep residual network of FHN neurons\.Li=BiT​Yi​BiL\_\{i\}=B\_\{i\}^\{T\}Y\_\{i\}B\_\{i\}
As a consequence of this form, our Laplacian matrices are real\-symmetric\. Now, we can add an input current to the network\. We will represent this input current as a vector of the same dimension as the activatoruu, and therefore we can inject current into any or all of theuunodes in the network, whether they be input, hidden, or output nodes\. While we restrict our input current to the activator variablesuufor simplicity, the inhibitor variablesvvcan be treated in the same way\. Our system is now given by:

T​ut=D​L​u\+f​\(u\)−I,T\\mathrm\{u\}\_\{t\}=DL\\mathrm\{u\}\+\\mathrm\{f\(u\)\}\-I,
where again,u\\mathrm\{u\}is defined as\[u,v\]T\[u,v\]^\{T\},LLis a block diagonal matrix containingL1L\_\{1\}andL2L\_\{2\}, andI=\[Iu,0\]TI=\[I\_\{u\},0\]^\{T\}is the input current, assumed constant in time\. Note thatd1d\_\{1\}andd2d\_\{2\}inDDare global scaling factors for the conductances \(diffusion terms\) foruuandvv, multiplying their respective LaplaciansL1L\_\{1\}andL2L\_\{2\}\. Throughout the rest of this section, we restrict our attention to stationary solutions of \(6\), given by:

D​L​us\+f​\(us\)=IDL\\mathrm\{u\_\{s\}\}\+\\mathrm\{f\(u\_\{s\}\)\}=I
Here,us\\mathrm\{u\_\{s\}\}is a stationary solution of \(6\)\. We can analyze how small changes in our \(steady\-state\) input currentIIaffect the steady state solution by linearizing \(7\) aroundusu\_\{s\}\.

\(D​L\+∂f​\(us\)∂u\)​δ​u=δ​I\\left\(DL\+\\frac\{\\partial\\mathrm\{f\(u\_\{s\}\)\}\}\{\\partial\\mathrm\{u\}\}\\right\)\\delta\\mathrm\{u\}=\\delta I
We obtain:

\[d1​L1\+f′​\(u\)𝕀−𝕀d2​L2−γ\]​\[δ​uδ​v\]=\[δ​Iu0\],\\begin\{bmatrix\}d\_\{1\}L\_\{1\}\+f^\{\\prime\}\(u\)&\\mathbb\{I\}\\\\ \-\\mathbb\{I\}&d\_\{2\}L\_\{2\}\-\\gamma\\end\{bmatrix\}\\begin\{bmatrix\}\\delta u\\\\ \\delta v\\end\{bmatrix\}=\\begin\{bmatrix\}\\delta I\_\{u\}\\\\ 0\\end\{bmatrix\},
wheref′f^\{\\prime\}andγ\\gammaare given by the derivative off​\(u\)\\mathrm\{f\(u\)\}with respect touuandvv, respectively\. Note that, as a whole, the full linearized system matrix in \(9\) is not symmetric\. This means that it is not a gradient system, and is therefore not governed by a Lyapunov function whose value strictly decreases with time\. Instead, it is askew\-gradientsystem, where one set of variables aims to minimize, and the other set of variables aims to maximize, the free energy functionFF\(Yanagida,[2002a](https://arxiv.org/html/2605.21568#bib.bib3)\)\.

We can solve \(9\) for the activator variableuu, yielding the activator’s linearized response to the current perturbation at the steady stateus\\mathrm\{u\_\{s\}\}\. In particular, we find:

δ​u=\(\(d1​L1\+f′​\(us\)\)\+\(d2​L2−γ\)−1\)−1​δ​Iu=M−1​δ​Iu\\delta u=\\left\(\\left\(d\_\{1\}L\_\{1\}\+f^\{\\prime\}\(u\_\{s\}\)\\right\)\+\\left\(d\_\{2\}L\_\{2\}\-\\gamma\\right\)^\{\-1\}\\right\)^\{\-1\}\\delta I\_\{u\}=M^\{\-1\}\\delta I\_\{u\}
SinceL1L\_\{1\}andL2L\_\{2\}are symmetric by construction, andf′​\(u\)f^\{\\prime\}\(u\)andγ\\gammaare diagonal matrices, and since the inverse of a symmetric matrix is again symmetric, the matrixM−1M^\{\-1\}must therefore be symmetric\. Note thatM−1M^\{\-1\}is asubmatrixof the full system response matrix, corresponding to current injection of the activator speciesIuI\_\{u\}and subsequent change in the steady state of the activatorusu\_\{s\}\.

This result tells us that the network’seffective response matrixfor the activator variables \(the Jacobian of the activator variablesusu\_\{s\}with respect to the injected currentIuI\_\{u\}\) is self\-adjoint for all nodes, and across all settings of the Fitzhugh\-Nagumo parameter values, as long as the network is at a steady state\. Thus, the methods of Equilibrium Propagation can be applied\. A summary of why EqProp\-like learning rules follow directly from the self\-adjointness of the underlying network is now given\.

Let us assume we wish to compute the gradient of a loss functionℒ\\mathcal\{L\}with respect to a parameter at layerii, saygikg\_\{i\}^\{k\}\. The lossℒ\\mathcal\{L\}is a function of the network’s output variables\. Let the output variables be denoteduyu\_\{y\}\. These output variables are in turn a function of the network’s hidden state variables at layerii, let these be denoteduhiu\_\{h\_\{i\}\}\. Finally, the hidden state at layeriiis a function of the parameter valuegikg\_\{i\}^\{k\}\. Then we have

d​ℒ=∂ℒ∂uy​∂uy∂uhi​∂uhi∂gik​d​gikd\\mathcal\{L\}=\\frac\{\\partial\\mathcal\{L\}\}\{\\partial u\_\{y\}\}\\frac\{\\partial u\_\{y\}\}\{\\partial u\_\{h\_\{i\}\}\}\\frac\{\\partial u\_\{h\_\{i\}\}\}\{\\partial g\_\{i\}^\{k\}\}dg\_\{i\}^\{k\}
This equation describes the differential change in the value of the lossℒ\\mathcal\{L\}via a change in the parameter of interest,gikg\_\{i\}^\{k\}\. In other words, this is the forward directional derivative of the loss with respect to a change in the parameter of interest\. Since there are many different parameters in the network, and only a single scalar loss function, the adjoint of this equation is typically used, which is also known as backward\-mode differentiation\. This is given by transposing the above relation to yield:

∂ℒ∂gik=\(∂uhi∂gik\)T​\(∂uy∂uhi\)T​\(∂ℒ∂uy\)T\\frac\{\\partial\\mathcal\{L\}\}\{\\partial g\_\{i\}^\{k\}\}=\\left\(\\frac\{\\partial u\_\{h\_\{i\}\}\}\{\\partial g\_\{i\}^\{k\}\}\\right\)^\{T\}\\left\(\\frac\{\\partial u\_\{y\}\}\{\\partial u\_\{h\_\{i\}\}\}\\right\)^\{T\}\\left\(\\frac\{\\partial\\mathcal\{L\}\}\{\\partial u\_\{y\}\}\\right\)^\{T\}\(3\)
We see that in backward mode differentiation, the network’s forward Jacobian∂uy∂uhi\\frac\{\\partial u\_\{y\}\}\{\\partial u\_\{h\_\{i\}\}\}is replaced by its transpose\. For a network to be self\-adjoint, this matrix must be symmetric, and therefore the transposed quantity is identical to the network’s forward Jacobian\. This means that we can performa single inference\-time perturbationto the output neurons, proportional to∂ℒ∂uy\\frac\{\\partial\\mathcal\{L\}\}\{\\partial u\_\{y\}\}, which travels backwards through∂uy∂uhi=\(∂uy∂uhi\)T\\frac\{\\partial u\_\{y\}\}\{\\partial u\_\{h\_\{i\}\}\}=\\left\(\\frac\{\\partial u\_\{y\}\}\{\\partial u\_\{h\_\{i\}\}\}\\right\)^\{T\}, carrying gradient information to the rest of the nodes in the network\. These perturbations can be locally measured at each hidden layer and used to assign credit to the parameters via Eq\. \(3\)\. The general mechanism of encoding gradient information in terms of perturbations or activity differences is known as NGRAD \(Neural Gradient Representation by Activity Differences\), and is discussed in general inLillicrapet al\.\([2020](https://arxiv.org/html/2605.21568#bib.bib23)\)\.

Because there is only a single perturbation required for the entire gradient estimation step, EqProp and related methods have much lower variance than node or weight perturbation methods, where the variance in the gradient estimate scales with the number of parameters or nodes in the network\(Lillicrapet al\.,[2020](https://arxiv.org/html/2605.21568#bib.bib23)\)\.

## IIITraining Deep Fitzhugh\-Nagumo Networks

In order to demonstrate that the above theoretical analysis holds experimentally, we demonstrate the training of deep networks of Fitzhugh\-Nagumo neurons, in an analogous fashion to deep EBMs, using Equilibrium Propagation\. We mirror the experiments done inScellier and Bengio \([2017](https://arxiv.org/html/2605.21568#bib.bib1)\)by training a deep Fitzhugh\-Nagumo network with 5 hidden layers on MNIST using EqProp\. Code is available on Github:[Self\-Adjoint\-Learning](https://github.com/jackdkendall/Self-Adjoint-Learning)

For the deep FHN networks we train, we initialize the FHN parameter valuesδ\\delta,ε\\varepsilon,α\\alpha,β\\betain the Turing pattern forming regime\. We keep these parameters fixed during training\. We initialize the weights normally around zero, as in a standard EBM\. Since the skew\-gradient formulation of the FHN nonlinearity already admits negative conductances, we do not clamp parameter conductances to be strictly positive, and we simply treat negative conductances as effective gain elements\. For Equilibrium Propagation, we find that a large value ofβn​u​d​g​e\\beta\_\{nudge\}, near 0\.9, works best\. We use the centered difference formulation of EqProp, and consistent with prior work on EqProp training, ensure that earlier layers have larger learning rates than later layers\.

TABLE I:Optimal hyperparameters and MNIST test error for a 5\-layer Fitzhugh\-Nagumo Network trained with Equilibrium PropagationFor inference, we follow the time\-dynamics defined by the FHN model to approximate convergence, which typically takes 50 iterations, and we use a fixed step size of 0\.1\. Our weight initialization scale \(scalar multiplier of standard normal initialization\) was chosen to be 0\.01, with 0\.014 showing the best results\. We noticed in some instances, the occurrence of loss spikes during training\. We speculate that these may be caused by the FHN model leaving the steady\-state dynamical regime it was initialized in, via accumulation of weight updates which increase the effective gain of the network\. More experiments are needed to investigate this behavior\.

## IVHamiltonian Inference

While Energy\-Based Models have a number of desirable properties, such as good sample diversity in generative models and natural compositionality properties\(Duet al\.,[2020](https://arxiv.org/html/2605.21568#bib.bib35)\), they have not seen widespread use in industry, because of their expensive inference process relative to feedforward networks\. Whereas feedforward networks require a single forward pass through the network to perform inference, EBMs require many forward passes to reach convergence to the energy minimum\.

We will show that deep EBMs possess a feedforward Hamiltonian formulation,which is equivalent to the original EBM up to specification of boundary conditions\. This Hamiltonian formulation is entirely feedforward in the inference process, but requires the specification of an additional "momentum" variable at the input to the network, as it turns the original EBM’s two\-point boundary problem into a purely initial value problem\. To begin with, we will generalize the Hamiltonian formulation of the FHN from the continuous setting, where it is well established\(Parra\-Rivaset al\.,[2025](https://arxiv.org/html/2605.21568#bib.bib24)\), to the discrete network setting\. Then, we will transfer these results to the Energy\-Based Models considered by Scellier and Bengio\.

### IV\-AHamiltonian Formulation of the Fitzhugh\-Nagumo Model

The fact that the Fitzhugh\-Nagumo model admits a Hamiltonian formulation was first noticed by Yanagida\(Yanagida,[2002b](https://arxiv.org/html/2605.21568#bib.bib7),[a](https://arxiv.org/html/2605.21568#bib.bib3)\)\. There is now a large body of literature which characterizes the stationary solutions of the Fitzhugh\-Nagumo model in terms ofspatial dynamics, i\.e\. they describe the spatial evolution of the solutionsus​\(x\)\\mathrm\{u\_\{s\}\}\(x\)in terms of a Hamiltonian which is conserved across a spatial dimensionxx, wherexxnow plays the role of "time" in a typical nonlinear Hamiltonian problem\(Kuwamura,[2004](https://arxiv.org/html/2605.21568#bib.bib25); Karasözenet al\.,[2017](https://arxiv.org/html/2605.21568#bib.bib26); Parra\-Rivaset al\.,[2025](https://arxiv.org/html/2605.21568#bib.bib24)\)\. However, these results hold for continuous 1D systems, and to the author’s knowledge, there have been no attempts at generalizing this framework to the setting of discrete graphs\. So, we will now apply the Hamiltonian formalism of the Fitzhugh\-Nagumo model to the discrete undirected graph Laplacian setting\. We will analyze the autonomous case, whereI=0I=0\.

In\(Kuwamura,[2004](https://arxiv.org/html/2605.21568#bib.bib25)\), the 1\-dimensional continuous Laplacian is considered, with the spatial coordinate given byxx\. Steady\-state solutions of the continuous 1\-dimensional Fitzhugh\-Nagumo model are described by:

D​∂2∂x2​u=−f​\(u\)D\\frac\{\\partial^\{2\}\}\{\\partial x^\{2\}\}\\mathrm\{u\}=\-\\mathrm\{f\(u\)\}
whereDDandf​\(u\)\\mathrm\{f\(u\)\}obey the constraints given in \(1\)\. For our analysis, we will replace the 1\-dimensional continuous Laplacian∂2∂x2\\frac\{\\partial^\{2\}\}\{\\partial x^\{2\}\}with a generic weighted, undirected graph LaplacianLL\. We will assume the Laplacian has been grounded, or reduced, such that its zero eigenvalue is removed\. SinceLLis positive definite, it possesses a unique non\-negative square root, which we will denoteL1/2L^\{1/2\}, such thatL1/2​\(L1/2\)T=L1/2​L1/2=LL^\{1/2\}\(L^\{1/2\}\)^\{T\}=L^\{1/2\}L^\{1/2\}=L\. We will define the "spatial velocity"eeto be the square root of the Laplacian acting onu\\mathrm\{u\}\. This will be our discrete analog ofux\\mathrm\{u\}\_\{x\}, the spatial gradient ofu\\mathrm\{u\}alongxx\.

e:=L1/2​u\\mathrm\{e\}:=L^\{1/2\}\\mathrm\{u\}
Now, followingKuwamura \([2004](https://arxiv.org/html/2605.21568#bib.bib25)\), we define a new variableZZas:

Z=\[ue\]Z=\\begin\{bmatrix\}\\mathrm\{u\}\\\\ \\mathrm\{e\}\\end\{bmatrix\}
We then define the spatial gradient ofZZas:

Zx=\[L1/2​uL1/2​e\]=\[eL​u\]Z\_\{x\}=\\begin\{bmatrix\}L^\{1/2\}\\mathrm\{u\}\\\\ L^\{1/2\}\\mathrm\{e\}\\end\{bmatrix\}=\\begin\{bmatrix\}\\mathrm\{e\}\\\\ L\\mathrm\{u\}\\end\{bmatrix\}
Using this notation, we can break the set of n second\-order difference equations into a set of 2n first\-order difference equations:

ux\\displaystyle\\mathrm\{u\}\_\{x\}=e=L1/2​u\\displaystyle=\\mathrm\{e\}=L^\{1/2\}\\mathrm\{u\}D​ex\\displaystyle D\\mathrm\{e\}\_\{x\}=D​L1/2​e=−f​\(u\)\\displaystyle=DL^\{1/2\}\\mathrm\{e\}=\-\\mathrm\{f\(u\)\}
Indeed, we see thatD​L1/2​\(L1/2​u\)=D​L​uDL^\{1/2\}\(L^\{1/2\}\\mathrm\{u\}\)=DL\\mathrm\{u\}\. Next, define the following \(spatial\) Hamiltonian:

H:=12​eT​D​Q​e\+F​\(u\)H:=\\frac\{1\}\{2\}\\mathrm\{e\}^\{T\}DQe\+F\(\\mathrm\{u\}\)
Where

F:=∫0uf​\(u′\)​𝑑u′=12​α​u2−14​α​u4−u​v\+12​β​v2F:=\\int\_\{0\}^\{\\mathrm\{u\}\}\\mathrm\{f\(u^\{\\prime\}\)\}d\\mathrm\{u^\{\\prime\}\}=\\frac\{1\}\{2\}\\alpha u^\{2\}\-\\frac\{1\}\{4\}\\alpha u^\{4\}\-uv\+\\frac\{1\}\{2\}\\beta v^\{2\}
This corresponds to the Fitzhugh\-Nagumo model Free Energy function given inYanagida \([2002b](https://arxiv.org/html/2605.21568#bib.bib7)\)\. Finally, we can write the Fitzhugh\-Nagumo equations in their Hamiltonian, or symplectic, form:

\[0Q​D−Q​D0\]​\[uxex\]=−\[∂H∂u∂H∂e\]=−∂H∂Z\\begin\{bmatrix\}0&QD\\\\ \-QD&0\\end\{bmatrix\}\\begin\{bmatrix\}\\mathrm\{u\}\_\{x\}\\\\ \\mathrm\{e\}\_\{x\}\\end\{bmatrix\}=\-\\begin\{bmatrix\}\\frac\{\\partial H\}\{\\partial\\mathrm\{u\}\}\\\\ \\frac\{\\partial H\}\{\\partial\\mathrm\{e\}\}\\end\{bmatrix\}=\-\\frac\{\\partial H\}\{\\partial Z\}
Or, more compactly, as:

K​Zx=−∂H∂ZKZ\_\{x\}=\-\\frac\{\\partial H\}\{\\partial Z\}\(4\)
The existence of a Hamiltonian formulation in the coordinatesZZhas several important implications for discrete 1\-dimensional Laplacians\. First, it means if one knows the pair\(u,e\)\(\\mathrm\{u,e\}\)at any nodeiiin the path graph, the values of\(u,e\)\(\\mathrm\{u,e\}\)on the entire graph can be constructed by recursively solving Eq\. \(4\) forZxZ\_\{x\}, starting fromii, and integrating\(u\+ux,e\+ex\)\(u\+u\_\{x\},e\+e\_\{x\}\)along the graph, rather than iteratively simulating the dynamics of the entire network to convergence, as is commonly done in inference of energy\-based models\. Second, the spatial HamiltonianHHis constant \(for stationary solutions\) along the path graph\. In other words, the sum of all currents into each node of the graph is equal to zero\. Third, the spatial evolution of\(u,e\)\(\\mathrm\{u,e\}\)along the graph is formally equivalent to a conservative nonlinear oscillator problem in discrete time, for Turing pattern steady\-state solutions\.

We will now use these ideas to develop a Hamiltonian approach for an example of a Fitzhugh\-Nagumo network with the structure of a deep residual neural network\. We will start with a discrete 1\-dimensional \(path graph\) LaplacianLLwithnnnodes given by222Note that we treat the boundary of the Laplacian as if it were ”impedance matched”, i\.e\. a ghost node is added to the boundary nodes which double the value of the Laplacian’s diagonal terms at the first and last row, making the diagonal constant\.:

L=\[2−10000−12−10000−12⋱0000⋱⋱−10000−12−10000−12\]L=\\begin\{bmatrix\}2&\-1&0&0&0&0\\\\ \-1&2&\-1&0&0&0\\\\ 0&\-1&2&\\ddots&0&0\\\\ 0&0&\\ddots&\\ddots&\-1&0\\\\ 0&0&0&\-1&2&\-1\\\\ 0&0&0&0&\-1&2\\end\{bmatrix\}
We will consider the result of coupling many of these 1D path graphs together in parallel\. To simplify the analysis in this section, we will use the following formulation of the Fitzhugh\-Nagumo model, used inParra\-Rivaset al\.\([2025](https://arxiv.org/html/2605.21568#bib.bib24)\)to study the spatial dynamics of the steady states\.

δ2​Δ​u\+u−u3−v\\displaystyle\\delta^\{2\}\\Delta u\+u\-u^\{3\}\-v=0\\displaystyle=0Δ​v\+ε​\(u−α​v−β\)\\displaystyle\\Delta v\+\\varepsilon\(u\-\\alpha v\-\\beta\)=0\\displaystyle=0
First, we can extend our 1\-dimensional network by expanding our path graph Laplacian to includemmseparate path graphs, initially uncoupled from one another, each with path lengthnn\. Then, we can introduce layer\-wise coupling between adjacent sets of nodes in the independent path graphs, i\.e\. by creating a Laplacian with the structure of a deep residual neural network, as shown in Fig\. 1C\. For stationary solutions corresponding to Turing patterns, the resulting system is formally equivalent to acoupledsystem of nonlinear conservative oscillators in discrete time\.

We will now derive a layer\-wise recursion for steady\-state solutions of this network\. This means that instead of simulating the \(skew\-gradient\) time dynamics of the network until convergence, which requires many forward passes over the network until the steady state is found, we can directly compute the steady state solution of each layer directly from the steady state of the previous layer\. That is, we can perform inference in a single forward pass, given an appropriate initial condition, which we will define as the initial position and momentum\(u0,p0\)\(u^\{0\},p^\{0\}\)\. Thus, instead of defining our network only in terms of theui\\mathrm\{u\}^\{i\}variables, we must describe it in terms of the phase\-space variables\(ui,pi\)\(\\mathrm\{u\}^\{i\},\\mathrm\{p\}^\{i\}\)\.

Next, since the network is at a steady state, we note that the total current through the capacitors is equal to zero\. We define the value of the activator at nodeiialong pathkkto beukiu\_\{k\}^\{i\}\. The current through the nonlinear tunnel diode branchfkif^\{i\}\_\{k\}is equal to the gradient of the free energyFFwith respect toukiu\_\{k\}^\{i\}\.

The momentum variable, we define to be the voltage drop between adjacent nodes in the same path,pi:=ui\+1−uip^\{i\}:=u^\{i\+1\}\-u^\{i\}\. Here,uiu^\{i\}andpip^\{i\}are vectors of node activator variables and edge drops, respectively, of all nodes at layeriiin their path graph\. The coupling matrices couple the nodes at layeriito the nodes at layeri\+1i\+1, and are given bygj​kig^\{i\}\_\{jk\}\.

Now, the condition of the spatial Hamiltonian being constant along the path in the abstract Laplacian analysis given in Section 3\.1 is equivalent to the statement that for all nodes in the network, the currents into that node sum to zero\. Or in other words, it is a consequence of Kirchhoff’s current law\. This can also be viewed as a homology relation on the network\(Smale,[1971](https://arxiv.org/html/2605.21568#bib.bib27)\)\.

For uncoupled path graphs, one path of which is shown in Fig\. 1B, the total current at nodei\+1i\+1in pathkkof the graph is given by:

δ2​pki\+1−δ2​pki−fki\+1=0\\delta^\{2\}p^\{i\+1\}\_\{k\}\-\\delta^\{2\}p^\{i\}\_\{k\}\-f^\{i\+1\}\_\{k\}=0
From this relation, the voltage drop to the next node in the path graph,pki\+1p^\{i\+1\}\_\{k\}, can be computed fromfki\+1f^\{i\+1\}\_\{k\}andpkip^\{i\}\_\{k\}\.

For the case of coupled paths, the total current is more complicated, but illuminates the connection between the iterative energy\-based dynamics and the direct Hamiltonian computation of the steady\-state solution\. The total current at nodei\+1i\+1for coupled paths is given by:

δ2​pki\+1−δ2​pki−fki\+1\+∑jm\(uji−uki\+1\)​gk​ji\+∑j\(uji\+2−uki\+1\)​gj​ki\+1=0\\delta^\{2\}p^\{i\+1\}\_\{k\}\-\\delta^\{2\}p^\{i\}\_\{k\}\-f^\{i\+1\}\_\{k\}\+\\sum\_\{j\}^\{m\}\\left\(u\_\{j\}^\{i\}\-u\_\{k\}^\{i\+1\}\\right\)g\_\{kj\}^\{i\}\+\\sum\_\{j\}\\left\(u^\{i\+2\}\_\{j\}\-u^\{i\+1\}\_\{k\}\\right\)g\_\{jk\}^\{i\+1\}=0\(5\)
Using the definitions ofpkip^\{i\}\_\{k\}andpki\+1p^\{i\+1\}\_\{k\}, we can rewrite the two sums above as:

∑jm\(uji−uki\+1\)​gk​ji=∑jm\(\(uji\+1−pji\)−uki\+1\)​gk​ji\\sum\_\{j\}^\{m\}\\left\(u\_\{j\}^\{i\}\-u\_\{k\}^\{i\+1\}\\right\)g\_\{kj\}^\{i\}=\\sum\_\{j\}^\{m\}\\left\(\(u^\{i\+1\}\_\{j\}\-p^\{i\}\_\{j\}\)\-u\_\{k\}^\{i\+1\}\\right\)g\_\{kj\}^\{i\}
∑j\(uji\+2−uki\+1\)​gj​ki\+1=∑j\(\(uji\+1\+pji\+1\)−uki\+1\)​gj​ki\+1\\sum\_\{j\}\\left\(u^\{i\+2\}\_\{j\}\-u^\{i\+1\}\_\{k\}\\right\)g\_\{jk\}^\{i\+1\}=\\sum\_\{j\}\\left\(\(u^\{i\+1\}\_\{j\}\+p^\{i\+1\}\_\{j\}\)\-u^\{i\+1\}\_\{k\}\\right\)g\_\{jk\}^\{i\+1\}
This reduces Eq\. \(5\) to the following equation, which is only a function ofpip^\{i\},pi\+1p^\{i\+1\}, andui\+1u^\{i\+1\}, and where\(Gi\)j​k:=gj​ki\(G^\{i\}\)\_\{jk\}:=g\_\{jk\}^\{i\},g^i\\hat\{g\}^\{i\}is the sum along the columns ofGiG^\{i\}, andg~i\\tilde\{g\}^\{i\}is the sum along the rows ofGiG^\{i\}\.

\(Gi\+1\+δ2​𝕀m\)​pi\+1=\(\(Gi\)T\+δ2​𝕀m\)​pi\+fi\+1−\(\(Gi\)T\+Gi\+1−d​i​a​g​\(g^i\+g~i\+1\)\)​ui\+1\\left\(G^\{i\+1\}\+\\delta^\{2\}\\mathbb\{I\}\_\{m\}\\right\)p^\{i\+1\}=\\left\(\(G^\{i\}\)^\{T\}\+\\delta^\{2\}\\mathbb\{I\}\_\{m\}\\right\)p^\{i\}\+f^\{i\+1\}\-\\left\(\(G^\{i\}\)^\{T\}\+G^\{i\+1\}\-diag\(\\hat\{g\}^\{i\}\+\\tilde\{g\}^\{i\+1\}\)\\right\)u^\{i\+1\}
LetMi=\(Gi\+1\+δ2​𝕀m\)−1M^\{i\}=\(G^\{i\+1\}\+\\delta^\{2\}\\mathbb\{I\}\_\{m\}\)^\{\-1\},Ni=\(\(Gi\)T\+δ2​𝕀m\)N^\{i\}=\(\(G^\{i\}\)^\{T\}\+\\delta^\{2\}\\mathbb\{I\}\_\{m\}\), andOi=\(\(Gi\)T\+Gi\+1−d​i​a​g​\(g^i\+g~i\+1\)\)O^\{i\}=\(\(G^\{i\}\)^\{T\}\+G^\{i\+1\}\-diag\(\\hat\{g\}^\{i\}\+\\tilde\{g\}^\{i\+1\}\)\)\. Then we can solve forpi\+1p^\{i\+1\}as

pi\+1=Mi​Ni​pi\+Mi​fi\+1−Mi​Oi​ui\+1p^\{i\+1\}=M^\{i\}N^\{i\}p^\{i\}\+M^\{i\}f^\{i\+1\}\-M^\{i\}O^\{i\}u^\{i\+1\}
Given auiu\_\{i\}and apip\_\{i\}at any layerii, we can then compute the steady\-state solution of the network at layeri\+1i\+1via the following layer\-wise recursion relation\. This yields the final Hamiltonian recurrence relations for the deep residual Fitzhugh\-Nagumo network for the activator variables:

ui\+1\\displaystyle u^\{i\+1\}=ui\+pi\\displaystyle=u^\{i\}\+p^\{i\}pi\+1\\displaystyle p^\{i\+1\}=Mi​Ni​pi\+Mi​fi\+1−Mi​Oi​ui\+1\\displaystyle=M^\{i\}N^\{i\}p^\{i\}\+M^\{i\}f^\{i\+1\}\-M^\{i\}O^\{i\}u^\{i\+1\}
And for the inhibitor variables, which can be treated similarly, but with no inter\-path coupling:

vi\+1\\displaystyle v^\{i\+1\}=vi\+qi\\displaystyle=v^\{i\}\+q^\{i\}qi\+1\\displaystyle q^\{i\+1\}=ε​\(ui\+1−α​vi\+1−β\)\\displaystyle=\\varepsilon\(u^\{i\+1\}\-\\alpha v^\{i\+1\}\-\\beta\)

### IV\-BHamiltonian Formulation of Deep EBMs

Consider the energy function defined inScellier and Bengio \([2017](https://arxiv.org/html/2605.21568#bib.bib1)\), which is a kind of modified Hopfield energy, originally used inBengio and Fischer \([2015](https://arxiv.org/html/2605.21568#bib.bib2)\)\.

E=12​∑iui2−12​∑i≠jWi​j​ρ​\(ui\)​ρ​\(uj\)−∑ibi​ρ​\(ui\)E=\\frac\{1\}\{2\}\\sum\_\{i\}u\_\{i\}^\{2\}\-\\frac\{1\}\{2\}\\sum\_\{i\\neq j\}W\_\{ij\}\\rho\(u\_\{i\}\)\\rho\(u\_\{j\}\)\-\\sum\_\{i\}b\_\{i\}\\rho\(u\_\{i\}\)
The gradient of the energy function with respect touiu\_\{i\}, which drives the dynamics of the network, is then given by:

∂E∂ui=ui−ρ′​\(ui\)​\(∑jWi​j​ρ​\(uj\)−bi\)\\frac\{\\partial E\}\{\\partial u\_\{i\}\}=u\_\{i\}\-\\rho^\{\\prime\}\(u\_\{i\}\)\\left\(\\sum\_\{j\}W\_\{ij\}\\rho\(u\_\{j\}\)\-b\_\{i\}\\right\)
We will now show how to perform the above Hamiltonian analysis for this energy\-based model\. First, we will perform a change of variables which will allow us to transfer the nonlinearity from the interaction term involvingWi​jW\_\{ij\}to the local term involvinguiu\_\{i\}\. Assuming the nonlinear activation functionρ​\(ui\)\\rho\(u\_\{i\}\)is strictly monotonic \(e\.g\. sigmoidal\), we can make the following change of variables:

vi:=ρ​\(ui\)⟹ui=ρ−1​\(vi\)v\_\{i\}:=\\rho\(u\_\{i\}\)\\implies u\_\{i\}=\\rho^\{\-1\}\(v\_\{i\}\)
The gradient of the energy function with respect touiu\_\{i\}can be written in terms ofviv\_\{i\}using the chain rule as:

∂E∂ui=∂E∂vi​∂vi∂ui=∂E∂vi​ρ′​\(ui\)\\frac\{\\partial E\}\{\\partial u\_\{i\}\}=\\frac\{\\partial E\}\{\\partial v\_\{i\}\}\\frac\{\\partial v\_\{i\}\}\{\\partial u\_\{i\}\}=\\frac\{\\partial E\}\{\\partial v\_\{i\}\}\\rho^\{\\prime\}\(u\_\{i\}\)
Substituting this into the LHS of \(43\) and the definition ofviv\_\{i\}into the RHS of \(43\):

∂E∂vi​ρ′​\(ui\)=ρ−1​\(vi\)−ρ′​\(ui\)​\(∑jWi​j​vj−bi\)\\frac\{\\partial E\}\{\\partial v\_\{i\}\}\\rho^\{\\prime\}\(u\_\{i\}\)=\\rho^\{\-1\}\(v\_\{i\}\)\-\\rho^\{\\prime\}\(u\_\{i\}\)\\left\(\\sum\_\{j\}W\_\{ij\}v\_\{j\}\-b\_\{i\}\\right\)
∂E∂vi=ρ−1​\(vi\)ρ′​\(ρ−1​\(vi\)\)−∑jWi​j​vj−bi\\frac\{\\partial E\}\{\\partial v\_\{i\}\}=\\frac\{\\rho^\{\-1\}\(v\_\{i\}\)\}\{\\rho^\{\\prime\}\(\\rho^\{\-1\}\(v\_\{i\}\)\)\}\-\\sum\_\{j\}W\_\{ij\}v\_\{j\}\-b\_\{i\}
We can absorb the nonlinear termρ−1​\(vi\)ρ′​\(ρ−1​\(vi\)\)\\frac\{\\rho^\{\-1\}\(v\_\{i\}\)\}\{\\rho^\{\\prime\}\(\\rho^\{\-1\}\(v\_\{i\}\)\)\}and the biasbib\_\{i\}into a single local nonlinearityfif\_\{i\}atviv\_\{i\}:

∂E∂vi=fi​\(vi\)−∑jWi​j​vj\\frac\{\\partial E\}\{\\partial v\_\{i\}\}=f\_\{i\}\(v\_\{i\}\)\-\\sum\_\{j\}W\_\{ij\}v\_\{j\}
We will now make the assumption that the network is structured as adeepenergy\-based model, that is, that the weights and activations can be partitioned into layers, and we will assume a constant layer width\. The gradient of the energy with respect to theiith neuron at layerl\+1l\+1in the network is then given by:

∂E∂vil\+1=fi​\(vil\+1\)−∑jWj​il​vjl−∑jWi​jl\+1​vjl\+2\\frac\{\\partial E\}\{\\partial v\_\{i\}^\{l\+1\}\}=f\_\{i\}\(v\_\{i\}^\{l\+1\}\)\-\\sum\_\{j\}W\_\{ji\}^\{l\}v\_\{j\}^\{l\}\-\\sum\_\{j\}W\_\{ij\}^\{l\+1\}v\_\{j\}^\{l\+2\}
In layerwise vector form, this can be written as:

∂E∂vl\+1=fi​\(vl\+1\)−\(Wl\)T​vl−Wl\+1​vl\+2\\frac\{\\partial E\}\{\\partial v^\{l\+1\}\}=f\_\{i\}\(v^\{l\+1\}\)\-\(W^\{l\}\)^\{T\}v^\{l\}\-W^\{l\+1\}v^\{l\+2\}
As in Section 3\.2, we will make the following definition:

pl:=vl\+1−vl⟹vl\+1=vl\+plp^\{l\}:=v^\{l\+1\}\-v^\{l\}\\implies v^\{l\+1\}=v^\{l\}\+p^\{l\}
This will allow us to express the termsvlv^\{l\}andvl\+2v^\{l\+2\}in terms ofvl\+1v^\{l\+1\},plp^\{l\}, andpl\+1p^\{l\+1\}:

∂E∂vl\+1=f​\(vl\+1\)−\(Wl\)T​\(vl\+1−pl\)−Wl\+1​\(vl\+1\+pl\+1\)\\frac\{\\partial E\}\{\\partial v^\{l\+1\}\}=f\(v^\{l\+1\}\)\-\(W^\{l\}\)^\{T\}\(v^\{l\+1\}\-p^\{l\}\)\-W^\{l\+1\}\(v^\{l\+1\}\+p^\{l\+1\}\)
For stationary solutions of the network, the above energy gradient sums to zero\. Rearranging terms, we get:

0=f​\(vl\+1\)−Wl\+1​pl\+1\+\(Wl\)T​pl−\(\(Wl\)T\+Wl\+1\)​vl\+10=f\(v^\{l\+1\}\)\-W^\{l\+1\}p^\{l\+1\}\+\(W^\{l\}\)^\{T\}p^\{l\}\-\(\(W^\{l\}\)^\{T\}\+W^\{l\+1\}\)v^\{l\+1\}
Wl\+1​pl\+1=f​\(vl\+1\)\+\(Wl\)T​pl−\(\(Wl\)T\+Wl\+1\)​vl\+1W^\{l\+1\}p^\{l\+1\}=f\(v^\{l\+1\}\)\+\(W^\{l\}\)^\{T\}p^\{l\}\-\(\(W^\{l\}\)^\{T\}\+W^\{l\+1\}\)v^\{l\+1\}
This allows us to directly obtain the following Hamiltonian recurrence relations onvlv^\{l\}andplp^\{l\}, which gives the equation for exact layer\-wise inference of steady\-state solutions in deep energy\-based models:

vl\+1\\displaystyle v^\{l\+1\}=vl\+pl\\displaystyle=v^\{l\}\+p^\{l\}pl\+1\\displaystyle p^\{l\+1\}=Ml​f​\(vl\+1\)\+Nl​pl−\(Nl\+𝕀\)​vl\+1\\displaystyle=M^\{l\}f\(v^\{l\+1\}\)\+N^\{l\}p^\{l\}\-\(N^\{l\}\+\\mathbb\{I\}\)v^\{l\+1\}
Where

Ml\\displaystyle M^\{l\}=\(Wl\+1\)−1\\displaystyle=\(W^\{l\+1\}\)^\{\-1\}Nl\\displaystyle N^\{l\}=\(Wl\+1\)−1​\(Wl\)T\\displaystyle=\(W^\{l\+1\}\)^\{\-1\}\(W^\{l\}\)^\{T\}
We remark that the analogy here is that for energy\-based models at their steady\-state solution, the state variables are stationary with respect to time\. As a result, the gradients of the energy with respect to each of the state variables \(which can be interpreted as "currents", if the state variables are interpreted as "voltages"\), satisfy a local conservation law: the components of the gradient sum to zero for each node in the network\. These additional conservation equations on the gradients can then be used to calculate the components of the gradient along the forward integration direction\. Given that the state variable and the momentum variable at one layer are known, the next layer’s activation can be calculated from this component of the gradient\.

## VHamiltonian Simulations

In this section, we focus on simulations of the Hamiltonian inference process for the deep Fitzhugh\-Nagumo networks described above\. For the Hamiltonian formulation of Deep EBMs, we defer this discussion to a separate paper which focuses on this topic\.

![Refer to caption](https://arxiv.org/html/2605.21568v1/Energy_Based_vs_Hamiltonian_Figure_2.png)Figure 2:Comparison between the result of a time\-dynamics simulation and a Hamiltonian spatial integration in depth for a deep Fitzhugh\-Nagumo residual network initialized in the Turing pattern regime\. A single neuron’s activator potentials are plotted with respect to depth\. The network simulated was 30 layers deep and 64 neurons in width\.∼\\sim30 layers was the maximum forward integration depth, beyond which the forward integration diverges\. We note that instability of the forward dynamics of deep EBMs has been observed during training of very deep modelsInnocentiet al\.\([2026](https://arxiv.org/html/2605.21568#bib.bib36)\)\.To show that the Hamiltonian spatial integration of the FHN network from the initial conditions\(u0,p0\)\(u\_\{0\},p\_\{0\}\)and\(v0,q0\)\(v\_\{0\},q\_\{0\}\)indeed computes the same solution as the iterative time\-based simulation of the steady state, we first perform a time\-based simulation of the above equations to calculate the steady\-state solution\. We then take the first layer’s node values\(u0,v0\)\(u\_\{0\},v\_\{0\}\)along with its residual branch values\(p0,q0\)\(p\_\{0\},q\_\{0\}\), and use those as the input to the following forward recursion:

ui\+1\\displaystyle u^\{i\+1\}=ui\+pi\\displaystyle=u^\{i\}\+p^\{i\}vi\+1\\displaystyle v^\{i\+1\}=vi\+qi\\displaystyle=v^\{i\}\+q^\{i\}pi\+1\\displaystyle p^\{i\+1\}=Mi​Ni​pi\+Mi​fi\+1−Mi​Oi​ui\+1\\displaystyle=M^\{i\}N^\{i\}p^\{i\}\+M^\{i\}f^\{i\+1\}\-M^\{i\}O^\{i\}u^\{i\+1\}qi\+1\\displaystyle q^\{i\+1\}=ε​\(ui\+1−α​vi\+1−β\)\\displaystyle=\\varepsilon\(u^\{i\+1\}\-\\alpha v^\{i\+1\}\-\\beta\)
To maintain the correct boundary condition at the initial input \(where we have a ghost node\), we must modify the first step in the recursion by multiplying by a factor of1/21/2:

p1\\displaystyle p^\{1\}=12​\[M0​N0​p0\+M0​f1−M0​O0​u1\]\\displaystyle=\\frac\{1\}\{2\}\\left\[M^\{0\}N^\{0\}p^\{0\}\+M^\{0\}f^\{1\}\-M^\{0\}O^\{0\}u^\{1\}\\right\]q1\\displaystyle q^\{1\}=12​\[ε​\(u1−α​v1−β\)\]\\displaystyle=\\frac\{1\}\{2\}\\left\[\\varepsilon\(u^\{1\}\-\\alpha v^\{1\}\-\\beta\)\\right\]
The results of these simulations are shown in Fig\. 2\. Note that for networks with a depth of under approximately 30 layers, the spatial integration stays close to the time\-based steady\-state\. However, for deeper networks, the solution eventually diverges, because the forward spatial integration is along an unstable manifold for the Fitzhugh\-Nagumo model\. More advanced spatial integration schemes or modifying the network architecture, for example to include layer\-wise normalization, may alleviate this issue\.

We note that the above analysis is performed simply to show that the time\-based and layer\-wise simulations are equivalent\. More practical methods for performing inference directly in the Hamiltonian framework are available, e\.g\. inPourcel and Ernoult \([2025](https://arxiv.org/html/2605.21568#bib.bib28)\)\. A discussion of the Hamiltonian approach to inference in deep Energy\-Based Models will be discussed in a separate work\.

## VIConclusion

We have shown that skew\-gradient systems, such as the Fitzhugh\-Nagumo model, are compatible with the Equilibrium Propagation algorithm at their steady state solutions\. Additionally, using techniques from the study of skew\-gradient systems, we have also shown how, in deep Energy\-Based Models, one can derive a Hamiltonian recurrence relation which allows for exact layer\-wise computation of the steady\-state solution, provided an additional initial condition on the momentum is given\. Thus, the Energy\-Based and Hamiltonian formulations differ only via the definition of the boundary conditions for the system: the Energy\-Based formulation is a two\-point boundary value problem, and the Hamiltonian formulation is a purely initial\-value problem\. Having an explicit bridge between the Energy\-Based \(spatial Lagrangian\) and Hamiltonian formulations of a network opens the door to applying many useful methods of analysis which rely on a Hamiltonian form to deep Energy\-Based Models\. Most notably, this includes bifurcation analysis of the uniform network’s spatial orbits\. We believe this work will provide a useful new tool for training deep networks of activator\-inhibitor neurons, and give additional insight into the widespread occurrence of activator\-inhibitor structure in pattern formation and biological networks\.

## References

- L\. E\. Altman, M\. Stern, A\. J\. Liu, and D\. J\. Durian \(2024\)Experimental demonstration of coupled learning in elastic networks\.Physical Review Applied22\(2\),pp\. 024053\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p2.1)\.
- Y\. Bengio and A\. Fischer \(2015\)Early inference in energy\-based models approximates back\-propagation\.arXiv preprint arXiv:1510\.02777\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p1.1),[§I\-C](https://arxiv.org/html/2605.21568#S1.SS3.p5.1),[§IV\-B](https://arxiv.org/html/2605.21568#S4.SS2.p1.1)\.
- R\. Bogacz \(2017\)A tutorial on the free\-energy framework for modelling perception and learning\.Journal of mathematical psychology76,pp\. 198–211\.Cited by:[§I\-A](https://arxiv.org/html/2605.21568#S1.SS1.p3.1)\.
- D\. Cebrián\-Lacasa, P\. Parra\-Rivas, D\. Ruiz\-Reynés, and L\. Gelens \(2024\)Six decades of the fitzhugh–nagumo model: a guide through its spatio\-temporal dynamics and influence across disciplines\.Physics Reports1096,pp\. 1–39\.Cited by:[§I\-C](https://arxiv.org/html/2605.21568#S1.SS3.p2.1)\.
- Y\. Du, S\. Li, and I\. Mordatch \(2020\)Compositional visual generation with energy based models\.Advances in Neural Information Processing Systems33,pp\. 6637–6647\.Cited by:[§IV](https://arxiv.org/html/2605.21568#S4.p1.1)\.
- R\. FitzHugh \(1961\)Impulses and physiological states in theoretical models of nerve membrane\.Biophysical journal1\(6\),pp\. 445–466\.Cited by:[§I\-C](https://arxiv.org/html/2605.21568#S1.SS3.p2.1)\.
- K\. Friston \(2003\)Learning and inference in the brain\.Neural networks16\(9\),pp\. 1325–1352\.Cited by:[§I\-A](https://arxiv.org/html/2605.21568#S1.SS1.p3.1)\.
- S\. Grossberg \(1987\)Competitive learning: from interactive activation to adaptive resonance\.Cognitive science11\(1\),pp\. 23–63\.Cited by:[§I](https://arxiv.org/html/2605.21568#S1.p2.1)\.
- R\. Høier, T\. Kerjan, and B\. Scellier \(2026\)Training a convergent energy transformer with equilibrium propagation\.InNew Frontiers in Associative Memories\-Workshop at ICLR 2026,Cited by:[§I\-A](https://arxiv.org/html/2605.21568#S1.SS1.p2.1)\.
- T\. W\. Hughes, M\. Minkov, Y\. Shi, and S\. Fan \(2018\)Training of photonic neural networks through in situ backpropagation and gradient measurement\.Optica5\(7\),pp\. 864–871\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p2.1)\.
- F\. Innocenti, E\. M\. Achour, and C\. L\. Buckley \(2026\)μ\\mu\-PC: scaling predictive coding to 100\+ layer networks\.Advances in Neural Information Processing Systems38,pp\. 87414–87455\.Cited by:[Figure 2](https://arxiv.org/html/2605.21568#S5.F2),[Figure 2](https://arxiv.org/html/2605.21568#S5.F2.2.1)\.
- B\. Karasözen, T\. Küçükseyhan, and M\. Uzunca \(2017\)Structure preserving integration and model order reduction of skew\-gradient reaction–diffusion systems\.Annals of Operations Research258\(1\),pp\. 79–106\.Cited by:[§IV\-A](https://arxiv.org/html/2605.21568#S4.SS1.p1.4)\.
- J\. Kendall, R\. Pantone, K\. Manickavasagam, Y\. Bengio, and B\. Scellier \(2020\)Training end\-to\-end analog neural networks with equilibrium propagation\.arXiv preprint arXiv:2006\.01981\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p2.1)\.
- D\. Krotov and J\. J\. Hopfield \(2016\)Dense associative memory for pattern recognition\.Advances in neural information processing systems29\.Cited by:[§I\-A](https://arxiv.org/html/2605.21568#S1.SS1.p1.1)\.
- M\. Kuwamura \(2004\)On the turing patterns in one\-dimensional gradient/skew\-gradient dissipative systems\.SIAM Journal on Applied Mathematics65\(2\),pp\. 618–643\.Cited by:[§IV\-A](https://arxiv.org/html/2605.21568#S4.SS1.p1.4),[§IV\-A](https://arxiv.org/html/2605.21568#S4.SS1.p2.1),[§IV\-A](https://arxiv.org/html/2605.21568#S4.SS1.p6.1)\.
- A\. Laborieux, M\. Ernoult, B\. Scellier, Y\. Bengio, J\. Grollier, and D\. Querlioz \(2021\)Scaling equilibrium propagation to deep convnets by drastically reducing its gradient estimator bias\.Frontiers in neuroscience15,pp\. 633674\.Cited by:[§I\-A](https://arxiv.org/html/2605.21568#S1.SS1.p2.1)\.
- J\. Laydevant, D\. Marković, and J\. Grollier \(2024\)Training an ising machine with equilibrium propagation\.Nature Communications15\(1\),pp\. 3671\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p2.1)\.
- J\. H\. Lee, T\. Delbruck, and M\. Pfeiffer \(2016\)Training deep spiking neural networks using backpropagation\.Frontiers in neuroscience10,pp\. 508\.Cited by:[§I](https://arxiv.org/html/2605.21568#S1.p1.1)\.
- T\. P\. Lillicrap, A\. Santoro, L\. Marris, C\. J\. Akerman, and G\. Hinton \(2020\)Backpropagation and the brain\.Nature Reviews Neuroscience21\(6\),pp\. 335–346\.Cited by:[§I](https://arxiv.org/html/2605.21568#S1.p2.1),[§II](https://arxiv.org/html/2605.21568#S2.p28.3),[§II](https://arxiv.org/html/2605.21568#S2.p29.1)\.
- V\. Lopez\-Pastor and F\. Marquardt \(2023\)Self\-learning machines based on hamiltonian echo backpropagation\.Physical Review X13\(3\),pp\. 031020\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p3.1)\.
- E\. Martin, M\. Ernoult, J\. Laydevant, S\. Li, D\. Querlioz, T\. Petrisor, and J\. Grollier \(2021\)Eqspike: spike\-driven equilibrium propagation for neuromorphic implementations\.Iscience24\(3\)\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p4.1)\.
- B\. Millidge, A\. Seth, and C\. L\. Buckley \(2021\)Predictive coding: a theoretical and experimental review\.arXiv preprint arXiv:2107\.12979\.Cited by:[§I\-A](https://arxiv.org/html/2605.21568#S1.SS1.p3.1)\.
- B\. Millidge, Y\. Song, T\. Salvatori, T\. Lukasiewicz, and R\. Bogacz \(2022\)Backpropagation at the infinitesimal inference limit of energy\-based models: unifying predictive coding, equilibrium propagation, and contrastive hebbian learning\.arXiv preprint arXiv:2206\.02629\.Cited by:[§I\-A](https://arxiv.org/html/2605.21568#S1.SS1.p3.1)\.
- R\. Muolo, L\. Giambagli, H\. Nakao, D\. Fanelli, and T\. Carletti \(2024\)Turing patterns on discrete topologies: from networks to higher\-order structures\.Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences480\(2302\)\.Cited by:[§II](https://arxiv.org/html/2605.21568#S2.p7.1)\.
- J\. Nagumo, S\. Arimoto, and S\. Yoshizawa \(1962\)An active pulse transmission line simulating nerve axon\.Proceedings of the IRE50\(10\),pp\. 2061–2070\.Cited by:[§I\-C](https://arxiv.org/html/2605.21568#S1.SS3.p2.1)\.
- P\. Parra\-Rivas, F\. A\. Saadi, and L\. Gelens \(2025\)Spatial localization in the fitzhugh\-nagumo model\.arXiv preprint arXiv:2501\.10271\.Cited by:[§IV\-A](https://arxiv.org/html/2605.21568#S4.SS1.p1.4),[§IV\-A](https://arxiv.org/html/2605.21568#S4.SS1.p23.1),[§IV](https://arxiv.org/html/2605.21568#S4.p2.1)\.
- G\. Pourcel and M\. Ernoult \(2025\)Learning long range dependencies through time reversal symmetry breaking\.arXiv preprint arXiv:2506\.05259\.Cited by:[§V](https://arxiv.org/html/2605.21568#S5.p5.1)\.
- B\. Scellier and Y\. Bengio \(2017\)Equilibrium propagation: bridging the gap between energy\-based models and backpropagation\.Frontiers in computational neuroscience11,pp\. 24\.Cited by:[§I\-A](https://arxiv.org/html/2605.21568#S1.SS1.p1.1),[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p1.1),[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p2.1),[§I\-C](https://arxiv.org/html/2605.21568#S1.SS3.p5.1),[§III](https://arxiv.org/html/2605.21568#S3.p1.1),[§IV\-B](https://arxiv.org/html/2605.21568#S4.SS2.p1.1)\.
- B\. Scellier, M\. Ernoult, J\. Kendall, and S\. Kumar \(2023\)Energy\-based learning algorithms for analog computing: a comparative study\.Advances in neural information processing systems36,pp\. 52705–52731\.Cited by:[§I\-A](https://arxiv.org/html/2605.21568#S1.SS1.p2.1)\.
- B\. Scellier \(2024\)Quantum equilibrium propagation: gradient\-descent training of quantum systems\.arXiv preprint arXiv:2406\.00879\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p2.1)\.
- S\. Smale \(1971\)On the mathematical foundations of electrical circuit theory\.InThe Collected Papers of Stephen Smale: Volume 2,pp\. 951–968\.Cited by:[§IV\-A](https://arxiv.org/html/2605.21568#S4.SS1.p28.1)\.
- M\. Stern and A\. Murugan \(2023\)Learning without neurons in physical systems\.Annual Review of Condensed Matter Physics14\(1\),pp\. 417–441\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p2.1)\.
- C\. C\. Wanjura and F\. Marquardt \(2025\)Quantum equilibrium propagation for efficient training of quantum systems based on onsager reciprocity\.Nature Communications16\(1\),pp\. 6595\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p2.1)\.
- T\. C\. Wunderlich and C\. Pehle \(2021\)Event\-based backpropagation can compute exact gradients for spiking neural networks\.Scientific Reports11\(1\),pp\. 12829\.Cited by:[§I](https://arxiv.org/html/2605.21568#S1.p1.1)\.
- E\. Yanagida \(2002a\)Mini\-maximizers for reaction\-diffusion systems with skew\-gradient structure\.Journal of Differential Equations179\(1\),pp\. 311–335\.Cited by:[§I\-C](https://arxiv.org/html/2605.21568#S1.SS3.p4.1),[§II](https://arxiv.org/html/2605.21568#S2.p1.1),[§II](https://arxiv.org/html/2605.21568#S2.p19.6),[§IV\-A](https://arxiv.org/html/2605.21568#S4.SS1.p1.4)\.
- E\. Yanagida \(2002b\)Standing pulse solutions in reaction\-diffusion systems with skew\-gradient structure\.Journal of Dynamics and Differential Equations14\(1\),pp\. 189–205\.Cited by:[§I\-C](https://arxiv.org/html/2605.21568#S1.SS3.p4.1),[§II](https://arxiv.org/html/2605.21568#S2.p1.1),[§IV\-A](https://arxiv.org/html/2605.21568#S4.SS1.p1.4),[§IV\-A](https://arxiv.org/html/2605.21568#S4.SS1.p16.1)\.
- S\. Yi, J\. D\. Kendall, R\. S\. Williams, and S\. Kumar \(2023\)Activity\-difference training of deep neural networks using memristor crossbars\.Nature Electronics6\(1\),pp\. 45–51\.Cited by:[§I\-B](https://arxiv.org/html/2605.21568#S1.SS2.p4.1)\.