Towards training digitally-tied analog blocks via hybrid gradient computation
Timothy Nest Maxence Ernoult



Quantum Algorithm For Testing Convexity of Function
    Tzu-Chieh Wei
    Received 16 July 2024; accepted 04 September 2024
=====================================================




§ ABSTRACT 
⧫Equal contributionfootnote
Power efficiency is plateauing in the standard digital electronics realm such that novel hardware, models, and algorithms are needed to reduce the costs of AI training. The combination of energy-based analog circuits and the Equilibrium Propagation (EP) algorithm constitutes one compelling alternative compute paradigm for gradient-based optimization of neural nets. Existing analog hardware accelerators, however, typically incorporate digital circuitry to sustain auxiliary non-weight-stationary operations, mitigate analog device imperfections, and leverage existing digital accelerators.This heterogeneous hardware approach calls for a new theoretical model building block. In this work, we introduce Feedforward-tied Energy-based Models (ff-EBMs), a hybrid model comprising feedforward and energy-based blocks accounting for digital and analog circuits. We derive a novel algorithm to compute gradients end-to-end in ff-EBMs by backpropagating and “eq-propagating” through feedforward and energy-based parts respectively, enabling EP to be applied to much more flexible and realistic architectures. We experimentally demonstrate the effectiveness of the proposed approach on ff-EBMs where Deep Hopfield Networks (DHNs) are used as energy-based blocks. We first show that a standard DHN can be arbitrarily split into any uniform size while maintaining performance. We then train ff-EBMs on ImageNet32 where we establish new SOTA performance in the EP literature (46 top-1 %). Our approach offers a principled, scalable, and incremental roadmap to gradually integrate self-trainable analog computational primitives into existing digital accelerators.





§ INTRODUCTION

Gradient-based optimization, the cornerstone and most energy greedy component of deep learning, fundamentally relies upon three factors: i) highly parallel digital hardware such as GPUs, ii) feedforward models and iii) backprop (BP). With skyrocketing demands of AI compute, cutting the energy consumption of AI systems, learning has become a economical, societal and environmental stake <cit.> and calls for the exploration of novel compute paradigms <cit.>.


One promising path towards this goal is analog in-memory computing  <cit.>: when mapping weights onto a crossbar of resistive devices, Kirchoff current and voltage laws inherently achieve matrix-vector multiplications in constant time complexity <cit.>.
Stacking multiple such crossbars, an entire neural network can be mapped onto a physical system. An important formalism for such systems is that of energy-based (EB) analog circuits <cit.> which are “self-learning” systems that compute loss gradients through two relaxations to equilibrium (i.e. two “forward passes”), a procedure falling under the umbrella of energy-based learning (EBL) algorithms <cit.>. 

One of these learning algorithms, Equilibrium Propagation (EP) <cit.>, particularly stands out with strong theoretical guarantees, relative scalability in the realm of backprop alternatives <cit.> and experimental demonstrations on small analog systems which are 10,000× more energy-efficient and substantially faster than their GPU-based counterpart <cit.>. This suggests an alternative triad as a new compute paradigm for gradient-based optimization: i) analog hardware, ii) EBMs, iii) EP. 
[22]l0.4
    

    
    < g r a p h i c s >

    

    Illustrating BP-EP backward gradient chaining through feedforward (red) and energy-based (yellow) blocks, accounting for digital and analog circuits respectively.

In this paper, we propose a theoretical framework to extend end-to-end gradient computation to a realistic setting where the system at use may or may not be fully analog. Such a setting is plausible in the near term, due to by two major limitations. First, analog circuits exhibit many non-ideal physical behaviors which affect both the inference pathway <cit.> and parameter optimization <cit.> ,  in-turn compromising performance. Second, owing to the latency and energy-consumption of resistive devices' write operations, such analog circuits should be fully weight stationary – weights must be written before the inference procedure begins – which excludes many operations used conventionally in machine learning such as activation functions, normalization, attention <cit.>. Therefore, analog systems are likely to be used in combination with auxiliary digital circuitry, resulting in hybrid mixed precision systems <cit.>. While the design of purely inferential engines made of analog and digital parts is nearing commercial maturity <cit.>, in-situ learning of such systems has barely been explored. One final challenge lies in proving EBL algorithms can scale in a manner comparable to backprop, given the requirement of simulating EB systems on GPUs. Given the necessity of convergence, this amounts in practice in performing lengthy root finding algorithms to simulate physical equilibrium, limiting proof-of-concepts thereof to relatively shallow models <cit.>.

Our work contends that the best of both worlds can be achieved with the following triad: i) hybrid digital and analog hardware, ii) feedforward and EB models, iii) BP and EP. Namely, by modelling digital and analog parts as feedforward and EB modules respectively, the core contribution of our paper is to show how backprop and EP error signals can be chained end-to-end through feedforward and EB blocks respectively in a principled fashion. Rather than opposing digital and analog, or backprop and “alternative” learning algorithms as often done in the literature, we propose a novel hardware-aware building block which can, in principle, leverage advances from both digital and analog hardware in the near-term. More specifically:


    
  * We propose Feedforward-tied Energy-based Models (ff-EBMs, Section <ref>) as high-level models of mixed precision systems whose inference pathway read as the composition of feedforward and EB modules (Eq. (<ref>), Alg. <ref>). 
    
  * We show that gradients in ff-EBMs can be computed in an end-to-end fashion (Section <ref>), backpropagating through feedforward blocks and “eq-propagating” through EB blocks (Theorem <ref>, Alg. <ref>) and that this procedure is rooted in a deeply-nested optimization problem (Section <ref>). 
    
    
  * Finally, we experimentally demonstrate the effectiveness of our algorithm on ff-EBMs where EBM blocks are Deep Hopfield Networks (DHNs) (Section <ref>). We show that i) gradient estimates computed by our algorithm (Alg. <ref>) near perfectly match gradients computed by end-to-end automatic differentiation (Section <ref>), ii) a standard DHN model can be arbitrarily split into a ff-DHN with the equivalent layers and architectural  layers while maintaining performance and remaining on par with automatic differentiation (Section <ref>), iii)  the proposed approach yields 46 % top-1 (70% top-5) validation accuracy on ImageNet32 when training a ff-EBM of 16 layers, thereby significantly beating EP current performance state-of-the-art by a large margin without using holomorphic transformations inside EBM blocks <cit.> 




§ BACKGROUND




  
Notations. Denoting A: ℝ^n →ℝ^m a differentiable mapping, we denote its total derivative with respect to s_j as d_s_j A(s) := d A(s)/ ds_j ∈ℝ^m, its partial derivative with respect to s_j as ∂_jA(s) := ∂ A(s)/ ∂ s_j ∈ℝ^m. When A takes scalar values (m=1), its gradient with respect to s_j is denoted as ∇_j A(s) := ∂_j A(s)^⊤.



 §.§ Energy-based models (EBMs)

For a given static input and set of weights, Energy-based models (EBMs) implicitly yield a  prediction through the minimization of an energy function – as such they are a particular kind of implicit model. Namely, an EBM is defined by a (scalar) energy function E: s, θ, x → E(s, θ, x) ∈ℝ where x, s, and θ respectively denote a static input, hidden and output neurons and model parameters, and each such tuple defines a configuration with an associated scalar energy value. Amongst all configurations for a given input x and some model parameters θ, the model prediction s_⋆ is implicitly given as an equilibrium state which minimizes the energy function:


    s_⋆ := min_s E(s, θ, x).



 §.§ Standard bilevel optimization

Assuming that ∇_s^2 E(x, s_⋆, θ) is invertible, note that the equilibrium state s_⋆ implicitly depends on x and θ by virtue of the implicit function theorem <cit.>. Therefore our goal when training an EBM, for instance in a supervised setting, is to adjust the model parameters θ such that s_⋆(x, θ) minimizes some cost function ℓ: s, y →ℓ(s, y) ∈ℝ where y is some ground-truth label associated to x. More formally, this learning objective can be stated with the following bilevel optimization problem <cit.>:

    min_θ𝒞(x, θ, y) := ℓ(s_⋆, y)   s.t.    s_⋆ = min_s E(s, θ, x) .

Solving Eq. (<ref>) in practice amounts to computing the gradient of its outer objective 𝒞(x, θ) with respect to θ (d_θ𝒞(x, θ)) and then perform gradient descent over θ. 



 §.§ Equilibrium Propagation (EP)

An algorithm used to train an EBM model in the sense of Eq. (<ref>) may be called an EBL algorithm <cit.>. Equilibrium Propagation (EP) <cit.> is an EBL algorithm which computes an estimate of d_θ𝒞(x, θ) with at least two phases. During the first phase, the model is allowed to evolve freely to s_⋆ = min_s E(s, θ, x). Then, the model is slightly nudged towards decreasing values of cost ℓ and settles to a second equilibrium state s_β. This amounts to augment the energy function E by an additional term βℓ(s, y) where β∈ℝ^⋆ is called the nudging factor. Then the weights are updated to increase the energy of s_⋆ and decrease that of s_β, thereby “contrasting” these two states. More formally,
<cit.> prescribe in the seminal EP paper: 

    s_β := min_s [ E(s, θ, x) + βℓ(s, y)],   Δθ ^ EP := α/β(∇_2 E(s_⋆, θ, x) - ∇_2 E(s_β, θ, x) ),


where α denotes some learning rate. EP comes in different flavours depending on the sign of β inside Eq. (<ref>) or on whether two nudged states of opposite nudging strengths (±β) are contrasted, a variant called Centered EP (C-EP) which was shown to work best in practice <cit.> and reads as:


    Δθ ^ C-EP := α/2β(∇_2 E(s_-β, θ, x) - ∇_2 E(s_β, θ, x) ),




§ TYING ENERGY-BASED MODELS WITH FEEDFORWARD BLOCKS

This section mirrors the background section by introducing a new model, the naturally associated optimization problem and a new learning algorithm. We first introduce Feedforward-tied EBMs (ff-EBMs, section <ref>) which read as composition of feedforward and EB transformations (Alg. <ref>). We then show how optimizing ff-EBMs amounts to solving a multi-level optimization problem (Section <ref>) and propose a BP-EP gradient chaining algorithm as a solution (Section <ref>, Theorem <ref>, Alg. <ref>). We highlight as an edge case that ff-EBMs reduce to standard feedforward nets (Lemma <ref>) and the proposed BP-EP gradient chaining algorithm to standard BP (Corollary <ref>) when each EB block comprises a single hidden layer.




 §.§ Feedforward-tied Energy-based Models (ff-EBMs)





  
Inference procedure. We define Feedforward-tied Energy-based Models (ff-EBMs) as compositions of feedforward and EB transformations. Namely, an data sample x is fed into the first feedforward transformation F^1 parametrized by some weights ω^1, which yields an output x^1_⋆. Then, x^1_⋆ is fed as a static input into the first EB block E^1 with parameters θ^1, which relaxes to an equilibrium state s^1_⋆. s^1_⋆ is in turn fed into the next feedforward transformation F^1 with weights ω^1 and the above procedure repeats until reaching the output layer ô. More formally, denoting F^k and E^k the k^ th feedforward and EB blocks parametrized by the weights ω^k and θ^k respectively, the inference pathway of a ff-EBM reads as:


    {[                                                                        s^0 := x; x^k_⋆ := F^k(s^k-1_⋆, ω^k),    s^k_⋆ := smin E^k(s, θ^k, x^k_⋆)   ∀ k = 1 ⋯ N-1;                                                        ô_⋆ := F^N(s^N-1_⋆, ω^N) ].


ff-EBM inference procedure is depicted more compactly inside Fig. <ref> (left) and Alg. <ref>. 





[c]0.46

    < g r a p h i c s >

figureDepiction of the forward (left) and backward (right) pathways through a ff-EBM, with yellow and pink blocks denoting EB and feedforward transformations.


[c]0.5





  
Form of the energy functions. We specify further the form of the energy of the k^ th EB block of a ff-EBM as defined per Eq. (<ref>). The associated energy function E^k takes some static input x^k from the output of the preceding feedforward transformation, has hidden neurons s^k and is parametrized by weights θ^k and more precisely defined as:


    E^k(s^k, θ^k, x^k) := G^k(s^k) - s^k^⊤· x^k + U^k(s^k, θ^k)


Eq. (<ref>) reveals three different contributions to the energy. The first term determines the non-linearity applied inside the EB block <cit.>: for a given invertible and continuous activation function σ, G is defined such that ∇ G = σ^-1 (see Appendix <ref>).

The second term inside Eq. (<ref>) accounts for a purely feedforward contribution from the previous feedforward block F^k. Finally, the third term accounts for internal interactions within the layers of the EB block. 



  
Recovering a feedforward net. When taking the gradient of E^k as defined in Eq. (<ref>) with respect to s^k and zeroing it out, it can be seen that s^k_⋆ is implicitly defined as:


    s^k_⋆ := σ(x^k - ∇_1 U^k(s^k_⋆, θ^k) )


An interesting edge case highlighted by Eq. (<ref>) is when U^k = 0 for all k's, i.e. when there are no intra-block layer interactions, or equivalently when the EB block comprises a single layer only. In this case, s_⋆^k is simply a feedforward mapping x^k through σ and in turn the ff-EBM is simply a standard feedforward architecture (see Lemma <ref> inside Appendix <ref>).



 §.§ Multi-level optimization of ff-EBMs



In the same way as learning EBMs can naturally be cast into a bilevel optimization problem, learning ff-EBMs can be inherently be mapped into a multi-level optimization problem where the variables being optimized over in the inner subproblems are the EB block variables s^1, ⋯, s^N-1. To make this clearer, we re-write the energy function of the k^ th block E^k from Eq. (<ref>) to highlight the dependence between two consecutive EB block states:


    E^k(s^k, θ^k, s^k-1_⋆, ω^k) := E^k(s^k, θ^k, F^k(s^k-1_⋆, ω^k-1))


It can be seen from Eq. (<ref>) that the equilibrium state s^k_⋆ obtained by minimizing E^k will be dependent upon the equilibrium state s^k-1_⋆ of the previous EB block, which propagates back through prior EB blocks. Denoting W := {θ^1, ⋯, θ^N-1, ω^1, ⋯, ω^N}, the learning problem for a ff-EBM can therefore be written as:










    min_W    𝒞(x, W, y) :=ℓ(ô_⋆ = F^N(s^N-1_⋆, ω^N), y) 
    s.t.      s^N-1_⋆ = smin E^N-1(s, θ^N-1, s^N-2_⋆, ω^N-1)   ⋯  s.t.   s^1_⋆ = smin E^1(s, θ^1, x, ω^1)


Here again and similarly to bilevel optimization, solving Eq. (<ref>) in practice amounts to computing g_θ^k := d_θ^k𝒞 and  g_ω^k := d_ω^k𝒞 and perform gradient descent on θ^k and ω^k.



 §.§ A BP–EP gradient chaining algorithm





  
Main result: explicit BP-EP chaining. Based on the multilevel optimization formulation of ff-EBMs learning in Eq. (<ref>), we state the main theoretical result of this paper in Theorem <ref> (see proof in Appendix <ref>). 



Assuming a model of the form Eq. (<ref>), we denote s^1_⋆, x^1_⋆, ⋯, s^N-1_⋆, ô_⋆ the states computed during the forward pass as depicted in Alg. <ref>. 
We define the nudged state of block k, denoted as s^k_β,  implicitly through ∇_1 ℱ^k(s^k_β, θ^k, x^k_⋆, δ s^k, β) = 0 with:

    ℱ^k(s^k, θ^k, x^k_⋆, δ s^k, β) := E^k(s^k, θ^k, x^k_⋆) + β s^k^⊤·δ s^k


Denoting δ s^k and Δ x^k the error signals computed at the input of the feedforward block F^k and of the EB block E^k respectively, then the following chain rule applies:


    δ s^N-1 := ∇_s^N-1ℓ(ô_⋆, y),    g_ω^N = ∇_ω^Nℓ(ô_⋆, y) 
       ∀ k=2 ⋯ N-1: 
       {[ Δ x^k = .d_β(∇_3E^k (s^k_β, θ^k, x^k_⋆))|_β=0,    g_θ^k = .d_β(∇_2E^k (s^k_β, θ^k, x^k_⋆))|_β=0;              δ s^k-1 = ∂_1 F^k(s^k-1_⋆, ω^k )^⊤·Δ x^k,    g_ω^k = ∂_2 F^k(s^k-1_⋆, ω^k)^⊤·Δ x^k ].





  
Proposed algorithm: implicit BP-EP chaining. Theorem <ref> reads intuitively: it prescribes an explicit chaining of EP error signals passing backward through E^k (δ s^k →Δ x^k) and BP error signals passing backward through ∂ F^k^⊤ (Δ x^k →δ s^k-1), which directly mirrors the ff-EBM inference pathway as depicted in Fig. <ref>. Yet noticing that:




    {[ δ s^k-1 =  ∂_1 F^k(s^k-1_⋆, ω^k )^⊤·Δ x^k = d_β.(∇_3 E^k(s^k_β, θ^k, s^k-1_⋆, ω^k))|_β=0,;   g_ω^k = ∂_2 F^k(s^k-1_⋆, ω^k )^⊤·Δ x^k =  d_β.(∇_4 E^k(s^k_β, θ^k, s^k-1_⋆, ω^k))|_β=0, ].


the same error signal can by passed through E^k (δ s^k →δ s^k-1) where BP and EP are implicitly chained inside E^k (see Appendix <ref>). This insight, along with a centered scheme to estimate derivatives with respect to β around 0 as done for the C-EP algorithm (Eq. (<ref>)), motivates the implicit BP-EP gradient chaining algorithm in Alg. <ref> we used for our experiments (see Alg. <ref> inside Appendix <ref> for its explicit counterpart). For simplicity for here onwards and as the proposed algorithm appears to be a generalization of EP, we may refer to Alg. <ref> as “EP” in the experimental section.





  
Recovering backprop. When the ff-EBM under consideration is purely feedforward (U^k=0), we show that Eqs. (<ref>)–(<ref>) reduce to standard BP through a feedforward net (Corollary <ref>, Alg. <ref> and Alg. <ref> in Appendix <ref>). Therefore, since this case is extremely close to standard BP through feedforward nets, we do not consider this setting in our experiments.



§ EXPERIMENTS


In this section, we first present the ff-EBMs at use in our experiments (Section <ref>) and carry out static gradient analysis – computing and analyzing ff-EBM parameter gradients for some x and y (Section <ref>). We extend the observation made by <cit.> to ff-EBMs that transient EP parameter gradients throughout the second phase match those computed by automatic differentiation through equilibrium and across blocks (Fig. (<ref>)), with resulting final gradient estimates near perfectly aligned (Fig. <ref>). We then show on the CIFAR-10 task that performance of ff-EBMs can be maintained across various block splits and on par with automatic differentiation while keeping the same number of layers (Section <ref>). Finally, we perform further ff-EBM training experiments on CIFAR-100 and ImageNet32 where we establish a new performance state-of-the-art in the EP literature (Section <ref>). 


 §.§ Setup




  
Model. Using the same notations as in Eq. (<ref>), the ff-EBMs at use in this section are defined through: 

    {[    U^k_ FC(s^k, θ^k) := -1/2s^k⊤·θ^k · s^k,; U^k_ CONV(s^k, θ^k) := -1/2s^k ∙(θ^k ⋆ s^k) ]. ,    F^k(s^k-1, ω^k) :=  BN(𝒫(ω^k_ CONV⋆ s^k-1_L); ω^k_α, ω^k_β)

with BN(·; ω_α^k  ω_β^k), 𝒫 and ⋆ the batchnorm, pooling and convolution operations, ∙ the generalized dot product for tensors and s^k := (s^k^⊤_1, ⋯ s^k^⊤_L)^⊤ the state of block k comprising L layers. The EBM blocks are usually called Deep Hopfield Networks (DHNs) and the weight matrix θ^k is symmetric and has a sparse, block-wise structure such that each layer s^k_ℓ is bidirectionally connected to its neighboring layers s^k_ℓ - 1 and s^k_ℓ + 1 through connections θ^k_ℓ -1 and θ^k^⊤_ℓ respectively (see Appendix <ref>), either with fully connected (U^k_ FC) and convolutional operations (U^k_ CONV). Finally, the non-linearity σ applied within EB blocks is σ(x) := min(max(x/2, 0), 1). 



  
Equilibrium computation. As depicted in Alg. <ref>, the steady states s_±β may be computed with any loss minimization algorithm. Here and as done in most past works on EP <cit.>, we employ a fixed-point iteration scheme to compute the EB blocks steady states. Namely, we iterate Eq. (<ref>) until reaching equilibrium (the same scheme is used for ff-EBM inference, Alg. <ref>, with β=0.):

    s^k_±β, t + 1σ(x^k - ∇_1 U^k(s^k_±β, t, θ^k) ∓βδ s^k )

We employ a scheme to asynchronously update even (s^k_2ℓ') and odd (s^k_2ℓ' + 1) layers <cit.> – see Appendix <ref>. 



  
Algorithm baseline. As an algorithmic baseline, we simply use automatic differentiation (AD) backward through the fixed-point iteration scheme Eq. (<ref>) with β=0 and directly initializing s^k_t=0 = s_⋆. This version of AD, where we backpropagate through equilibrium, is known as “Recurrent Backpropagation” <cit.> or Implicit Differentiation (ID). 



 §.§ Static comparison of EP and ID on ff-EBMs

In order to study the transient dynamics of ID and EP, we define, with W^k := {θ^k, ω^k}:

    {[                                             g^ ID_W^k(t) := ∑_k=0^T d_W^k(T - k)𝒞(x, W, y),; g^ EP_W^k(t) := 1/2β(∇_W^kE^k(s^k_β(t), W^k, s^k-1_⋆) - ∇_W^kE^k(s^k_-β(t), W^k, s^k-1_⋆)), ].

where s^±β(t) is computed from Eq. (<ref>) with the nudging error current δ s^k computed with Alg. <ref>, and T is the total number of iterations used for both ID and EP in the gradient computation phase.



 For a given block k, d_W^k(T - k)𝒞(x, W, y) is the “sensitivity” of the loss 𝒞 to parameter W^k at timestep T - k so that g^ AD_W^k(t) is a ID gradient truncated at T-t. Similarly, g^ EP_W^k(t) is an EP gradient truncated at t steps forward through the nudged phase. When T is sufficiently large, g^ ID_W^k(T) and g^ EP_W^k(T) converge to d_W^k𝒞(x, W, y). Fig. <ref> displays (g^ ID_W^k(t))_t ≥ 0 and (g^ EP_W^k(t))_t ≥ 0 on an heterogeneous ff-EBM of 6 blocks and 15 layers with blocks comprising 2 or 3 layers for a randomly selected sample x and its associated label y – see caption for a detailed description. It can be seen EP and ID error weight gradients qualitatively match very well throughout time, across layers and blocks. More quantitatively, we display the cosine similarity between the final EP and ID weight gradient estimate g^ ID_W^k(T) and g^ EP_W^k(T) for each layer and observe that EP and ID weight gradients are near perfectly aligned.



 §.§ Splitting experiment



For a given EBM (standard, single block) and a fixed number of layers, we ask whether block splitting of this EBM into a ff-EBM with multiple EB blocks affects training performance. We address this question with two different depths (L=6 and L=12 layers in total) and various block splits maintaining a total number of layers (e.g. for L=6,  1 block of 6 layers, 2 blocks of 3 layers, etc.) and display the results obtained on the CIFAR-10 task inside Table <ref>. We observe that the performance achieved by EP on the 6-layers deep EBM is maintained across 4 different block splits between 89% and 90% and is consistently on par with the ID baseline for each ff-EBM and with the literature on EBMs of same depth <cit.>. Similarly, we observe that the performance achieved by EP on ff-EBMs with a total of 12 layers is maintained around 92% with three different block sizes, matches ID performance on each of these and surpasses EP state-of-the art on CIFAR-10 <cit.>. Overall these results suggest the agnosticity of ff-EBMs to EB block sizes and therefore appear to be flexible in design. 



























 §.§ Scaling experiment


Finally, unlike Section <ref>, we now consider ff-EBMs of fixed block size 2 and train them with two different depths (L=12 and L=15) on CIFAR-100 and ImageNet32 by EP and ID and show the results obtained in Table <ref>. Here again we observe that EP matches ID performance on all models and tasks, ff-EBMs benefit from depth, and the performance obtained by training the 16-layers deep ff-EBM by EP exceeds state-of-the-art performance on ImageNet32 by around 10% top-1 validation accuracy <cit.> and by around 5% the best performance reported on this benchmark among backprop alternatives <cit.>. 





§ DISCUSSION




  
Related work.  Since fixed-point iteration schemes were proposed to facilitate EP experiments <cit.>, there is a growing body of work revolving around algorithmic extensions of EP and assessments of its scalability on vision tasks. Most notably, <cit.> introduced a holomorphic version of EP where loss gradients are computed with adiabatic oscillations of the model through nudging in the complex plane, and
was very recently extended to more general implicit models <cit.>. Moving further towards physical implementations of EP, <cit.> proposed a fully black-box version of EP where details about the system may not be known. All these advances could be readily applied inside our EP-BP chaining algorithm to EB blocks. The work closest to ours, albeit with a purely theoretical motivation and without clear algorithmic prescriptions, is that of <cit.> where feedforward model learning is cast into a deeply nested optimization where consecutive layers are tied by elemental pair-wise energy functions, which more recently inspired the Dual Propagation algorithm <cit.>. This setting can be construed as a particular case of ff-EBM learning by EP where each EB block comprises a single layer (U^k = 0 inside Eq. (<ref>) which, however, remains extremely similar to BP (see last paragraph of Section <ref>).



  
Limitations and future work. Since our recipe advocates EP–BP chaining by construction, it is fair to say that ff-EBM learning partially inherits the pitfalls of BP.

Fortunately, nothing prevents feedforward modules inside ff-EBMs to be trained by any BP alternative to mitigate specific issues. For instance: BP can be parameterized by feedback weights to obviate weight transport from the inference circuit to the gradient computation circuit <cit.>; BP gradients can be approximated as finite differences of feedback operators <cit.>; or computed via implicit forward-mode differentiation by applying random weight perturbations in the inference circuit <cit.>; local layer-wise self-supervised or supervised loss functions can be used to prevent “backward locking” <cit.>. This insight may help exploring many variants of ff-EBM training.





Pursuing the core motivation of this work, one natural extension of this study is to incorporate more hardware realism into ff-EBMs. Beyond Deep Hopfield networks, Deep Resistive Nets (DRNs) – concurrently developed by <cit.> and strongly inspired by <cit.> – are exact models of idealized analog circuits, are fast to simulate and were shown to be trainable by EP. As such, using DRNs as EB blocks inside ff-EBMs is an exciting research direction. Yet, going further into analog hardware modeling for ff-EBMs comes with new challenges when taking into account device non-idealities which may affect the inference pathway, such as analog-to-digital and digital-to-analog noise <cit.>.


Finally, considerable work is needed to prove ff-EBM further at scale on more difficult tasks (e.g. standard ImageNet), considerably deeper architectures and beyond vision tasks. One other exciting research direction would be the design of ff-EBM based transformers, with attention layers being chained with energy-based fully connected layers inside attention blocks.



  
Concluding remarks and broader impact. We show that ff-EBMs constitute a novel framework for deep-learning in heterogeneous hardware settings. We hope that the  algorithm proposed can help to move beyond the typical division between digital versus analog or BP versus BP-free algorithms and that the greater energy-efficiency afforded by this framework provides a more pragmatic, near-term blueprint to mitigate the dramatic carbon footprint of AI training <cit.>. Being still a long way from fully analog training accelerators at commercial maturity, we believe this work offers an incremental and sustainable roadmap to gradually integrate analog, energy-based computational primitives as they are developed into existing digital accelerators. 





§ ACKNOWLEDGEMENTS AND DISCLOSURE OF FUNDING

The authors warmly thank Irina Rish, Jack Kendall and Suhas Kumar for their support of the project idea from the very start, as well as Gregory Kollmer and Mohammed Fouda for useful feedback on the manuscript. TN acknowledges the support from the Canada Excellence Research Chairs Program, as well as CIFAR and Union Neurosciences et Intelligence Artificielle Quebec (UNIQUE). This research was enabled by the computational resources provided by the Summit supercomputer, awarded through the Frontier DD allocation and INCITE 2023 program for the project "Scalable Foundation Models for Transferable Generalist AI" and SummitPlus allocation in 2024. These resources were supplied by the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, with support from the Office of Science of the U.S. Department of Energy. ME acknowledges funding from Rain AI which commercializes technologies based on brain-inspired learning algorithms, as well as Constance Castres Saint-Martin for her unwavering support at the maternity hospital where most of this manuscript was written. 

abbrvnat








§ APPENDIX






 §.§ Further insights on ff-EBMs



In this section:


    
  * We formally define Feedforward-tied Energy-based Models (ff-EBMs) with precise assumptions on the energy-based and feedforward blocks (Def. <ref>).
    
  * We show that when energy-based blocks comprise a single layer only, the ff-EBM becomes purely feedforward (Lemma <ref>).




A Feedforward-tied Energy-based Model (ff-EBM) of size N comprises N twice differentiable feedforward mapping F^1, ⋯, F^N and N-1 twice differentiable energy functions E^1, ⋯, E^N-1. For a given x, the inference procedure reads as:
    
    {[                                                                        s^0 := x; x^k_⋆ := F^k(s^k-1_⋆, ω^k),    s^k_⋆ := smin E^k(s, θ^k, x^k_⋆)   ∀ k = 1 ⋯ N-1;                                                        ô_⋆ := F^N(s^N-1_⋆, ω^N) ].


Finally, we assume that ∀ k = 1 ⋯ N-1, ∇_1^2 E^k(s^k_⋆, θ^k, ω^k) is invertible.





We consider ff-EBM per Def. (<ref>) where the energy functions E^k have the form:


    E^k(s^k, θ^k, x^k) := G^k(s^k) - s^k^⊤· x^k + U^k(s^k, θ^k).


We assume that U^k = 0 for k = 1 ⋯ N-1, s →∇ G(s) is invertible and we denote σ := ∇ G^-1. Then, the resulting model is a feedforward model described by the following recursive equations:
    
    {[                                                      s^0_⋆ = x; x^k_⋆ = F^k(s^k-1_⋆, ω^k),    s^k_⋆ = σ(x^k_⋆)   ∀ k = 1 ⋯ N-1;                                       ô_⋆ := F^N(s^N-1_⋆, ω^N) ].




    Let k ∈ [1, N-1]. By definition of s^k_⋆ and x^k_⋆:

    
    ∇_1 E^k(s^k_⋆, θ^k, x^k_⋆) = 0 
    ⇔     ∇ G^k(s^k_⋆) - x^k_⋆ + ∇_1 U^k(s^k_⋆, θ^k) = 0 
    ⇔      s^k_⋆ = σ(x^k_⋆ - ∇_1 U^k(s^k_⋆, θ^k) )


    Therefore Eq. (<ref>) is immediately obtained from Eq. (<ref>) with U^k=0.
    



 §.§ Proof of Theorem <ref>



The proof of Theorem <ref> is structured as follows:


    
    
  * We directly solve the multilevel problem optimization defined inside Eq. (<ref>) using a Lagrangian-based approach (Lemma <ref>), yielding optimal Lagrangian multipliers, block states and loss gradients.
    
  * We show that by properly nudging the blocks, EP implicitly estimates the previously derived Lagrangian multipliers (Lemma <ref>).
    
  * We demonstrate Theorem <ref> by combining Lemma <ref> and Lemma <ref>.
    
  * Finally, we highlight that when a ff-EBM is a feedforward net (Lemma <ref>), then the proposed algorithm reduces to BP (Corollary <ref>).






Assuming a ff-EBM (Def. <ref>), we denote s^1_⋆, x^1_⋆, ⋯, s^N-1_⋆, ô_⋆ the states computed during the forward pass as prescribed by Eq. (<ref>). Then, the gradients of the objective function 𝒞:=ℓ(ô(s^N-1_⋆), y) as defined in the multilevel optimization problem (Eq. (<ref>)), where it is assumed that ℓ is differentiable, read:


    {[                         d_ω^N𝒞 = ∂_2F^N(s_⋆^N-1, ω^N)^⊤·∂_1 ℓ(ô_⋆, y),; d_θ^k𝒞 = ∇^2_1, 2E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k) ·λ_⋆^k   ∀ k = 1 ⋯ N-1,; d_ω^k𝒞 = ∇^2_1, 4E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k) ·λ_⋆^k   ∀ k = 1 ⋯ N-1, ].

where λ^1_⋆, ⋯, λ^N-1_⋆ satisfy the following conditions:

    {[                   ∇_s^N-1ℓ(ô(s^N-1_⋆), y) + ∇^2_1 E^N-1(s^N-1_⋆, θ^N-1, s^N-2_⋆, ω^N-1) ·λ_⋆^N-1 = 0;                                                                                      ∀ k =N-2, ⋯, 1:; ∇^2_1, 3E^k+1(s^k+1_⋆, θ^k+1, s^k_⋆, ω^k+1) ·λ_⋆^k+1 + ∇_1^2 E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k)·λ_⋆^k = 0 ].




Denoting s := (s^1, ⋯, s^N-1)^⊤ the state variables of the energy-based blocks, λ := (λ^1, ⋯, λ^N-1)^⊤ the Lagrangian multipliers associated with each of these variables, W := {θ_1, ω_1, ⋯, θ_N-1, ω_N-1} the energy-based and feedforward parameters and ô(s^N-1) := F^N(s^N-1, ω^N-1) the logits, the Lagrangian of the multilevel optimization problem as defined in Eq. (<ref>) reads:


    ℒ(s, λ, W) := ℓ(ô(s^N-1), y) + ∑_k=1^N-1λ^k^⊤·∇_1 E^k(s^k, θ^k, s^k-1, ω^k),    s^0 := x


Writing the associated Karush-Kuhn-Tucker (KKT) conditions ∂_1, 2ℒ(s_⋆, λ_⋆, W) := 0 satisfied by optimal states and Lagrangian multipliers s_⋆ and λ_⋆, we get :


    ∇_1 E^k (s^k_⋆, θ^k, s^k-1_⋆, ω^k) = 0   ∀ k = 1, ⋯, N- 1
       ∇_s^N-1ℓ(ô(s^N-1_⋆), y) + ∇^2_1 E^N-1(s^N-1_⋆, θ^N-1, s^N-2_⋆, ω^N-1) ·λ_⋆^N-1 = 0  
       ∇^2_1, 3E^k+1(s^k+1_⋆, θ^k+1, s^k_⋆, ω^k+1) ·λ_⋆^k+1 + ∇_1^2 E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k)·λ_⋆^k = 0   ∀ k =N-2, ⋯, 1


Eq. (<ref>) governs the bottom-up block-wise relaxation procedure (as depicted in Alg. <ref>), while Eq. (<ref>) and  Eq. (<ref>) governs error propagation in the last block and previous blocks respectively. Given s_⋆ and λ_⋆ by Eq. (<ref>) – Eq. (<ref>), the total derivative of the loss function with respect to the model parameters read:


    d_W ℓ(ô_⋆, y)    = d_W (ℓ(ô_⋆, y) + ∑_k=1^N-1λ^k^⊤_⋆·∇_1 E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k)_=0   ) 
       = d_W ℒ(s_⋆, λ_⋆, W) 
       = d_W s_⋆^⊤·∂_1 ℒ(s_⋆, λ_⋆, W)_=0   )) + d_W λ_⋆^⊤·∂_2 ℒ(s_⋆, λ_⋆, W)_=0    + ∂_3  ℒ(s_⋆, λ_⋆, W) 
       = ∂_3  ℒ(s_⋆, λ_⋆, W)


More precisely, applying Eq. (<ref>) to the feedforward and energy-based block parameters yields:


    d_ω^Nℓ(ô_⋆, y)    = ∂_2F^N(s_⋆^N-1, ω^N)^⊤·∇_1 ℓ(ô_⋆, y), 
    
        d_θ^kℓ(ô_⋆, y)    = ∇^2_1, 2E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k) ·λ_⋆^k   ∀ k = 1 ⋯ N-1 
    
        d_ω^kℓ(ô_⋆, y)    = ∇^2_1, 4E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k) ·λ_⋆^k   ∀ k = 1 ⋯ N-1

    





Under the same hypothesis as Lemma <ref>, we define the nudged state of block k, denoted as s^k_β,  implicitly through ∇_1 ℱ^k(s^k_β, θ^k, x^k_⋆, δ s^k, β) = 0 with:

    ℱ^k(s^k, θ^k, x^k_⋆, δ s^k, β) := E^k(s^k, θ^k, x^k_⋆) + β s^k^⊤·δ s^k.


Defining (δ s^k)_k=1 ⋯ N-1 recursively as:


    δ s^N-1 := ∇_s^N-1ℓ(ô_⋆, y),   δ s^k := d_β.(∇_3 E^k+1(s^k+1_β, θ^k+1, s^k_⋆, ω^k+1))|_β=0 ∀ k=1⋯ N-2,


then we have:


    λ_⋆^k = d_β(s_β^k )|_β=0  ∀ k = 1 ⋯ N - 1,


where (λ_k)_k=1 ⋯ N-1 are the Lagrangian multipliers associated to the multilevel optimization problem defined in Eq. (<ref>).


 We prove this result by backward induction on k. 


Initialization (k = N-1). By definition, s^N-1_β satisfies :


    β∇_s^N-1ℓ(ô_⋆, y) + ∇_1 E^N-1(s^N-1_β, θ^N-1, s^N-2_⋆, ω^N-1) = 0


Differentiating Eq. (<ref>) with respect to β and evaluating the resulting expression at β = 0, we obtain:


    ∇_s^N-1ℓ(ô_⋆, y) + ∇_1^2 E^N-1(s^N-1_⋆, θ^N-1, s^N-2_⋆, ω^N-1)· d_β s^N-1_β|_β=0 = 0


Substracting out Eq. (<ref>) defining the Lagrangian multiplier λ^N-1_⋆ and Eq. (<ref>), we obtain:


    ∇_1^2 E^N-1(s^N-1_⋆, θ^N-1, s^N-2_⋆, ω^N-1)·( d_β s^N-1_β|_β=0 - λ_⋆^N-1) = 0


By invertibility of ∇_1^2 E^N-1(s^N-1_⋆, θ^N-1, s^N-2_⋆, ω^N-1), we therefore have that:


    λ_⋆^N-1=d_β s^N-1_β|_β=0


Backward induction step (k + 1 → k). Let us assume that λ^k+1_⋆ = d_β s^k+1_β |_β=0. We want to prove that λ^k_⋆ = d_β s^k_β |_β=0. Again, s^k+1_β satisfies by definition:


    βδ s^k + ∇_1 E^k(s^k_β, θ^k, s^k-1_⋆, ω^k) = 0,   δ s^k := d_β.(∇_3 E^k+1(s^k+1_β, θ^k+1, s^k_⋆, ω^k+1))|_β=0.


On the one hand, proceeding as for the initialization step, differentiating Eq. (<ref>) with respect to β and taking β=0 yields:


    δ s^k + ∇_1^2 E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k)· d_β s^k_β |_β=0 = 0.


On the other hand, note that δ s^k rewrites :


    δ s^k    = d_β.(∇_3 E^k+1(s^k+1_β, θ^k+1, s^k_⋆, ω^k+1))|_β=0
       =∇_1, 3^2 E^k+1(s^k+1_⋆, θ^k+1, s^k_⋆, ω^k+1) ·.d s^k+1_β|_β=0
       = ∇_1, 3^2 E^k+1(s^k+1_⋆, θ^k+1, s^k_⋆, ω^k+1) ·λ^k+1_⋆,


where we used at the last step the recursion hypothesis. Therefore combining Eq. (<ref>) and Eq. (<ref>), we get:


    ∇_1, 3^2 E^k+1(s^k+1_⋆, θ^k+1, s^k_⋆, ω^k+1) ·λ^k+1_⋆ + ∇_1^2 E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k)· d_β s^k_β |_β=0 = 0.


Finally, we substract out Eq. (<ref>) and Eq. (<ref>) to obtain:


    ∇_1^2 E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k)·( d_β s^k_β |_β=0  - λ_⋆^k)= 0.


We conclude again by invertibility of ∇_1^2 E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k) that λ^k_⋆ = d_β s^k_β |_β=0.






Assuming a model of the form Eq. (<ref>), we denote s^1_⋆, x^1_⋆, ⋯, s^N-1_⋆, ô_⋆ the states computed during the forward pass as prescribed by Alg. <ref>. 
We define the nudged state of block k, denoted as s^k_β,  implicitly through ∇_1 ℱ^k(s^k_β, θ^k, x^k_⋆, δ s^k, β) = 0 with:

    ℱ^k(s^k, θ^k, x^k_⋆, δ s^k, β) := E^k(s^k, θ^k, x^k_⋆) + β s^k^⊤·δ s^k.


Denoting δ s^k and Δ x^k the error signals computed at the input of the feedforward block F^k and of the energy-based block E^k respectively, g_θ^k and g_ω^k the gradients of the loss function:


    ∀ k = 1, ⋯, N-1:  g_θ^k := d_θ^k𝒞,     ∀ k =1 ⋯ N:  g_ω^k := d_ω^k𝒞,



then the following chain rule applies:


    δ s^N-1 := ∇_s^N-1ℓ(ô_⋆, y),    g_ω^N = ∂_2 F^N(s^N-1_⋆, ω^N )^⊤·∇_1ℓ(ô_⋆, y) 
       ∀ k=1 ⋯ N-1: 
       {[ Δ x^k = .d_β(∇_3E^k (s^k_β, θ^k, x^k_⋆))|_β=0,    g_θ^k = .d_β(∇_2E^k (s^k_β, θ^k, x^k_⋆))|_β=0;              δ s^k-1 = ∂_1 F^k(s^k-1_⋆, ω^k )^⊤·Δ x^k,    g_ω^k = ∂_2 F^k(s^k-1_⋆, ω^k)^⊤·Δ x^k ].






Combining Lemma <ref> and Lemma <ref>, the following chain rule computes loss gradients correctly:


    δ s^N-1 := ∇_s^N-1ℓ(ô_⋆, y),    g_ω^N = ∂_2 F^N(s^N-1_⋆, ω^N )^⊤·∇_1ℓ(ô_⋆, y)
       ∀ k=1 ⋯ N-1: 
       {[ Δ s^k-1 = d_β.(∇_3 E^k(s^k_β, θ^k, s^k-1_⋆, ω^k))|_β=0,    g_θ^k = ∇^2_1, 2E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k) · d_β s^k_β |_β=0;                                                            g_ω^k = ∇_1, 4^2 E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k)· d_β s^k_β |_β=0 ].


Therefore to conclude the proof, we need to show that ∀ k = 1, ⋯, N-1:


    d_β.(∇_3 E^k(s^k_β, θ^k, s^k-1_⋆, ω^k))|_β=0   =  ∂_1 F^k(s^k-1_⋆, ω^k )^⊤·.d_β(∇_3E^k (s^k_β, θ^k, x^k_⋆))|_β=0
    ∇^2_1, 2E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k) · d_β s^k_β |_β=0   =.d_β(∇_2E^k (s^k_β, θ^k, x^k_⋆))|_β=0
    ∇_1, 4^2 E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k)· d_β s^k_β |_β=0   = ∂_2 F^k(s^k-1_⋆, ω^k)^⊤·.d_β(∇_3E^k (s^k_β, θ^k, x^k_⋆))|_β=0


Let k ∈ [1, N-1]. We prove Eq. (<ref>) as:


    d_β.(∇_3 E^k(s^k_β, θ^k, s^k-1_⋆, ω^k))|_β=0   =d_β.(∇_s^k-1 E^k(s^k_β, θ^k, F^k(s^k-1_⋆, ω^k)))|_β=0
       = ∂_1 F^k(s^k-1_⋆, ω^k)^⊤·.d_β(∇_3E^k (s^k_β, θ^k, x^k_⋆))|_β=0


Eq. (<ref>) can be obtained as:


    ∇^2_1, 2E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k) · d_β s^k_β |_β=0   = .d_β(∇_2 E^k(s^k_β, θ^k, s^k-1_⋆, ω^k) )|_β=0
       = .d_β(∇_2E^k (s^k_β, θ^k, x^k_⋆))|_β=0


Finally and similarly, we have:


    ∇_1, 4^2 E^k(s^k_⋆, θ^k, s^k-1_⋆, ω^k)· d_β s^k_β |_β=0   = . d_β( ∇_4E^k(s^k_β, θ^k, s^k-1_⋆, ω^k)) |_β=0
       =  . d_β( ∇_ω^k E^k(s^k_β, θ^k, F^k(s^k-1_⋆, ω^k))) |_β=0
       = . d_β(∂_2 F(s^k-1_⋆, ω^k )^⊤·∇_3 E^k(s^k_β, θ^k, x^k_⋆)) |_β=0
       = ∂_2 F(s^k-1_⋆, ω^k )^⊤·. d_β(∇_3 E^k(s^k_β, θ^k, x^k_⋆)) |_β=0





Under the same hypothesis as Theorem <ref> and Lemma <ref>, then the following chain rule applies to compute error signals backward from the output layer:

    {[                             δ s^N-1 := ∇_s^N-1ℓ(ô_⋆, y),    g_ω^N = ∇_ω^Nℓ(ô_⋆, y);                                                           Δ x^k = σ'(x^k_⋆) ⊙δ s^k; δ s^k-1 = ∂_1 F^k(s^k-1_⋆, ω^k )^⊤·Δ x^k,    g_ω^k = ∂_2 F^k(s^k-1_⋆, ω^k)^⊤·Δ x^k ].




Let k ∈ [1, N-1]. As we can directly apply Theorem <ref> here, proving the result simply boils down to showing that:


    Δ x^k = σ'(x^k_⋆) ⊙δ s^k


First, we notice that when E^k is of the form of Eq. (<ref>), then Δ x^k reads as:


    Δ x^k = .d_β(∇_3E^k (s^k_β, θ^k, x^k_⋆))|_β=0 = -. d_β(s^k_β) |_β=0.


s^k_β satisfies, by definition and when U^k = 0:


    σ^-1(s^k_β) - x^k_⋆ + βδ s^k = 0 
    ⇔      s^k_β = σ( x^k_⋆ - βδ s^k)


Combining Eq. (<ref>) and Eq. (<ref>)
yields Eq. (<ref>), and therefore, along with Theorem <ref>, the chain-rule Eq. (<ref>).


We showcase in Alg. <ref> the resulting algorithm with finite β and implicit BP-EP chaining, with lines in blue highlighting differences with the general algorithm Alg. <ref>. 






 §.§ Explicit BP-EP chaining



We presented in Alg. <ref> a “pure” EP algorithm where the BP-EP gradient chaining is implicit. We show below, inside Alg. <ref>, an alternative implementation (equivalent in the limit β→ 0) where the use of BP through feedforward modules is explicit and which is the direct implementation of Theorem <ref>. We also show the resulting algorithm when the ff-EBM reduces to a feedforward net (Lemma <ref>) inside Alg. <ref>, highlight in blue the statements which differ from the general case presented inside Alg. <ref>.









 §.§ Model and algorithm details





  
Equilibrium computation. As mentioned in Section <ref>, the energy function of the k^ th  EB block has the form:


    E^k(s^k, θ^k, x^k) := G^k(s^k) - s^k^⊤· x^k + U^k(s^k, θ^k)
        ,


where x^k is the output of the preceding feedforward block. For a given choice of a continuously invertible activation function, G^k_σ is defined as:


    G^k_σ(s^k) :=∑_i=1^(s^k)∫^s_iσ_i^-1(u_i)du_i     ∇ G^k_σ(s^k)_i = σ^-1_i(s^k_i)   ∀ i = 1 ⋯(s^k).


To be more explicit and as we did previously, we re-write the augmented energy-function which encompasses both the k^ th EB block and the feedforward module that precedes it:


    E^k(s^k, θ^k, s^k-1_⋆, ω^k) := E^k(s^k, θ^k, F^k(s^k-1_⋆, ω^k))
        
        .


We showed that when is chosen such that ∇ G = σ^-1 for some activation function σ, then the steady state of the k^ th block reads:


    s^k_⋆ := σ(x^k - ∇_1 U^k(s^k_⋆, θ^k) ),


which justifies the following fixed-point iteration scheme, when the block is influenced by some error signal δ s with nudging strength β:


    s^k_±β, t + 1σ(x^k - ∇_1 U^k(s^k_±β, t, θ^k) ∓βδ s^k ).


The dynamics prescribed by Eq. <ref> are also used for the inference phase with β=0. To further refine Eq. (<ref>), let us re-write Eq. (<ref>) with a layer index ℓ where ℓ∈ [1, L_k] with L_k being the number of layers in the k^ th block, and replacing x^k by its explicit expression:


    ∀ℓ = 1 ⋯ L_k:     s^k_ℓ, ±β, t + 1σ(F^k(s^k-1_⋆, ω^k-1) - ∇_s^k_ℓ U^k(s^k_±β, t, θ^k) ∓βδ s^k ).


As done in past EP works <cit.> and for notational convenience, we introduce the primitive function of the k^ th block as:


    Φ^k(s^k, θ^k, s^k-1_⋆, ω^k) := s^k^⊤· F^k(s^k-1_⋆, ω^k) - U^k(s^k, θ^k)


such that Eq. (<ref>) re-writes:


    ∀ℓ = 1 ⋯ L_k:  s^k_ℓ, ±β, t + 1σ(∇_s^k_ℓΦ(s^k_±β, t, θ^k, s^k-1_⋆, ω^k) ∓βδ s^k ).


Eq. (<ref>) depicts a synchronous scheme where all layers are simultaneously updated at each timestep. Another possible scheme, employed by <cit.>, instead prescribes to asynchronously update odd and even layers and was shown to speed up convergence in practice:


    {[                                                         ∀  ℓ∈{1, ⋯, L_k}:     s^k_ℓ, ±β, t + 1/2σ(∇_s^k_ℓΦ(s^k_±β, t, θ^k, s^k-1_⋆, ω^k) ∓βδ s^k ),;                                                         ∀  ℓ∈{1, ⋯, L_k}: s^k_ℓ, ±β, t + 1σ(∇_s^k_ℓΦ(s^k_±β, t + 1/2, θ^k, s^k-1_⋆, ω^k) ∓βδ s^k ). ].


We formally depict this procedure as the subroutine  inside Alg. <ref>. In practice, we observe that it was more practical to use a fixed number of iterations rather than using a convergence criterion with a fixed threshold.






  
Nudging the last block. From looking at the procedure prescribed by Theorem <ref> and algorithms thereof (Alg. <ref>, Alg. <ref>), all the error signals used to nudge the EB blocks are stationary, including the top-most block where the loss error signal is fed in. Namely, the augmented energy function of the last block reads as:


    ℱ^N-1(s^N-1, θ^N-1, x^N-1_⋆, β) := E^N-1(s^N-1, θ^N-1, x^N-1_⋆) + β s^N-1^⊤·∇_s^N-1ℓ(ô_⋆, y),


where ô_⋆ := F^N(s^N-1_⋆, ω^N) is constant. Up to a constant, Eq. (<ref>) uses the cost function linearized around s^N-1_⋆ instead of the cost function itself. This is, however, in contrast with most EP implementations where the nudging force acting upon the EB block is usually elastic, i.e. the nudging depends on the current state of the EB block. More precisely, instead of using Eq. (<ref>), we instead use:


    ℱ^N-1(s^N-1, θ^N-1, x^N-1_⋆, β) := E^N-1(s^N-1, θ^N-1, x^N-1_⋆) + βℓ(F^N(s^N-1, ω^N), y),


This results in the following asynchronous fixed-point dynamics for the last block:


    {[                                                                        ∀  ℓ∈{1, ⋯, L_k}:     s^k_ℓ, ±β, t + 1/2σ(∇_s^k_ℓΦ(s^k_±β, t, θ^k, s^k-1_⋆, ω^k) ∓β∇_s^kℓ(s^k_±β, t, y) ),;                                                                        ∀  ℓ∈{1, ⋯, L_k}: s^k_ℓ, ±β, t + 1σ(∇_s^k_ℓΦ(s^k_±β, t + 1/2, θ^k, s^k-1_⋆, ω^k) ∓β∇_s^kℓ(s^k_±β, t, y) ). ].


The resulting  subroutine, applying for the last block, is depicted inside Alg. <ref>.







  
Readout. <cit.> introduced the idea of the “readout” whereby the last linear layer computing the loss logits is not part of the EB free block dynamics but simply “reads out” the state of the penultimate block. In all our experiments we use such a readout in combination with the cross entropy loss function. Using our formalism, our readout is simply the last feedforward transformation used inside ℓ, namely F^N(·, ω^N).



  
Deep Hopfield Networks (DHNs). In our experiments, we used weight matrices of the form:


    θ^k = [       0 θ_1^k^⊤       0                ;   θ_1^k       0 θ_2^k^⊤                ;       0   θ_2^k       ⋱       ⋱        ;                       ⋱       0 θ_L^k^⊤;                           θ_L^k       0 ],


whereby each layer ℓ is only connected to its adjacent neighbors. Therefore, fully connected and convolutional DHNs with L layers have an energy function of the form:


    U^k_ FC(s^k, θ^k)    := -1/2s^k⊤·θ^k · s^k = -1/2∑_ℓ s^k^⊤_ℓ + 1·θ_ℓ^k · s^k_ℓ
    
        U^k_ CONV(s^k, θ^k)    := -1/2s^k ∙(θ^k ⋆ s^k)  = -1/2∑_ℓ s^k_ℓ + 1∙( θ_ℓ^k ⋆ s^k_ℓ)




  
Detailed inference algorithm. With the aforementioned details in hand, we re-write the inference algorithm Alg. <ref> presented in the main as a  subroutine. 







  
Detailed implicit EP-BP chaining algorithm.
We provide a detailed implementation of our algorithm presented in the main (Alg. <ref>) in Alg. <ref>. As usually done for EP experiments, we always perform a “free phase” to initialize the block states ( subroutine, Alg. <ref>). Then, two nudged phases are applied to the last block and parameter gradients subsequently computed, as done classically ( subroutine for the last block, Alg. <ref>), with an extra computation to compute the error current to be applied to the penultimate block (δ s^N-2). Then, the same procedure is recursed backward through blocks (Alg. <ref>), until reaching first block.













  
Details about the implicit differentiation baseline.
We describe our implementation of Implicit Differentiation (ID) inside Alg. <ref>. First, we relax all blocks sequentially to equilibrium following Alg. <ref> and we do not track gradients throughout this first phase, using T_ free fixed-point iteration steps per block. Then, initializing the block states with those computed at the previous step, we re-execute the same procedure (still with Alg. <ref>), this time tracking gradients and using T_ nudge steps fixed-point iteration steps for each block. Then, we use automatic differentiation to backpropagate through the last T_ nudge steps for each block, namely backpropagating, backward in time, through equilibrium. 







 §.§ Details about the static gradient analysis



  
Single block case.  Let us consider a single block with parameters θ and state s. Earlier, we showed that the (free) dynamics of a block read in the general form as:


    s_t + 1∇_1 Φ(s_t, θ, x),


where we have ignored the activation function, as done by <cit.>. Therefore, the computational graph for Eq. (<ref>) is that of a recurrent neural network with a static input (x) where the parameters θ are shared across time steps. Therefore, one could view Eq. (<ref>) as:


    s_t + 1∇_1 Φ(s_t, θ(t) = θ, x ).


At the end of the recurrent dynamics Eq. (<ref>), a loss is computed based on the output of the last block, namely:


    𝒞 = ℓ(s_T, y),


where s_T is the state of the last block after T iterations of the recurrent dynamics Eq. (<ref>). In principle, the gradient of the loss 𝒞 with respect to θ reads as:


    g_θ = d_θ(1)𝒞 + d_θ(2)𝒞 + ⋯ + d_θ(T)𝒞.


In other words,  the gradient of the loss 𝒞 with respect to θ is the sum of all “sensitivities” of the loss to the same parameter taken at all timesteps where it is involved. One can define the truncated gradient at T-t, i.e. t steps backward in time from T as:


    ĝ^ AD_θ(t) := ∑_k=0^t-1 d_θ(T - k)𝒞,


such that ĝ^ AD_θ(T) = g_θ. AD stands for “automatic differentiation” because this is how these gradients may be computed in practice – it can also be viewed as an instantiation of backpropagation through time (BPTT) with a static input. Then, writing the nudged dynamics prescribed by EP as:


    s_t + 1^±β := ∇_1 Φ(s^±β_t, θ, x ) ∓β∇_1 ℓ(s_t^β, y),


we define the truncated EP gradient estimate at timestep t through the second phase as:


    ĝ^ EP_θ(β, t) := 1/2β(∇_2 Φ(s^β_t, θ, x) -  ∇_2 Φ(s^-β_t, θ, x)).


Furthermore, we assume that s is already at its steady state s_⋆ over the last K steps, i.e. T-K, T-K-1, ..., T, so that AD here occurs through equilibrium, then the following equality holds <cit.>:


    ∀ t = 1 ⋯ K:   lim_β→ 0ĝ^ EP_θ(β, t) = ĝ^ AD_θ(t)


This property was called the Gradient Descending Updates (GDU) property. In practice, the GDU property is a tremendously helpful sanity check to debug code to make sure EP gradients are correctly computed. Our coding work heavily relied on this property.



  
Multiple blocks. In ff-EBMs, the GDU property may immediately hold for the last block. For penultimate block backwards, the sole difference compared to the previous paragraph is simply the nudging force that is applied to the model, for which Eq. (<ref>) may simply be changed into:



    s_t + 1^±β∇_1 Φ(s^±β_t, θ, x ) ∓βδ s,


where δ s is a constant error signal. Therefore, so long as the error signal δ s computed by EP matches that computed by AD (which we proved inside Theorem <ref>), the GDU property should still hold.












 §.§ Experimental Details



  §.§.§ Datasets


Simulations were run on CIFAR-10, CIFAR-100 and Imagenet32 datasets, all consisting of color images of size  32 × 32 pixels. CIFAR-10  <cit.> includes 60,000 color images of objects and animals . Images are split into 10 classes, with 6,000 images per class. Training data and test data include 50,000 images, and 10,000 images respectively. Cifar-100 <cit.> likewise comprises 60,000, and features a diverse set of objects and animals split into 100 distinct classes. Each class contains 600 images. Like CIFAR-10, the dataset is divided into a training set with 50,000 images and a test set containing the remaining 10,000 images. The ImageNet32 dataset <cit.> is a downsampled version of the original ImageNet dataset <cit.> containing 1,000 classes with 1,281,167 training images, 50,000 validation images,  100,000 test images and 1000 classes.



  §.§.§ Data Pre-processing

All data were normalized according to statistics shown in <ref>and augmented with 50 % random horizontal flips. Images were randomly cropped and padded with the last value along the edge of the image. 





  §.§.§ Simulation Details



  
Weight Initialization
Equilibrium propagation, similar to other machine learning paradigms reliant on fixed-point iteration <cit.>, is highly sensitive to initialization statistics <cit.>, hence conventionally difficult to tune, and requiring many iterations for the three relaxation phases. Initialization of weights as Gaussian Orthogonal Ensembles (GOE) ensures better stability  (reduced variance) and, combined with other stabilizing measures discussed below, empirically yields faster convergence. 


According to GOE, weights are intialized as:



    W_ij∼𝒩(0, V/N),    if  i ≠ j 
    𝒩(0, 2V/N),    if  i = j


where 𝒩(μ, σ^2) denotes a Gaussian (normal) distribution with mean μ and variance σ^2.




  
State Initialization
All states are initialized as zero matrices.


  
Activation Functions
An important detail for faithful reproduction of these experiments is the choice and placement of activation functions applied during the iterative root-finding procedure. In the literature, clamping is conventionally applied at each layer, with the exception of the final layer, where it is sometimes included e.g. <cit.>, and sometimes omitted <cit.>. For these experiments we use both the standard hard activation employed in <cit.> and <cit.>, and the more conservative one given in <cit.>. Details are given in <ref>. 


In practice, we find that the laborieux activation, in conjunction with GOE, and the omission of clamping at the output of each block significantly enhances gradient stability and speeds convergence in the multiple block setting. In the interest of multi-scale uniformity and consistency with previous literature <cit.>  <cit.>, however, we apply clamping using the ernoult activation on all layers in our 6-layer architecture. 

 For the scaling experiments, we apply laborieux activation at every layer except the output of each block.  For the 12-layer splitting experiment, we do the same, omitting clamping from the output of the final layer of each block in the block-size=4 and block-size=3 experiments. However, in the block-size=2 case we clamp the output of the second and fourth blocks to preserve dynamics of the block-size=4 split. Such consistency is not possible for the block-size=3 experiment, constituting a possible discrepancy in the scaling dynamics.




  
Cross-entropy Loss and Softmax Readout

Following <cit.>, all models were implemented such that the output y is removed from the system (e.g. not included in the relaxation dynamics) but is instead the function of a weight matrix: W_out of size dim(y) ×dim(s), where s is the state of the final layer. For each time-step t, ŷ_t = softmax(W_outs_t).

The cross-entropy cost function associated with the softmax readout is then:


    l(s, y, W_out) = -∑_c=1^C y_c log(softmax_c(W_out· s)).


 It is important to note that by convention we refer to architectures throughout this text to the exclusion of the softmax readout, which is technically an additional layer, though not involved in the relaxation process.




  
Architecture
Details of architectures used in the experiments reported in <ref>  and <ref> are given in table <ref>. All convolutional layers used in experiments are of kernel size 3 and stride and padding 1. Max-pooling was applied with a window of 2 × 2 and stride of 2. For the 6-layer model used in the <ref> batch norm was applied after the first layer convolution and pooling operation. All other models in both experiments use batch-normalization on the first layer of each block after convolution and pooling (where applied). These details exclude the Linear Softmax readout of size n classes.






  
Hyperparameters
Hyperparameters and implementation details for <ref> and <ref> are given in table <ref>. Note that all architectural details for the 12-layer models are identical across splitting and scaling experiments. Additionally, identical hyperparameters were used for Cifar-100 and Imagenet experiments of <ref>. Unlike previous literature, the use of GOE intialization eliminates the need for separate layerwise learning rates and initialization parameters. One noteworthy detail is that only 100 epochs were used for the larger model for <ref> compared with 200 epochs for the smaller 12-layer model. This was due to prohibitively long run-time of training the larger model. Noteably performance still significantly improves with decreased overall runtime.




  
Other details
All experiments were run using Adam optimizer <cit.>and Cosine Annealing scheduler<cit.> with minimum learning rates indicated in <ref> and maximum T equal to epochs (i.e. no warm restarts). Code was implemented in Pytorch 2.0 and all simulations were run on A-100 SXM4 40GB GPUs.






§ NEURIPS PAPER CHECKLIST
























 







  * Claims
    
   Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
    
   Answer: 
    
    
    
    
    
    
    
    


  * Limitations
    
   Question: Does the paper discuss the limitations of the work performed by the authors?
    
   Answer:  
    
   Justification: we have a dedicated paragraph in the “Discussion” section of the paper which explicitly mentions limitations and future work.
    
    
    
    
    
    
    
    
    
    
    


  * Theory Assumptions and Proofs
    
   Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?
    
   Answer: 
    
   Justification: All proofs are included in the appendix.
    
    
    
    
    
    
    
    
    

    
  * Experimental Result Reproducibility
    
   Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
    
   Answer: 
     
   Justification: detailed information to reproduce experiment is provided in the appendix.
    
    
    
    
    
    
    
    
    
    
    
    
    
    



  * Open access to data and code
    
   Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?
    
   Answer:  
    
   Justification: We only employed public datasets and will release our code at the end of the reviewing process. In the meanwhile, we provided very detailed pseudo-algorithms in the appendix.
    
    
    
    
    
    
    
    
    
    
    



  * Experimental Setting/Details
    
   Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
    
   Answer: 
     
   Justification: all these details are provided in the appendix. 
    
    
    
    
    
    


  * Experiment Statistical Significance
    
   Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
    
   Answer:  
    
   Justification: we perform each of our training simulations on 3 different seeds and reported mean and standard deviation of the resulting performance. 
    
    
    
    
    
    
    
    
    
    
    
    


  * Experiments Compute Resources
    
   Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
    
   Answer:  
    
   Justification: this information is also inside our appendix.
    
    
    
    
    
    
    
    

  * Code Of Ethics
    
   Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics <https://neurips.cc/public/EthicsGuidelines>?
    
   Answer:  
    
    
    
    
    
    
    



  * Broader Impacts
    
   Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
    
   Answer:  
    
   Justification: we wrote a dedicated paragraph inside our “Discussion” section.
    
    
    
    
    
    
    
    
    
    

  * Safeguards
    
   Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
    
   Answer:  
    
   Justification: Our work poses no such risk at present as it only provides proof-of-concepts for systems which do not yet exist.
    
    
    
    
    
    
    


  * Licenses for existing assets
    
   Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
    
   Answer:  
    
   Justification: this work does not use existing assets.
    
    
    
    
    
    
    
    
    
    
    


  * New Assets
    
   Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
    
   Answer:  
 
   Justification: this paper does not release new assets.
    
    
    
    
    
    
    


  * Crowdsourcing and Research with Human Subjects
    
   Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)? 
    
   Answer:  
    
   Justification: our paper does not involve crowdsourcing nor research with human subjects.
    
    
    
    
    
    


  * Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
    
   Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
    
   Answer:  
    
   Justification: our paper does not involve crowdsourcing nor research with human subjects.
    
    
    
    
    
    
    





