Casting hybrid digital-analog training into hierarchical energy-based learning

Venue: NeurIPS_cc_2024_Workshop_MLNCP
Published: 2024-09-11 19:24:12

Abstract: Deep learning requires new approaches to combat the rising cost of training large models. The combination of energy-based analog circuits and the Equilibrium Propagation (EP) algorithm offers one compelling alternative to BP for gradient-based optimization of neural nets. In this work, we introduce a hybrid framework comprising feedforward (FF) and energy-based (EB) blocks housed on digital and analog circuits. We derive a novel algorithm to compute gradients end-to-end via BP \emph{and} EP, through FF and EB parts respectively, enabling EP to be applied to much more flexible and realistic architectures as  analog units are incorporated into digital circuitry over time. We  demonstrate the effectiveness of the proposed approach, showing that a standard Hopfield Network can be split into any shape while maintaining automatic differentiation performance. We apply it to ImageNet32 where we establish new SOTA in the EP and BP-alternative literature (46\% top-1).

Main paper: 
 introduction: 
 Deep learning today relies on three factors: i) GPUs, ii) feedforward (FF) models and iii) backprop (BP) . With skyrocketing demands of AI compute, exploration of new compute paradigms has become an economic, and environmental need [Strubell et al., 2020] . One path towards this goal is analog inmemory computing [Sebastian et al., 2020] , promising constant time complexity as well as significant reduction in energy consumption [Cosemans et al., 2019] . By mapping a neural network onto a physical system formalized as a "self-learning" energy-based (EB) analog circuit [Kendall et al., 2020 , Stern et al., 2023 , Dillavou et al., 2023 , Scellier, 2024] loss gradients can be derived via two physical relaxations to equilibrium [Scellier et al., 2024] . Equilibrium propagation (EP) [Scellier and Bengio, 2017] , is a suitable such EB-algorithm notable for its strong theoretical guarantees, relative scalability (in the realm of BP-alternatives) [Laborieux and Zenke, 2022, 2023] and for being 10, 000× more energy-efficient and faster on analog devices than the GPU counterpart [Yi et al., 2023] . Nevertheless, end-to-end EP-training is confronted with significant challenges when executed on analog hardware, including non-ideal physical behaviors which affect both the inference pathway [Wang et al., 2023 , Ambrogio et al., 2023] and parameter optimization [Nandakumar et al., 2020 , Spoon et al., 2021 , Lammie et al., 2024] , as well as incompatibility of common non-stationary operations with analog hardware such as activation functions, batchnorm, attention [Spoon et al., 2021 , Jain et al., 2022 , Liu et al., 2023 , Li et al., 2023] . One possibility is hybrid systems incorporating both FF and EB components. Indeed, the design of inferential engines made of analog and digital parts is nearing commercial maturity [Ambrogio et al., 2023 ], yet in-situ learning of such systems remains unexplored.Here, we propose a theoretical framework to extend end-to-end gradient computation to a setting where the system at use may or may not be fully analog. Our work contends that a combination Submitted to the Second Workshop on Machine Learning with New Compute Paradigms at NeurIPS (MLNCP 2024) . Do not distribute.of digital and analog hardware, with FF and EB parts trained via BP and EP respectively, can leverage advances from both digital and analog hardware in the near-term. Specifically, we explore Feedforward-tied Energy-based Models (ff-EBMs, Section 3.1) whose inference pathway is a composition of FF and EB modules (Alg. 1). We show that gradients in ff-EBMs can be computed in an end-to-end fashion (Section 3.3), via BP in FF blocks and EP in EB blocks (Theorem 3.1, Alg. 2) (Section 3.2), and that when each analog block comprises a single layer, the ff-EBM is purely FF (Lemma A.2), whose algorithm is equivalent to BP (Corollary A.5.1). Finally, we demonstrate the effectiveness of our algorithm on ff-EBMs where EBM blocks are DHNs (Section 4). We show that i) gradient estimates computed by our algorithm (Alg. 2) near perfectly match gradients computed by end-to-end automatic differentiation (Fig. 2 ), ii) a standard DHN model can be split into a ff-EBM with equivalent layers and architecture without compromising performance and remaining on par with automatic differentiation, with significant reductions in run-time due to the use of smaller EB blocks (Section 4.2), iii) the proposed approach scales, yielding 46 % top-1 (70% top-5) validation accuracy on ImageNet32, beating the current SOTA for BP alternatives by a large margin. 
background: 
 Notations. Given a differentiable mapping A : R n → R m , we denote its total derivative wrt s j as d sj A(s) := dA(s)/ds j ∈ R m , its partial derivative wrt s j as ∂ j A(s) := ∂A(s)/∂s j ∈ R m . WhenA takes scalar values (m = 1), its gradient wrt s j is denoted as ∇ j A(s) := ∂ j A(s) ⊤ .For a given static input and set of weights, 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) ∈ R where x, s, and θ respectively denote a static input, hidden and output neuron states, and model parameters (weights). Each such tuple defines a configuration with an associated scalar energy value. Amongst all configurations for an input x and some model parameters θ, the model prediction s ⋆ is given as an equilibrium state which minimizes the energy function:s ⋆ := arg min s E(s, θ, x).(1)Assuming that ∇ 2 s E(x, s ⋆ , θ) is invertible, note that the equilibrium state s ⋆ implicitly depends on x and θ by virtue of the implicit function theorem [Dontchev et al., 2009] . Therefore our goal when training an EBM is to adjust the model parameters θ such that s ⋆ (x, θ) minimizes some cost function ℓ : s, y → ℓ(s, y) ∈ R where y is some ground-truth associated to x. More formally, this objective can be stated with the following bilevel optimization problem [Zucchet and Sacramento, 2022] :min θ C(x, θ, y) := ℓ(s ⋆ , y) s.t. s ⋆ = arg min s E(s, θ, x).Solving Eq. ( 2) in practice amounts to computing the gradient of its outer objective C(x, θ, y) with respect to θ (d θ C(x, θ, y)) and then perform gradient descent over θ.An algorithm used to train an EBM model in the sense of Eq. ( 2) may be called an EBL algorithm [Scellier et al., 2024] . EP [Scellier and Bengio, 2017] is an EBL algorithm which computes an estimate of d θ C(x, θ, y) with at least two phases. In its first phase, the model evolves freely to s ⋆ = arg min s E(s, θ, x). Then, the model is slightly nudged towards decreasing values of cost ℓ and settles to a second equilibrium state s β . In practice, we augment the energy function E by an additional term βℓ(s, y) where β ∈ R ⋆ is called the nudging factor. Weights are then updated to increase the energy of s ⋆ and decrease that of s β , thereby "contrasting" these two states. More formally, Scellier and Bengio [2017] prescribe in the seminal EP paper:EQUATIONwhere α denotes some learning rate. EP comes in different flavours depending on the form of β inside Eq. ( 2). Centered EP (C-EP), where two nudged states of opposite nudging strengths (±β) are contrasted, is shown to work best in practice [Laborieux et al., 2021 , Scellier et al., 2024] and reads:EQUATION3 Tying energy-based models with feedforward blocksHere we introduce our model, its associated optimization problem and a new learning algorithm. We show how learning amounts to solving a multi-level optimization problem (Section 3.2), and propose a BP-EP gradient chaining algorithm as a solution (Section 3.3, Theorem 3.1, Alg. 2). We highlight that ff-EBMs reduce to standard FF nets (Lemma A.2) and the proposed BP-EP gradient chaining algorithm to standard BP (Corollary A.5.1) when each EB block comprises a single hidden layer.3.1 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 FF transformation F 1 parameterized by some weights ω 1 , yielding an output x 1 ⋆ . 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 FF transformation F 1 with weights ω 1 and the above procedure repeats until reaching the output layer ô.Algorithm 1 ff-EBM inference1: s ← x 2: for k = 1 • • • N -1 do 3:x ← F k s, ω k Form of the energy functions. Given the k th EB block of a ff-EBM, the associated energy function E k takes some static input x k from the output of the preceding FF transformation, has hidden neurons s k and is parameterized by weights θ k . More precisely:E k (s k , θ k , x k ) := G k (s k ) -s k ⊤ • x k + U k (s k , θ k ) (4)Eq. ( 4) reveals three contributions to the energy. The first determines the non-linearity applied inside the EB block, the second accounts for a purely FF contribution from the previous FF block F k , while the third accounts for internal interactions within the layers of the EB block. An interesting edge case highlighted by Eq. ( 11) is when U k = 0 for all k's, i.e. no intra-block layer interactions, i.e. the EB block comprises a single layer. 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 A.2 inside Appendix A.1. 
multi-level optimization of ff-ebms: 
 In the same way that learning EBMs can naturally be cast as a bilevel optimization problem, learning ff-EBMs is a multi-level optimization problem where variables being optimized in the inner subproblems are 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. ( 4) to highlight the dependence between two consecutive EB block states.Namely, by writingE k (s k , θ k , s k-1 ⋆ , ω k ) := E k s k , θ k , F k s k-1 ⋆, ω k-1 , it can be seen 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. Therefore, the learning problem for a ff-EBM writes, omitting parameter dependence of the loss and the energy functions for readability:min θ ℓ(s N -1 ⋆ , y) s.t. s N -1 ⋆ = arg min s E N -1 (s, s N -2 ⋆ ), • • • , s.t. s 1 ⋆ = arg min s E 1 (s, x) (5)Main result: explicit BP-EP chaining. Based on the multilevel optimization problem in Eq. ( 5) ,we state the main theoretical result of this paper in Theorem 3.1 (see proof in Appendix A.2) and the resulting algorithm (Alg. 2), which also comes in another flavor (see Appendix A.3).Theorem 3.1 (Informal). Assuming a ff-EBM as depicted in Alg. 1, we denotes 1 ⋆ , x 1 ⋆ , • • • , s N -1 ⋆ , ô⋆the states computed during the forward pass as depicted in Alg. 1. We define the nudged state of block k, denoted ass k β , implicitly through ∇ 1 F k (s k β , θ k , x k ⋆ , δs k , β) = 0 with: F k (s k , θ k , x k ⋆ , δs k , β) := E k (s k , θ k , x k ⋆ ) + βs k ⊤ • δs k (6)Denoting δs k and ∆x k the error signals computed at the input of the FF 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) (7) ∀k = 2 • • • N -1 :    ∆x k = d β ∇ 3 E k (s k β , θ k , x k ⋆ ) β=0 , g θ k = d β ∇ 2 E 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 (8)Algorithm 2 Implicit BP-EP gradient chaining (Theorem (3.1))1: δs, g ω N ← ∇ s N -1 ℓ(ô ⋆ , y), ∇ ω N ℓ(ô ⋆ , y) ▷ Single backprop step 2: for k = N -1 • • • 1 do 3: s β ← Optim s E k (s, θ k , s k-1 ⋆ , ω k ) + βs ⊤ • δs ▷ EP through E k 4: s -β ← Optim s E k (s, θ k , s k-1 ⋆ , ω k ) -βs ⊤ • δs 5: g θ k ← 1 2β ∇ 2 E k (s β , θ k , s k-1 ⋆ , ω k ) -∇ 2 E k (s -β , θ k , s k-1 ⋆ , ω k ) 6: g ω k ← 1 2β ∇ 4 E k (s β , θ k , s k-1 ⋆ , ω k ) -∇ 4 E k (s -β , θ k , s k-1 ⋆ , ω k ) ▷ i-BP through F k 7: δs ← 1 2β ∇ 3 E k (s β , θ k , s k-1 ⋆ , ω k ) -∇ 3 E k (s -β , θ k , s k-1 ⋆ , ω k ) 8: end for 4 ExperimentsIn this section, we present the ff-EBMs used in our experiments (Section 4.1) and conduct a static gradient analysis, where we observe that the gradients computed by our algorithm are near perfectly aligned with those computed by automatic differentiation (Fig. 2 ). We then show on CIFAR-10 that performance of ff-EBMs does not degrade on block splits of decreasing size (despite much lower wall clock time) and remains on par with automatic differentiation (Section 4.2). 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 4.3). 
setup & static gradient analysis: 
 Model. Using the same notations as in Eq. ( 4), the ff-EBMs at use in this section are defined through:U k FC (s k , θ k ) := -1 2 s k⊤ • θ k • s k , U k CONV (s k , θ k ) := -1 2 s k • θ k ⋆ s k , F k (s k-1 , ω k ) := BN P ω k CONV ⋆ s k-1 L ; ω k α , ω k β (9) with BN(•; ω k α ω k β), P and ⋆ the batchnorm, pooling and convolution operations, • the generalized dot product for tensors and s k := s k ⊤ 1 , • • • s k ⊤ L ⊤ theFor a given EBM and a fixed number layers, we explore two models of depth (L = 6 and L = 12)with various block splits maintaining an equivalent architecture (e.g. for L = 6, 1 block of 6 layers, 2 blocks of 3 layers, etc.). Two variants of this experiment were run, one employing a tolerance-based (TOL) approach to relaxations to equilibrium, and the other utilizing a fixed number of iterations in all splits. Results of CIFAR-10 with a TOL-based scheme are shown in Fig 3 . We observe that performance achieved by EP on 6-layer EBM is improved with overall reduction in WCT.Furthermore, when using fixed iterations (see Table 6 ),performance is maintained across 4 distinct splits between 89% and 90% and is on par with both the ID baseline for each, and with the literature (in EBMs of roughly equivalent size) [Scellier et al., 2024, Laborieux and Zenke, 2022] . The same trend is on ff-EBMs with L = 12 where an accuracy of 92% across 3 splits, matches ID, and (in the fixed iteration setting) surpasses EP state-of-the art on CIFAR-10 [ Scellier et al., 2024] . The significant reduction in WCT observed is due to the fact that by design training time for ff-EBMs scales linearly with number of blocks rather than supralinearly with number of layers.EP ID Top-1 (%) WCT Top-1 (%) WCT L =6bs=6 88.5 ±0.2 13:30 88.3 ±0.2 13:25 bs=3 88.7 ±0.2 6:01 88.7 ±0.2 5:54 bs=2 90.1 ±0.2 3:10 90.1 ±0.1 3:09EQUATION88.3 ±0.1 11:00 88.5 Finally, we consider ff-EBMs of fixed block size 2 and train them with two different depths (L = 12and L = 15) on CIFAR-100 and ImageNet32 by EP and ID 1. We observe that EP matches ID performance on all models and tasks, and the performance obtained by training the 15-layer models by exceeds ImageNet32 SOTA for EP by 10% [Laborieux and Zenke, 2022] and by around 5% among all BP-alternatives [Høier et al., 2023] . Algorithm baseline. As an algorithmic baseline, we simply use automatic differentiation (AD) backward through the fixed-point iteration scheme Eq. ( 56) with β = 0 and directly initializing s k t=0 = s ⋆ . This version of AD, where we backpropagate through equilibrium, is known as "RecurrentBackpropagation" [Almeida, 1987 , Pineda, 1987] or Implicit Differentiation (ID).Equilibrium computation. As mentioned in Section 3.1, the energy function of the k th EB block has the form:EQUATIONwhere 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:EQUATION)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:EQUATION)We showed that when is chosen such that ∇G = σ -1 for some activation function σ, then the steady state of the k th block reads:EQUATIONwhich justifies the following fixed-point iteration scheme, when the block is influenced by some error signal δs with nudging strength β:EQUATION)The dynamics prescribed by Eq. 52 are also used for the inference phase with β = 0. To further refine Eq. ( 52), let us re-write Eq. ( 52) 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:EQUATION)As done in past EP works [Ernoult et al., 2019 , Laborieux et al., 2021 , Laborieux and Zenke, 2022 , 2023 , Scellier et al., 2024] 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 ) (54)such that Eq. ( 53) re-writes:EQUATION)Eq. ( 55) depicts a synchronous scheme where all layers are simultaneously updated at each timestep.Another possible scheme, employed by Scellier et al. [2024] , instead prescribes to asynchronously update odd and even layers and was shown to speed up convergence in practice:EQUATION)We formally depict this procedure as the subroutine Asynchronous inside Alg. 6. 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.Algorithm 6 Asynchronous (for all blocks until penultimate)Input: T , θ k , ω k , s k-1 ⋆ , β, δs k Output: s k β 1: s k ← 0 2: for t = 1 • • • T do 3: ∀ odd ℓ ∈ {1, • • • , L k } : s k ℓ,β ← σ ∇ s k ℓ Φ s k β , θ k , s k-1 ⋆ , ω k -βδs k 4: ∀ even ℓ ∈ {1, • • • , L k } : s k ℓ,β ← σ ∇ s k ℓ Φ s k β , θ k , s k-1 ⋆ , ω k -βδs k 5: end forNudging the last block. From looking at the procedure prescribed by Theorem 3.1 and algorithms thereof (Alg. 2, Alg. 4), 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:EQUATIONwhereô⋆ := F N s N -1 ⋆ , ω N is constant.Up to a constant, Eq. ( 58) 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. ( 57), we instead use:EQUATION)This results in the following asynchronous fixed-point dynamics for the last block:   ∀ odd ℓ ∈ {1, • • • , L k } : s k ℓ,±β,t+ 1 2 ← σ ∇ s k ℓ Φ s k ±β,t , θ k , s k-1 ⋆ , ω k ∓ β∇ s k ℓ(s k ±β,t , y) , ∀ even ℓ ∈ {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 Asynchronous subroutine, applying for the last block, is depicted inside Alg. 7.Readout. Laborieux et al. [2021] 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:EQUATIONAlgorithm 7 Asynchronous (for last block)Input: T , θ N -1 , ω N -1 , ω N , s k-1 ⋆ , β, ℓ (cost function), y Output: s N -1 β 1: s N -1 ← 0 2: for t = 1 • • • T do 3: ∀ odd ℓ ∈ {1, • • • , L N }: 4: s N -1 ℓ,β ← σ ∇ s N -1 ℓ Φ s N -1 β , θ N -1 , s N -2 ⋆ , ω N -1 -β∇ s N -1 ℓ ℓ(F N (s N -1 , ω N ), y) 5: ∀ even ℓ ∈ {1, • • • , L N }: 6: s N -1 ℓ,β ← σ ∇ s N -1 ℓ Φ s N -1 β , θ N -1 , s N -2 ⋆ , ω N -1 -β∇ s N -1 ℓ ℓ(F N (s N -1 , ω N ), y)7: end for 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 2 s k⊤ • θ k • s k = - 1 2 ℓ s k ⊤ ℓ+1 • θ k ℓ • s k ℓ (60) U k CONV (s k , θ k ) := - 1 2 s k • θ k ⋆ s k = - 1 2 ℓ s k ℓ+1 • θ k ℓ ⋆ s k ℓ (61)Detailed inference algorithm. With the aforementioned details in hand, we re-write the inference algorithm Alg. 1 presented in the main as a Forward subroutine.Algorithm 8 ForwardInput: T , x, W = {θ 1 , ω 1 , • • • ω N } Output: s 1 , • • • , s N -1 or ô depending on the context 1: s 0 ← x 2: for k = 1 • • • N -1 do 3: s k ← Asynchronous T, θ k , ω k , s k-1 ▷ Alg. 6 4: end for 5: ô ← F N s, ω NDetailed implicit EP-BP chaining algorithm. We provide a detailed implementation of our algorithm presented in the main (Alg. 2) in Alg. 11. As usually done for EP experiments, we always perform a "free phase" to initialize the block states (Forward subroutine, Alg. 8). Then, two nudged phases are applied to the last block and parameter gradients subsequently computed, as done classically (BlockGradient subroutine for the last block, Alg. 9), 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. 10), until reaching first block.Details about the implicit differentiation baseline. We describe our implementation of Implicit Differentiation (ID) inside Alg. 12. First, we relax all blocks sequentially to equilibrium following Alg. 8 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. 8), 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.Algorithm 9 BlockGradient (for last block)Input: T , s N -2 ⋆ , θ N -1 , ω N -1 , ω N , β, ℓ, y Output:δs N -2 1: s N -1 β ← Asynchronous T, θ N -1 , ω N -1 , ω N , β, ℓ, y ▷ Alg. 7 2: s N -1 -β ← Asynchronous T, θ N -1 , ω N -1 , ω N , -β, ℓ, y 3: g ω N ← 1 2 ∇ s N -1 ℓ(F N s N -1 β , ω N ) + ∇ s N -1 ℓ(F N s N -1 -β , ω N ) 4: g θ N -1 ← 1 2β ∇ 2 E N -1 (s N -1 β , θ N -1 , s N -2 ⋆ , ω N -1 ) -∇ 2 E N -1 (s N -1 -β , θ N -1 , s N -2 ⋆ , ω N -2 ) 5: g ω N -1 ← 1 2β ∇ 4 E N -1 (s N -1 β , θ N -1 , s N -2 ⋆ , ω N -1 ) -∇ 4 E N -1 (s N -1 -β , θ N -1 , s N -2 ⋆ , ω N -2 ) 6: δs N -2 ← 1 2β ∇ 3 E N -1 (s N -1 β , θ N -1 , s N -2 ⋆ , ω N -1 ) -∇ 3 E N -1 (s N -1 -β , θ N -1 , s N -2 ⋆ , ω N -2 )Algorithm 10 BlockGradient (for all blocks until penultimate) Input: T , s k-1 ⋆ , θ k , ω k , β, δs Output:δs k-1 1: s k β ← Asynchronous T, θ k , ω k , β, δs ▷ Alg. 6 2: s k -β ← Asynchronous T, θ k , ω k , -β, δs 3: g θ k ← 1 2β ∇ 2 E k (s k β , θ k , s k-1 ⋆ , ω k ) -∇ 2 E k (s k -β , θ k , s k-1 ⋆ , ω k ) 4: g ω k ← 1 2β ∇ 4 E k (s k β , θ k , s k-1 ⋆ , ω k ) -∇ 4 E k (s k -β , θ k , s k-1 ⋆ , ω k ) 5: δs k-1 ← 1 2β ∇ 3 E k (s k β , θ k , s k-1 ⋆ , ω k ) -∇ 3 E k (s k -β , θ k , s k-1 ⋆ , ω k ) Algorithm 11 Detailed implicit BP-EP gradient chaining 1: s 1 ⋆ , • • • , s N -1 ⋆ ← Forward (T free , x, W ) ▷ Alg. 8 2: δs ← BlockGradient T nudge , s N -2 ⋆ , θ N -1 , ω N -1 , ω N , β, ℓ, y ▷ Alg. 9 3: for k = N -2 • • • 1 do 4: δs ← BlockGradient T nudge , s k-1 ⋆ , θ k , ω k , βs 1 ⋆ , • • • , s N -1 ⋆ ← Forward (T free , x, W ) ▷ Alg. 8 3: Initialize block states at s 1 ⋆ , • • • , s N -1 ⋆ 4: ô⋆ ← Forward (T nudge , x, W )▷ This time gradients are tracked 5: C ← ℓ(ô ⋆ , y) 6: Backpropagate C backward through the last T nudge steps for each block ▷ e.g. C.backward()We herein present a more fine-grained static gradient analysis to complement the main which was instrumental for our coding work and which builds upon Ernoult et al. [2019] which we briefly summarize below for self-containedness.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:EQUATIONwhere we have ignored the activation function, as done by Ernoult et al. [2019] . Therefore, the computational graph for Eq. ( 62) 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. (62) as:EQUATIONAt the end of the recurrent dynamics Eq. ( 63), a loss is computed based on the output of the last block, namely:EQUATIONwhere s T is the state of the last block after T iterations of the recurrent dynamics Eq. ( 63). In principle, the gradient of the loss C with respect to θ reads as:EQUATIONIn other words, the gradient of the loss C 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:EQUATIONsuch 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:EQUATIONwe define the truncated EP gradient estimate at timestep t through the second phase as:EQUATIONFurthermore, 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 [Ernoult et al., 2019] :EQUATIONThis 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. ( 67) may simply be changed into:EQUATION)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 A.5), the GDU property should still hold.Figure 4 : EP and ID partially computed gradients (( g EP w (t)) t≥0 in black dotted curves and ( g ID w (t)) t≥0 in plain colored curves) going backward through equilibrium for ID and forward through the nudging phase for EP [Ernoult et al., 2019] for a random sample x and associated label y. The ff-EBM comprises 6 blocks and 15 layers in total, with block sizes of either 2 or 3 layers. Each subpanel represents a layer (labelled on the y-axis) with each curve corresponding to a randomly selected weight. "Backward" time is indexed from t = 0 to T = 120, starting from block 6 backward to block 1, with 20 fixed-point iteration dynamics (Eq. ( 56)) being used for both EP and ID within each EB block.Static comparison of EP and ID on ff-EBMs In the light of the above and in order to study the transient dynamics of ID and EP, we define, with W k := {θ k , ω k }:EQUATIONwhere s ±β (t) is computed from Eq. ( 56) with the nudging error current δs k computed with Alg. 2, 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 [Bai et al., 2019] , is highly sensitive to initialization statistics [Agarwala and Schoenholz, 2022] , 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.W k (T -k) C(x, W, y) is the "sensitivity" of the loss C 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)According to GOE, weights are intialized as:W ij ∼ N (0, V N ), if i ̸ = j N (0, 2V N ), if i = jwhere N (µ, σ 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. Scellier et al. [2024] , and sometimes omitted Laborieux et al. [2021] . For these experiments we use both the standard hard activation employed in Ernoult et al.[2019] and [Scellier et al., 2024] , and the more conservative one given in [Laborieux et al., 2021] .Details are given in 3.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 [Laborieux et al., 2021 ] Ernoult et al. [2019] , however, we apply clamping using the ernoult Ernoult et al., 2019 ] laborieux σ(x) = max(min(0.5 × x, 1), 0) [Laborieux et al., 2021] 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 [Laborieux et al., 2021] , 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 out s t ).The cross-entropy cost function associated with the softmax readout is then:l(s, y, W out ) = - C c=1 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 Figure 3 and table 1 are given in table 4. 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 6 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. [1-1-1-1-0-0] [0-0-1-0-1-0-0-0-1-0-0-0-] [0-0-1-0-0-0-1-0-0-1-0-0-0-0]Hyperparameters Hyperparameters and implementation details for 1 and 6 are given in table 5 . 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 1. 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 1 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. Root-finding Algorithm While in principle any root-finding algorithm may be used for the two relaxation phases of our EP implementation (free for inference, and nudged for gradient computation), our implementation utilized a simple fixed-point iteration procedure, in which neuron states are initialized as zero vectors with values updated asynchronously on each iteration to that of the gradient of the total system energy with respect to current state. An approximate illustration of this procedure is found in algorithm 6 As indicated in 4.2, two variants of the convergence procedure were employed, one in which the average value of current state is compared to that of the previous state for each layer, and relaxation is truncated when values for all layers have a difference of less than 1e12.This was known as the tolerance-based (TOL) procedure. Notably, tolerance based convergence criteria were applied on the free phase only, with nudging computed with a fixed value of iterations.While performance overall was seen to increase with smaller TOL threshold, it was observed that 1e12 constitutes an optimal trade-off between performance and runtime without additional tools for accelerating convergence such as Anderson Acceleration.Algorithm 13 Asynchronous with Tolerance (for all blocks until penultimate)Input: T , θ k , ω k , s k-1 ⋆ , β, δs k Output: s k β 1: s k ← 0 2: c ← ∞ 3: for t = 1 • • • T do 4: ∀ odd ℓ ∈ {1, • • • , L k } : s k ℓ,β,temp ← σ ∇ s k ℓ Φ s k β , θ k , s k-1 ⋆ , ω k -βδs k 5: ∀ even ℓ ∈ {1, • • • , L k } : s k ℓ,β,temp ← σ ∇ s k ℓ Φ s k β , θ k , s k-1 ⋆ , ω k -βδs k 6:if t ≥ 2 then s k ℓ,β ← s k ℓ,β,temp 13: end forIn addition to the TOL-based procedure, we obtained results for 4.2 using the more conventional approach of [Scellier and Bengio, 2019 ][Laborieux et al., 2021 ][Ernoult et al., 2019] , applying fixed number of iterations on the first and second relaxation phases (see 6). This approach was also the default used for our scaling experiments in 4.3. Importantly, with the TOL procedure described above in Alg 6 becomes Alg 13Other details All experiments were run using Adam optimizer [Kingma and Ba, 2014] and CosineAnnealing scheduler [Loshchilov and Hutter, 2017 ] with minimum learning rates indicated in 5 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 40BG GPUs.Here we present the results of the splitting experiment as obtained with a fixed iteration relaxation scheme unlike the one presented in section 4.2.Here we find closer agreement between final results obtained from the various block splits, without any significant reduction in WCT.  
discussion: 
 Related work. There is a growing body of work demonstrating scalability of EP on vision tasks.Most notably, Laborieux and Zenke [2022] introduced a holomorphic version of EP where loss gradients are computed with adiabatic oscillations of the model through nudging in the complex plane, and extended to more general implicit models [Laborieux and Zenke, 2023] . Scellier et al.[2022] proposed a fully black-box version of EP where details about the system may not be known.These advances could be readily applied inside our EP-BP chaining algorithm to EB blocks. The work closest to ours, albeit without clear algorithmic prescriptions, is that of Zach [2021] where FF model learning is cast as a deeply nested optimization whose consecutive layers are tied by elemental pair-wise energy functions, as does [Høier et al., 2023] . Such settings 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. ( 4).Limitations and future work. ff-EBM learning partially inherits the pitfalls of BP. Still, nothing prevents FF modules inside ff-EBMs from being trained by any BP alternative. For instance BP can be parameterized by feedback weights to obviate weight transport from the inference circuit to the gradient computation circuit [Akrout et al., 2019] ; BP gradients can be approximated as finite differences of feedback operators [Ernoult et al., 2022] ; or computed via implicit forward-mode differentiation by applying random weight perturbations in the inference circuit [Hiratani et al., 2022 , Fournier et al., 2023 , Malladi et al., 2023] ; local layer-wise loss functions can be used to prevent "backward locking" [Belilovsky et al., 2019 , Ren et al., 2022 , Hinton, 2022] .One natural extension of this study is to incorporate more realism into ff-EBMs. Beyond DHNs, DeepResistive Nets (DRNs) Scellier [2024] Kendall et al. [2020] 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 direction, which nevertheless brings new challenges including non-idealities affecting the inference pathway, such as analog-to-digital and digital-to-analog noise [Rasch et al., 2023 , Lammie et al., 2024] .Finally, considerable work is needed to prove ff-EBM further at scale on more difficult tasks (e.g. standard ImageNet), deeper architectures and novel data.Concluding remarks and broader impact. We show that ff-EBMs constitute a novel framework for deep-learning in heterogeneous hardware settings. We hope that this can help overcome 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 pragmatic, near-term blueprint to mitigate the dramatic carbon footprint of AI training [Strubell et al., 2020] .As fully analog training accelerators remain far from commercial maturity, we believe this work offers an incremental and sustainable plan for gradually integrating analog, energy-based computational primitives as they are integrated into existing digital accelerators. In this section:• We formally define Feedforward-tied Energy-based Models (ff-EBMs) with precise assumptions on the energy-based and feedforward blocks (Def. A.1).• We show that when energy-based blocks comprise a single layer only, the ff-EBM becomes purely feedforward (Lemma A.2).Definition A.1 (ff-EBMs). 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 ⋆ := arg min s E k (s, θ k , x k ⋆ ) ∀k = 1 • • • N -1 ô⋆ := F N (s N -1 ⋆ , ω N ) (10) Finally, we assume that ∀k = 1 • • • N -1, ∇ 2 1 E k (s k ⋆ , θ k , ω k ) is invertible.Lemma A.2. When taking the gradient of E k as defined in Eq. ( 4) 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 ) (11)We consider ff-EBM per Def. (A.1) where the energy functions E k have the form:EQUATION)We assume thatU 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 ) (13) Proof of Lemma A.2. 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 ) (14)Therefore Eq. ( 13) is immediately obtained from Eq. ( 14) with U k = 0.A.2 Proof of Theorem 3.1The proof of Theorem 3.1 is structured as follows:• We restate more formally the multilevel optimization problem introduced in Eq. ( 5).• We directly solve the multilevel problem optimization defined inside Eq. ( 15) using aLagrangian-based approach (Lemma A.3), 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 A.4).• We demonstrate Theorem 3.1 by combining Lemma A.3 and Lemma A.4.• Finally, we highlight that when a ff-EBM is a feedforward net (Lemma A.2), then the proposed algorithm reduces to BP (Corollary A.5.1).Formal statement of multi-level optimization of ff-EBMs.Denoting W := {θ 1 , • • • , θ N -1 , ω 1 , • • • , ω N },the learning problem for a ff-EBM can therefore be written as:min W C(x, W, y) := ℓ(ô ⋆ = F N (s N -1 ⋆ , ω N ), y) (15) s.t. s N -1 ⋆ = arg min s E N -1 (s, θ N -1 , s N -2 ⋆ , ω N -1 ) • • • s.t. s 1 ⋆ = arg min s E 1 (s, θ 1 , x, ω 1 )Here again and similarly to bilevel optimization, solving Eq. ( 5) in practice amounts to computingg θ k := d θ k C and g ω k := d ω k Cand performing gradient descent on θ k and ω k .Lemma A.3 (Lagrangian-based approach). Assuming a ff-EBM (Def. A.1), we denotes 1 ⋆ , x 1 ⋆ , • • • , s N -1 ⋆, ô⋆ the states computed during the forward pass as prescribed by Eq. ( 10). Then, the gradients of the objective function C := ℓ(ô(s N -1 ⋆ ), y) as defined in the multilevel optimization problem (Eq. ( 5)), where it is assumed that ℓ is differentiable, read:EQUATIONwhereλ 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,3 E k+1 s k+1 ⋆ , θ k+1 , s k ⋆ , ω k+1 • λ k+1 ⋆ + ∇ 2 1 E k s k ⋆ , θ k , s k-1 ⋆ , ω k • λ k ⋆ = 0 (17) Proof of Lemma A.3. 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. ( 5) reads:L(s, λ, W ) := ℓ ô(s N -1 ), y + N -1 k=1 λ k ⊤ • ∇ 1 E k (s k , θ k , s k-1 , ω k ), s 0 := x (18)Writing the associated Karush-Kuhn-Tucker (KKT) conditions ∂ 1,2 L(s ⋆ , λ ⋆ , W ) := 0 satisfied by optimal states and Lagrangian multipliers s ⋆ and λ ⋆ , we get :EQUATIONEq. ( 19) governs the bottom-up block-wise relaxation procedure (as depicted in Alg. 1), while Eq. ( 20)and Eq. ( 21) governs error propagation in the last block and previous blocks respectively. Given s ⋆ and λ ⋆ by Eq. ( 19) -Eq. ( 21), the total derivative of the loss function with respect to the model parameters read: 20)-( 21))EQUATIONEQUATION)) +d W λ ⊤ ⋆ • ∂ 2 L(s ⋆ , λ ⋆ , W ) =0 (Eq. (EQUATIONMore precisely, applying Eq. ( 22) to the feedforward and energy-based block parameters yields:d ω N ℓ(ô ⋆ , y) = ∂ 2 F N (s N -1 ⋆ , ω N ) ⊤ • ∇ 1 ℓ(ô ⋆ , y), d θ k ℓ(ô ⋆ , y) = ∇ 2 1,2 E k (s k ⋆ , θ k , s k-1 ⋆ , ω k ) • λ k ⋆ ∀k = 1 • • • N -1 d ω k ℓ(ô ⋆ , y) = ∇ 2 1,4 E k (s k ⋆ , θ k , s k-1 ⋆ , ω k ) • λ k ⋆ ∀k = 1 • • • N -1Lemma A.4 (Computing Lagrangian multipliers by EP). Under the same hypothesis as Lemma A.3, we define the nudged state of block k, denoted ass k β , implicitly through ∇ 1 F k (s k β , θ k , x k ⋆ , δs k , β) = 0with:EQUATION)EQUATION)then we have:EQUATIONwhere (λ k ) k=1•••N -1 are the Lagrangian multipliers associated to the multilevel optimization problem defined in Eq. ( 5).Proof of Lemma A.4. We prove this result by backward induction on k.Initialization (k = N -1). By definition, s N -1 β satisfies :EQUATIONDifferentiating Eq. ( 26) with respect to β and evaluating the resulting expression at β = 0, we obtain:EQUATIONSubstracting out Eq. ( 20) defining the Lagrangian multiplier λ N -1 ⋆ and Eq. ( 27), we obtain:∇ 2 1 E N -1 s N -1 ⋆ , θ N -1 , s N -2 ⋆ , ω N -1 • d β s N -1 β | β=0 -λ N -1 ⋆ = 0 (28) By invertibility of ∇ 2 1 E N -1 s N -1 ⋆ , θ N -1 , s N -2 ⋆, ω N -1 , we therefore have that:λ N -1 ⋆ = d β s N -1 β | β=0 (29) 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:EQUATIONOn the one hand, proceeding as for the initialization step, differentiating Eq. ( 30) with respect to β and taking β = 0 yields:EQUATION)On the other hand, note that δs k rewrites :EQUATIONwhere we used at the last step the recursion hypothesis. Therefore combining Eq. ( 31) and Eq. ( 32), we get:EQUATIONFinally, we substract out Eq. ( 21) and Eq. ( 33) to obtain:EQUATION)We conclude again by invertibility of∇ 2 1 E k (s k ⋆ , θ k , s k-1 ⋆ , ω k ) that λ k ⋆ = d β s k β | β=0 .Theorem A.5 (Formal). Assuming a model defined by Def. (A.1), we denote s1 ⋆ , x 1 ⋆ , • • • , s N -1 ⋆ , ô⋆the states computed during the forward pass as prescribed by Alg. 1. We define the nudged state of block k, denoted as s k β , implicitly throughEQUATION)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:EQUATIONthen the following chain rule applies:δs N -1 := ∇ s N -1 ℓ(ô ⋆ , y), g ω N = ∂ 2 F N s N -1 ⋆ , ω N ⊤ • ∇ 1 ℓ(ô ⋆ , y) (37) ∀k = 1 • • • N -1 :    ∆x k = d β ∇ 3 E k (s k β , θ k , x k ⋆ ) β=0 , g θ k = d β ∇ 2 E 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 (38)Proof of Theorem A.5. Combining Lemma A.3 and Lemma A.4, the following chain rule computes loss gradients correctly:EQUATIONTherefore 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 β ∇ 3 E k (s k β , θ k , x k ⋆ ) β=0 (41) ∇ 2 1,2 E k (s k ⋆ , θ k , s k-1 ⋆ , ω k ) • d β s k β | β=0 = d β ∇ 2 E k (s k β , θ k , x k ⋆ ) β=0 (42) ∇ 2 1,4 E k (s k ⋆ , θ k , s k-1 ⋆ , ω k ) • d β s k β | β=0 = ∂ 2 F k s k-1 ⋆ , ω k ⊤ • d β ∇ 3 E k (s k β , θ k , x k ⋆ ) β=0 (43) Let k ∈ [1, N -1].We prove Eq. ( 41) 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 β ∇ 3 E k (s k β , θ k , x k ⋆ ) β=0Eq. ( 42) can be obtained as:∇ 2 1,2 E 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 β ∇ 2 E k (s k β , θ k , x k ⋆ ) β=0Finally and similarly, we have:∇ 2 1,4 E k (s k ⋆ , θ k , s k-1 ⋆ , ω k ) • d β s k β | β=0 = d β ∇ 4 E 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 ⋆ ) β=0Corollary A.5.1. Under the same hypothesis as Theorem A.5 and Lemma A.2, 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 (44) Proof of Corollary A.5.1. Let k ∈ [1, N -1].As we can directly apply Theorem A.5 here, proving the result simply boils down to showing that:∆x k = σ ′ (x k ⋆ ) ⊙ δs k (45)First, we notice that when E k is of the form of Eq. ( 12), then ∆x k reads as:EQUATIONs k β satisfies, by definition and when U k = 0:EQUATIONCombining Eq. ( 46) and Eq. ( 47) yields Eq. ( 45), and therefore, along with Theorem A.5, the chain-rule Eq. ( 44).We showcase in Alg. 3 the resulting algorithm with finite β and implicit BP-EP chaining, with lines in blue highlighting differences with the general algorithm Alg. 2.Algorithm 3 Implicit BP-EP gradient chaining with U k = 01: δs, g ω N ← ∇ s N -1 ℓ(ô ⋆ , y), ∇ ω N ℓ(ô ⋆ , y) ▷ Single backprop step 2: for k = N -1 • • • 1 do 3: s β , s -β ← σ x k ⋆ -βδs k , σ x k ⋆ + βδs k ▷ EP through E k 4: g ω k ← 1 2β ∇ 4 E k (s β , θ k , s k-1 ⋆ , ω k ) -∇ 4 E k (s -β , θ k , s k-1 ⋆ , ω k ) ▷ i-BP through F k 5: δs ← 1 2β ∇ 3 E k (s β , θ k , s k-1 ⋆ , ω k ) -∇ 3 E k (s -β , θ k , s k-1 ⋆ , ω k ) 6: end for A.3 Explicit BP-EP chainingWe presented in Alg. 2 a "pure" EP algorithm where the BP-EP gradient chaining is implicit. We show below, inside Alg. 4, 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 A.5. We also show the resulting algorithm when the ff-EBM reduces to a feedforward net (Lemma A.2) inside Alg. 5, highlight in blue the statements which differ from the general case presented inside Alg. 4.Algorithm 4 Explicit BP-EP gradient chaining (Theorem (3.1))1: δs, g ω N ← ∇ s N -1 ℓ(ô ⋆ , y), ∇ ω N ℓ(ô ⋆ , y) ▷ Single backprop step 2: for k = N -1 • • • 1 do 3: s β ← Optim s E k (s, θ k , x k ⋆ ) + βs ⊤ • δs ▷ EP through E k 4: s -β ← Optim s E k (s, θ k , x k ⋆ ) -βs ⊤ • δs 5: g θ k ← 1 2β ∇ 2 E k (s β , θ k , x k ⋆ ) -∇ 2 E k (s -β , θ k , x k ⋆ ) 6: ∆x ← 1 2β ∇ 3 E k (s β , θ k , x k ⋆ ) -∇ 3 E k (s -β , θ k , x k ⋆ ) 7: g ω k ← ∂ 2 F k s k-1 ⋆ , ω k ⊤ • ∆x ▷ Explicit BP through F k 8: δs ← ∂ 1 F k s k-1 ⋆ , ω k ⊤ • ∆x 9: end for Algorithm 5 Explicit BP-EP gradient chaining with U k = 0 1: δs, g ω N ← ∇ s N -1 ℓ(ô ⋆ , y), ∇ ω N ℓ(ô ⋆ , y) ▷ Single backprop step 2: for k = N -1 • • • 1 do 3: ∆x ← -1 2β σ x k ⋆ -βδs k -σ x k ⋆ + βδs k 4: g ω k ← ∂ 2 F k s k-1 ⋆ , ω k ⊤ • ∆x ▷ Explicit BP through F k 5: δs ← ∂ 1 F k s k-1 ⋆ , ω k ⊤ • ∆x 6: end for A. 
 
 Appendix:  



Reviews: review: title: title: Interesting Hybrid Approach in Simulation
main_review: main_review: This paper describes an algorithm for training mixed feed forward (FF) and energy based (EB) networks. The methods are clear, and the results are encouraging. The observation that implementations of EB algorithms are far from mature, and therefore there is room for hybrid approaches, is a good one.

A few minor comments and questions:

Do the authors have any way of estimating the type of energy savings such algorithms could provide? That seems to be a main motivation of EB algorithms, and presumably the savings cannot be higher than the fraction of the network that is EB instead of FF. Or is there another motivation for including EB that is less straightforward?  

The reference to Eq (2) should be to the second (unnumbered) equation, I believe. 

The scale of Fig 2 makes it impossible to tell if there is any variation between layers – could the authors rescale or plot on a log scale 1- similarity?

The authors hint at but do not explicitly confirm that the advantage their system has over purely EB methods are nonlocal operations like batchnorm – is this the source of improved performace?
title: title: The paper introduces a way to combine feed-forward and energy-based blocks in models for neuromorphic computing
main_review: main_review: This paper derives expressions and algorithms for performing inference and weight update on models of deep networks that contain alternating blocks of feed forward layers and energy-based models. This was previously hard to do because, while feed forward layers are simple functions, energy-based models require internal optimization and an efficient algorithm to train (like Equilibrium Propagation). 

This is a nice paper, though I am not completely convinced by its motivation. Since this combining of feed-forward and energy-based components, both backprop and equilibrium propagation are required, it is unclear what benefit is computationally obtained by this combination. EP and its variants make sense for physical analog implementations precisely because only such physical local rules work for them. I can imagine combining physical self-learning elements using EP with computational models that use BP, but the advantage for doing that is not apriori clear. The authors partially address this by noting that you can break down the model into ‘independent’ EB blocks that can significantly reduce computations time, as shown in Table 3. This should be communicated in the introduction, and possibly also in the abstract to help motivate the work.

I find this paper interesting and definitely meriting acceptance to the MLNCP workshop. Some minor comments are appended below:
1) Why is EP incompatible with common non-stationary operations such as activation functions. Is this because such functions cannot be simply written as “energy-based” elements? There are most likely ways of fixing this in the majority of cases.
2) Eq. 2 is noted incorrectly in the text. It probably should refer to the minimization of C.
3) What does “L” in EBL stand for?
3) Indexing of F,w on ln 93 may be incorrect. Didn’t you mean F^2,w^2?
4) In Figure 2, given the y-scale from 0-1 there seems to be no difference and all layers do equally well. It might be more informative to show log⁡(1-cos⁡(g^EP - g^AD ) ).