Towards training digitally-tied analog blocks via hybrid gradient computation
Authors: Timothy Nest
Maxence Ernoult
Venue: NeurIPS 2024 spotlight (2024)
Keywords: implicit differentiation
equilibrium propagation
bilevel optimization
hopfield networks
analog computing
hardware-aware training
backprop
energy-based models
physical learning
Primary area: machine_learning_for_other_sciences_and_fields
TL;DR: We extend Equilibrium Propagation (EP) to a novel hardware-realistic model which comprises feedforward and energy-based blocks, resulting in chaining backprop and EP gradients backward through these blocks and a new SOTA performance on ImageNet32

Abstract: Power efficiency is plateauing in the standard digital electronics realm such that new 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 a 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 platforms. Such heterogeneous hardware lacks a supporting theoretical framework. In this work, we introduce \emph{Feedforward-tied Energy-based Models} (ff-EBMs), a hybrid model comprised of feedforward and energy-based blocks housed on 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 flexibly on realistic architectures. We experimentally demonstrate the effectiveness of this approach on ff-EBMs using Deep Hopfield Networks (DHNs)  as energy-based blocks, and show that a standard DHN can be arbitrarily split into any uniform size while maintaining or improving performance with increases in simulation speed of up to four times. We then train ff-EBMs on ImageNet32 where we establish a new state-of-the-art performance for the EP literature (46 top-1 \%). Our approach offers a principled, scalable, and incremental roadmap for the gradual integration of self-trainable analog computational primitives into existing digital accelerators.

Towards training digitally-tied analog blocks via hybrid gradient computation
Timothy Nest"* timothy. nest @milal quebec
Maxence Ernoult maxencedrain.c ai
Abstract
Power efficiency is plateauing in the standard digital electronics realm such that new hardware, models, and algorithms are needed to reduce the costs of Al training. The combination of energy-based analog circuits and the Equilibrium Propagation (EP) algorithm constitutes a 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 son--eeght-statoonary operations, mitigate analog device imperfections, and leverage existing digital platforms. Such heterogeneous hardware lacks a supporting theoretical framework. this work, we introduce Feedforward hybrid model comprised of feedforward and energy-based blocks housed on digital and analog circuits. We derive novel algorithm compute gradients end-to-end in ff-EBMs by backpropagating and "ot-propagating" through feedforward and energy- parts respectively, enabling EP to be applied flexibly on realistic architectures. We experimentally demonstrate the effectiveness of this approach on ff-EBMs using Deep Hopfield Networks (DHNs) as energy-based blocks, and show that standard DHN can be arbitrarily split into any uniform size while maintaining or improving performance with increases in simulation speed of up to four times. We then train ff- FBMM ImageNet32 where we establish a new state-of-the-an performance for the EP literature (46 top-1 %) F Our approach offers a principled, scalable, and incremental roadmap for the gradual integration self-trainable analog computational primitives into existing digital accelerators.
Energy-based Models (ff-EBMs), a
1 Introduction
Gradient-basee optimization, the cornerstone and most energy greedy component of deep learning, fundamentally relies upon three factors: i) highly parallel digital hardware such GPUs, feedfor- ward models and iii) backprop (BP). With skyrocketing demands of Al compute, reducing the energy consumption of AI systems has become matter of great economic, societal and environmental urgency [Strubell et al., 22000, calling for the exploration of novel compute paradigms (Thompson
al. 2020 Scellier 2021 Stern and Murugan 2023
One promising path towards this goal analog in-memory computing [Sebastian et al. 2020 by mapping weights onto crossbar of resistive devices, Kirchoff current and voltage laws in- herently perform matrix-vector multiplications with constant time complexity [Cosemans et al. 2019]. By stacking multiple such crossbars, entire neural network be mapped onto physical system. An important formalism for such a system is that of energy-based (EB) analog circuits
"Montreal Institute of Learning Algorithms (MILA)
Rain AI
Equul contribution
"Ou code is available on ottppp//ghaaaa.cc..aaaaaaapppppp iitp:///rara..ppp....ffiiiilllll official
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

[Kendall et al., 2020. Stern et al., 2023 Dillavou et al. 2023 Scellier 2024], which are "self- learning" systems that can compute loss gradients through two relaxations to equilibrium (i.e. two "forward passes"). Such a procedure falls under the umbrella of energy-based learning (EBI.) al- gorithms [Scellier et al., [024]. One such algorithm, Equilibrium Propagation (EP) [Scellier and Bengio 2017]. particularly stands out for its strong theoretical guarantees, relative scalability in the realm of backprop alternatives [Laborieux and Zenke 2022 2023] and proven application and substantial speedups com- pared to its GPU-based counterpart (Yi et al., 2023]. This suggests a new alternative compute paradigm for gradient-basse optimization consisting of: i) analog hardware, ii) EBMs, and iii) EP.
on small analog systems with 10, x00x greater energy-efficiene
In this paper, we propose a theoretical framework for ex- tending end-to-end gradient computation to a realistic set- ting where the system in question may or may not be fully analog. Such a setting is plausible in the near term, for two major reasons. First, analog circuits exhibit non-ideal physical behaviors which affect both the inference pathway [Wang et al.,
Feerforward
P. Let)
VV L8))
Energy based
2023, Ambrogio et al. 2023] and parameter optimization Nandakumar et al. 2020. Spoon et al. 2021
Backward pais
2024]. compromising performance. Second, owing to the latency and energy cefergy-oonuumptiin of resistive devices' write operations, analog circuits should be fully weight stationary weights must be written before the inference procedure begins which excludes many oper- ations used conventionally in machine learning including activation functions, normalization, and attention [Spoon et al. 2021 Jain et al. 2022 Liu et al., 2023 Li et al.
Lammie et al.
4,
Figure 1: Illustrating BP-EP backward gradient chaining through feedforward (red) and energy-based (blue) blocks, ac- counting for digital and analog circuits respectively.
2023 Therefore, analog systems are likely to be used in combination with auxiliary digital circuitry, resulting in hybrid mixed precision systems (Haensch et al. 2018 While the design of purely inferential engines made up of analog and digital parts is nearing commercial maturity AAmbbogio et al. 20233 in-situ learning of such systems has barely been explored. An important challenge remains in proving EBL algorithms can scale in a manner comparable to backprop, given the requirement of simulating EB systems on GPUs. Because of 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 (5-6 layer) models (Scellier et al.) 2024 Scellier 2024].
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 modeling digital and analog parts as feedforward and EB modules respectively, we show how backprop and EP error signals can be chained end-to-end via feedforward and EB blocks respectively in a principled fashion. Rather than opposing digital and analog, or backprop and "alternative" learning algorithms, as is 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 3.1 as high-level models of mixed precision systems whose inference pathway read as the composition of feedforward and EB modules (Eq. (5) Alg. D
We show that gradients in ff-EBMs be computed in an end-to-end fashion (Section 33 backpropagating through feedforward blocks and "eq- "ro--rropagatin" through EB blocks (Theorem 3.1
Alg. 2) and that this procedure is rooted in deeply-nested optimization problem (Section) B.2).
Finally, we experimentally demonstrate the effectiveness of our algorithm on ff-EBMs where EBM blocks are Deep Hopfield Networks (DHNs) (Section 4).In particular we show that i) final
and transient gradient estimates computed by our algorithm (Alg. 2 near perfectly match gradients computed by end-to-end automatic differentiation (Section 4.2) which also prove mathematically
(Theorem 4. ii) a standard DHN model can be arbitrarily split into a ff-DHN with the equivalent
2

---

layers and architectural layers while maintaining improving performance, remaining on par with automatic differentiation and being up to four times faster to simulate depending on the convergence criterion at use to compute equilibrium (Section 4.3). iii) the proposed approach yields 46 % top-1 (70% top-5) validation accuracy on ImageNet32 when training a ff-EBM of 15 layers, beating current state-of-the-art for EP by large margin, without relying on holomorphic transformations inside EBM blocks [Laborieux and Zenke 2022 2023
Background
Notations. Given A 12 + RM differentiable mapping, we denote its total derivative with respect to Sj as de == d.A(s)/ds) E IR",, its partial derivative with respect to sy as As)) == dA(ss)0ss E R.. When A takes scalar values (m 1), its gradient with respect to 8y is denoted as v,A(s) : A
2.1
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 energy function. As such they particular kind of implicit model. Namely, an EBM is defined by (scalar) energy function : 8, 0, x + E(s, 0, E IR where x, 8, and 0 respectively denote static input, hidden and output neurons and model parameters, and each such tuple defines a configuration with an associated scalar energy value. Among all configurations for a given input and some model parameters 0, the model prediction s, implicitly given as an equilibrium state which minimizes the energy function:
arg min (s,
2.2 Standard bilevel optimization
Assuming that V Exr, S., 0) is invertible, note that the equilibrium state 8, implicitly depends on x and 0 by virtue of the implicit function theorem (Dontchev al. 2009]. Therefore our goal when training an EBM-in a supervised setting, for instance - is to adjust the model parameters such that s. (z, 0) minimizes some cost function é s + {(s, y) E where y is some ground-trute label associated to x. More formally, this learning objective can be stated with the following bilevel optimization problem [uccheet and Sacramento 2022]:
min C(x, 0, y) == ((s.,)) v) s.t. 8. arg min E(s, 0, 2)
(2)
0
Solving Eq. (2) in practice amounts to computing the gradient of its outer objective C(r,, 0) with respect (doC((1 0)) and then performing gradient descent over 0.
2.3 Equilibrium Propagation (EP)
An algorithm used to train an EBM model in the sense of Eq. may be called an EBL algorithm [Scellier et al. 2024). Equilibrium Propagation (EP) (Scellier and Bengio 2017) is an EBL algorithm which computes an estimate of doC 8) with at least two phases. During the first phase, the model is allowed to evolve freely s. arg min, E(s, z). Then, the model is slightly nudged towards decreasing values of cost l and settles to second equilibrium state s3. This amounts to augmenting the energy function E by an additional term BE(x, y) where 3 IR" is called the nudging factor. Next, the weights are updated to increase the energy of 8, and decrease that of 83, thereby "contrasting" these two states. More formally, (Scellier and Bengio [017] show in the seminal EP paper:
8/g = arg min (x .r) + (8(x,
APEP
(2 0,0) - V F2 0, R)),
(3)
where a denotes some learning rate. EP comes in different flavors depending on the sign of / inside
Eq. (3) or on whether two nudged states of opposite nudging strengths (+,3) are contrasted, a variant called Centered EP (C-EP) which was shown to work best in practice [Laborieux et al. 2021 Scellier
ct al. 2024 and reads as:
A C-P
a
-
(4)
23
3

---

3 Tying energy-based models with feedforward blocks
In the present section introduce new model, Feedforward-tied EBMs (ff-EBMs, section 3.1]. which read as composition of feedforward and EB transformations (Alg. T We show how optimizing
ff-EBMs amounts solving multi-level optimization problem (Section 3.2) and propose a BP-EP
gradient chaining algorithm as a solution(Section [.3 Theorem 3.1 Alg. 2). We highlight as an
edge case the fact that ff-EBMs reduce to standard feedforward nets (Lemma A. 1 and the proposed
BP-EP gradient chaining algorithm to standard BP (Corollary) [A.I] when each EB block comprises a single hidden layer. We highlight in red and blue the parts the model and associated algorithms performed inside feedforward (digital) and EB (analog) blocks respectively.
3.1 Feedforward-tied Energy- based Models (ff-EBMs)
Inference procedure. We define Feedforward-tied Energy-based Models (ff-EBMs) compo- sitions of feedforward and EB transformations. Namely, an data sample x is fed into the first feedforward transformation F* parametrized by some weights w which yields an output æ!. Then, x. is fed as a static input into the first EB block E with parameters H¹, which relaxes to an equilib- rium state sy. s is in turn fed into the next feedforward transformation F' with weights wt and the above procedure repeats until reaching the output layer ô. More formally, denoting F and E* the kth feedforward and EB blocks parametrized by the weights w and 6² respectively, the inference pathway of ff-EBM reads as:
x = F w :- arg a- min E (s, x!) = N -
(5)
02 = p. (x N )
The ff-EBM inference procedure depicted more compactly inside Fig. (left) and Alg.
V,Fr( m0
F, On St =0
Algorithm 1 ff-EBM inference (Eq. (5)
X+
1: 8 2 x
2: for k 1...N do
P/a mn)
aF (e as,"
3: I t- pk ws,
4: - Optim (s,
5: end for
Figure 2: Depiction of the forward (left) and backward (right) pathways through a ff-EBM, with blue and pink blocks denoting EB and feedforward transformations.
6: Ô t FN (s,
Form of the energy functions. We further specify the form of the energy of the kth EB block of a ff-EBM as defined per Eq. (5). The associated energy function E takes some static input x² from the output of the preceding reedforward transformation, has hidden neurons & and is parametrized by weights 0b More precisely:
P := (s**) - th + Uh (s*,
Eq. (6) reveals three different contributions to the energy. The first term determines the non- ooneaineay applied inside the EB block (Zhang and Brand 2017 Høier et al. 2023) for given invertible and continuous activation function a, G is defined such that VG o (see Appendix A.1.3)
The second term inside Eq. (6) accounts for a purely feedforward contribution from the previous feedforward block f.r. Finally, the third term accounts for internal interactions within the layers of the EB block.
Recovering feedforward net. When taking the gradient of Ex as defined in Eq. (6) with respect to s and zeroing it out, it can be seen that s is implicitly defined as:
(a*
(sk (st,0*))

---

A noteworthy edge case highlighted by Eq. 7 is when U = 0 for all k's, i.e. when there are no intra-block layer interactions, equivalently when the EB block comprises single layer only. In this case, Sf. is simply feedforward mapping x* through and in turn the ff-EBM is simply a standard feedforward architecture (see Lemma A. inside Appendix A.I.I).
3.2 Multi-level optimization of ff-EBMs
Just as learning EBMs can be naturally cast as bilevel optimization problem, learning ff-EBMs equates multi-level optimization problem where the variables being optimized in the inner subproblems are comprised of EB block variables s .... . .N-1 To make this clearer, we re-write the energy function of the tth block E from Eq. (6) highlight the dependence between two consecutive EB block states:
(s 0² s ,+=+,,=))=== E (s* 0h, pk (6f-1w(-1)
(8)
It can be seen from Eq. (8) that the equilibrium state of obtained by minimizing E* will be dependent upon the equilibrium state sf-1 of the previous EB block, which propagates back through prior EB blocks. Denoting W = {0) .... 10-11..' ,w!,... ".....""" w }. the learning problem for ff-EBM can therefore be written as:
min C(s,, W, = E(o. F (s ).
(9)
s.t. 2 N-1 arg min E. n (s. NN sh W-
s.t. $ arg min E'(s, 0 w²)
Here again and similarly to bilevel optimization, solving Eq. (9) in practice amounts to computing Ook == dgs and Jar == d.xC and performing gradient descent on 6 and wt.
3.3 A BP-EP gradient chaining algorithm
Main result: explicit BP-EP chaining. Based on the multilevel optimization formulation of ff- EBMs learning Eq. (9. we state the main theoretical result of this paper Theorem 3.1 (see proof in Appendix II
Theorem 3.1 (Informal). Assuming a model of the form Eq.
NN ô, the states computed during the forward pass as depicted in Alg. We define the nudged state of block k. denoted as 8 implicitly through NrJ* (sty 8h .. 5s (3) with:
55) we denote s x.,
F* 0² x ,,) ;= (s 0h, (xx
Bs
ōs
(10)
Denoting k and A the error signals computed at the input of the feedforward block E and of the EB block E² respectively, then the following chain rule applies:
N = VoN-il(ōn.. u), IWN VuN ((0x.y)
(11)
Vk = 2 ....N -
k d3 (VaL' (s 0h, .att))
0=0" Fot = dg (V2 VV (st 0h (x
8-0
(12)
8sk-1 Fb (st-1, ww)) Aak, Dut = Fb (s ww)) T Art
Proposed algorithm:
implicit BP-EP chaining. Theorem 3.1 reads intuitively. It prescribes an explicit chaining of EP error signals passing backward through E (5s - A and BP error signals
passing backward through a (Ax -> 88k-1) which directly mirrors the ff-EBM inference pathway depicted in Fig.
Yet noticing that:
{ 85k Ft (sh- wot=),
T Art = dg (v3E* (st, sk-1 wat)) .aat))
B=0
Iut F, (s ) Art dg (V.Et (st HK sh-1 .aat)) 3=0"
the same error signal can by passed through Ee (8s" + ksk-1) where BP and EP are implicitly chained inside E (see Appendix ). This insight, along with centered scheme to estimate
5

---

derivatives with respect to 8 around O as done for the C-EP algorithm (Eq. 4 motivates the implicit BP-EP gradient chaining algorithm in Alg. 2 we used for our experiments (see Alg. insiiee Appendix A.3.1 for its explicit counterpart). Given that the proposed algorithm appears as a generalization EP, we refer to Alg. Las "EP" in the experimental section, for simplicity.
Algorithm Implicit BP-EP gradient chaining (Theorem (3.1)
1: ös, INN + V
N ( y), V [..((ó
D Single BP step
2: for k N
do
3: $3 Optim (s 8h sh-1 ,w + Bs ôs
EP in E
4: 8-8 Optim
E (s, 0h sh w Bs&T -8s
5: D00 23
V2 V2 (88, 0+, 8 - w V2 F²F² (8 B 0h sh- w V ($3. - wh) V,E" (s 8- Qh s ,/-1,,0)) ww))
6: Iuk +
Implicit BP in F
2
7: òs
V F (83, s4 wh F3 ( (s-B. a. Q #-1 (.-1w(
23
8: end for
Recovering backprop. When the ff-EBM under consideration is purely feedforward (U* 0), we
show that Eqs. I1] [2) reduce to standard BP through feedforward net (Corollary |A. Alg. 6 and Alg. | Appendix I)2.I). Since this case is extremely close to standard BP through feedforward nets, do not consider this setting in our experiments.
4
Experiments
In this section, we present the ff-EBMs used in our experiments (Section 4.1 and carry out static gra- dient analysis computing and analyzing ff ff-EBM parameter gradients for some and y (Section (42). We extend the observation made by Ernoult et al. [2019] -that transient EP parameter gradients obtained during the second phase match those computed by automatic differentiation through equilib- rium and across blocks- to ff-EBMs (Fig. (3) of automatic differentiation and EP in our framework, are near perfectly aligned (Fig- 55) We then show on the CIFAR-10 10 task that performance of ff-EBMs can be maintained improved across various block splits maintaining the same number of layers, while remaining on par with automatic differentiation
Theorem[4.] showing that gradient estimates
4.3). We show furthermore that blocks of smaller size are up to four times faster to simulate depending on the convergence criterion at use for computing equilibrium inside EB CIFAR- -100 and ImageNet32
blocks. Finally, perform further ff- ff-EBM training experiments where we establish a new state-of satee--fftheaaaee f-tthe-art art performance in the EP literature (Section
4.1 Setup
Models.
Using the same notations as in Eq. (6). the ff-EBMs at use in this section are defined:
{ Ukolot, OK) =- ATA 0# (s gh) = 1 (eh * )
Flx (s A- ( = BN (P (uČONV =²-1) swi w Pkk (s-1)) =
(13)
with (N: wi wt)). P and the batchnorm, pooling and convolution operations, the generalized
dot product for tensors and s² (sfT. st the state of block l: comprising L layers. Such EBM blocks are known as Deep Hopfield Networks (DHNs). DHNs are comprised of fully connected (Ufc) and convolutional operations (U&oNy) forming a symmetric weight matrix 6 with sparse, block-wise structure such that each layer s is bidirectionally connected its neighboring layers
*t= and */41 through connections oi_, and 0%* respectively (see Appendix A.1.3). To empirically ensure convergence, the non-linearity a applied within EB blocks is (x) == min (max (az, 0), with a (0, 1). Finally, two design choices instrumental to the success of ff-EBM gradient computation and subsequent training. First, we initialized the weights of U and U using

---

Gaussian Orthogonal Ensembles (GOE) AAgawwaaa and Schoenholz, 2022 to enable faster equilib- rium computation (see next paragraph). Second, while the last layer of given block was simply passed as an input to the next block (i.e. using Fib in Eq. (13) for small enough models (7. = 6 inside the experiment depicted in Section (43), the use of batchnorm layers in between blocks (i.e.
using FBN in Eq. 13 becomes essential for deeper models.
Equilibrium computation. As depicted in Alg. 2 the steady states 8 may be computed with any loss minimization algorithm. Here, in past works on EP (Eroolt et al., 2019 Laborieux et al., 2021 Laborieux and Zenke, 2022] [Scellie et al. 2024], we employ a fixed-point iteration scheme to compute the EB blocks steady states. Namely, iterate Eq. 7 until reaching equilibrium (the same
scheme is used for ff-EBM inference, Alg. [ with B = O.):
LB, +- a (x* VIU 2, 0*)
F 36s*))
(14)
Important details about how Eq. (4) is executed in practice have to be highlighted. First, we employ scheme to asynchronously update even (s)g) and odd (shg layers (Scellier al. 2024 see Appendix A.1.3] Second, Eq. (I4) were either executed for a fixed and predetermmmed number of steps as done in the aforemmentoned EP literature, or using an C- polerance-basee convergence criterion (TOL) which stops executing Eq. [4] when (s / s S on average see Appendix A.5.3 for details.
Algorithm baseline. As algorithmic baseline, we simply use automatic differentiation (AD) backward through the fixed-point iteration scheme Eq. (14) with 3 = and directly initializing s/=0 8. (Fig. 4). This version of AD, where we backpropagate through equilibrium, is known as Pineda 1987] or Implicit Differentiatinn (ID).
"Recurrent Backpropagation" [Almeida 1987
4.2 Static comparison of EP and ID on ff-EBMs
In order to study the transient dynamics of ID and EP, define, with W = {0t, w
ID T
(t) = - dw Cw)(r-s)((rr, W*, v),
k
EP Tip (B.t) = 2,9 V Ek (5 Wh s--1) - V V E Ba Wh s 1)-
(15)
where ah. is computed from Eq. 140 with the nudging error current 8s computed with Alg. T is the total number of iterations used for both ID and EP in the gradient computation phase.
and
Block
Block
Block
Block
Block
Block
30
90
100
100
110
120
10 20
40
50
60
60
70
80
Backward time
Figure 3: EP and ID partially computed gradients (( (t)) in black dotted curves and (G (20 in plain colored curves) going backward through equilibrium for ID and forward through the nudging phase for EP [Emoult et al., 2019] for a random sample x and associated label y. The ff- ff-EBM comprises 6 blocks and 15 layers in total, with block sizes of either 2 or 3 layers. Each sub-panel represents layer (labeled on the y-axis) with each curve corresponding to a randomly selected weight. "Backward" time is indexed from = to T 120, starting from block 6 backward to block 1. with 20 txxettpoint iteration dynamics (Eq. [4] being used for both EP and ID within each block.
For a given block k; dw W. y) is the "sensitivity" of the loss C to parameter W* at timestep T 1 so that T. (t) is ID gradient truncated at T f. Fig. depicts the computational graph that is differentiated through when using ID and shows where 9 (t) are obtained correspondingly. Similarly, 3 (1) is an EP gradient truncated at t steps forward through the nudged phase. When T is sufficiently large, IID (T)) and @PP. (T)) converge to Cz, W, y). Fig. 3 displays (o (t))
and (GEP, (t)) an heterogeneous ff-EBM of 6 blocks and 15 layers (16 if counting the last linear
7

---

"readout" layer computing the logits) with blocks comprising 2 or 3 layers for a randomly selected sample and its associated label y - see caption for detailed description. It can be seen EP and ID error weight gradients qualitatively match very well throughout time, across layers and blocks. We also display the cosine similarity between the final EP and ID weight gradient estimate Iw* (T) EPP. (T) for each layer and observe that EP and ID weight gradients are near perfectly aligned.
and gw
Theorem 4.1 generalizes the equivalence between EP and ID to ff-EBMs [Ernoult al. 2019
Theorem 4.1 (Informal). Assuming vk = 1 N s
8
Vk- = 1. N - 4t 00.... (t = Ww = lim EP (3,t)
(16)
B-0
W
W
W>
Block
Block (k
Figure 4: Light grey: computational graph associated with f--BM inference (Alg. D when applying fixed-point iteration to compute equilibrium states within each block (Eq. (4)) where the node st denotes the state of block & (comprising several layers) at timestep f. Blue arrows: backward auto- matic differentiation (AD) through the computational graph where @W (t) is the partially computed gradient truncated at T t. Since the states which are differentiated through are taken at equilibrium (sf = $ vtt = 0 . r) this instantiation of AD can be viewed as Implicit Differentiation (ID).
4.3 Splitting experiment
For a given (standard, single block) EBM with a fixed number of layers, we ask how block splitting of this EBM into ff-EBM with multiple EB blocks affects training performance and Wall Clock (simulation) Time (WCT). We address this question with two different depths (L = and L = 12 layers in total) and various block sizes (bs), maintaining fixed total number of layers (e.g. for L= 6, I block of 6 layers, blocks of layers, etc.). Additionally, to ensure the fairest comparison in terms of WCTs across different splits, adopt the aforementioned TOL approach to execute the fixed-point dynamics Eq. (14) within each EB block. We display the results obtained on the CIFAR-10 task inside Table We observe that EP performance improves with smaller block sizes
Table 1: Validation accuracy and Wall Clock Time (WCT) obtained on CIFAR-10 by EP (Alg. 2 and ID models with different number of layers (T.) and block sizes ("bs"). 3 seeds are used.
cost
0.5
0.0
EP
ID Top-1 (%)
12
11 S7
Top-1 1 (%)
WCT
WCT
L=
Figure 5: Cosine similarity between EP and ID weight gradients sample r and associated label y in the same setting Fig. using the same color code to label the layers. We observe near-perfect alignment between EP and ID gradients. See Fig. 7 for a precise depiction of the model at use.
randomly selected
bs=6
89.2 #0.2
7:01
87.3
#0.4
6:51
bs=3
89.8
+0.2
5:17
89.3
#0.2
5:10
bs=2
90.1
10.1
3:57
90.0
10.1
4:05
L:
bs=4
89.4 #0.7
11:59
89.5 #0.2
8:28
bs=3
92.5
+0.1
7:33
92.0
10.1
4:16
bs=2
92.0
#0.2
3:14
91.5
#0.2
3:07
(reaching 90.1% and 92.5% for = and T. = 12 respectively with bs 2 and bs 3) with overall WCT reduction (up to s x4) while remaining on par with ID. This significant reduction in
8

---

WCT is due to the fact that inference time for ff-EBMs with DHN blocks by construction scales linearly with the number of blocks rather than supralinearly with the number of layers as has been empirically observed in the EP literature [Ernoult et al. 2019]. When instead using a fixed number of iterations to execute Eq. (14) inside EB blocks (Table in Appendix (A.5.3) EP performance maintained across all splits (90.1% and 92.5% for I, 6 and I. = 12 resp.) and is still on par with ID. However, there is no advantage in terms of WCTs in this case as the number of iterations is kept the same across all block splits and is much larger than necessary for smaller block sizes. Results for = 6 are consistent with the existing literature and those for L = 12 surpass EP state-of-the-an on CIFAR- 10 [Scellier et al., 2024 Laborieux and Zenke, 2022]. Overall these results suggest that: i) ff-EBM performance is agnostic to EB block sizes and are therefore flexible in design, ii) ff-EBMs much faster to simulate that EBM counterparts of equivalent depth.
4.4 Scaling experiment
We now consider ff-EBMs of fixed block size 2, and relatively small number of iterations. We train two models of depth (L = 12 and L 15) on CIFAR-100 and ImageNet32 by EP and ID and show the results obtained in Table 2 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 15-layer deep ff-EBM by EP exceeds state-of-the-ar performance on ImageNet32 by around 10% top-1 validation accuracy [Laborieux and Zenke, 2022] and by around 5% the best performance reported on this benchmark among all backprop alternatives [Höier et al. 2023]
Table 2: Validation accuracy and Wall Clock Time (WCT) obtained on CIFARI00 and ImageNet32 by EP and Autodiff on models with different number of layers (L) and a block size of 2 (bs=2). 3 seeds are used. We compare our results against best published results on ImageNet32 by EP [Laborieux and Zenke 2022) and against all backprop alternatives [Höier et al. 2023].
EP
ID
Top-1 I (%)
Top-5 (%)
WCT
Top- I (%)
Top-5 (%)
WCT
L=12
69.3 #0.2
89.9 10.5
4:33
69.2.0.1
90.0 10.2
4:16
CIFAR100
L=15
71.2
90,2+1.2
2:54
71.1
90.9 #0.1
2:44
L=12
44.7 #0.1
61:00
+0.1
65:23
44.7
#0.6
68.9
57:00
ImageNet32
[.=15
46.0 #0.1
70.0 #0.2
46:00
45.5 #0.1
69.0 #0.1
40:01
Laborieux and Zenke 2022
36.5
60.8
Høier et al. 20233
41.5
64.9
Discussion
EP literature. Ever since fixed-point iteration schemes were first proposed to facilitate EP experiments [Ernoult et al. 2019 Laborieux et al. [202]], there has been growing body of work assessing scalability of EP and algorithmic extensions on standard vision tasks. Most notably, Laborieux and Zenke [202) introduced a holomorphic version of EP where loss gradients are computed with adiabatic oscillations of the model by nudging in the complex plane, which was very recently extended to more general implicit models [Laborieux and Zenke, [023). Moving further towards physical implementations of EP, [Sellier et al. [002]] proposed 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 purely theoretical motivation and without clear algorithmic prescriptions, is that of Zach [02T]] where feedforward model learning is cast as a deeply nested optimization in which consecutive layers are tied by elemental pair-wise energy functions. This work more recently inspired the Dual Propagation algorithm [Hoier et al. 2023) J. Such a setting can be construed as a particular case of ff-EBM learning by EP where each EB block comprises single layer (U* 0) in- side Eq. 6 which, as we have shown, is tantamount to standard BP(see last paragraph of Section 3.3).
Forward-only learning beyond EP. Given that zeroth-order (7.)) optimization and "forward- forward" (FF) algorithms [Dellaferrera and Kreiman 2022 Hinton 2022 can be applied to any

---

model, and-like EP- compute a learning rule through multiple inference steps, one may wonder why it is important that our models should be aeergy-based, While mechanistically appealing for analog hardware [Oguz et al., 2023 Momeni et al. 2023 2024, Xue et al. 2024), these forward-only approaches do not match the performance of automatic differentiation on equivalent models, even they are roughly the same size as those studied in our work. On the one hand, weight perturbation [Fiete et al. 2007] (WP or "SPSA" [Spall] [998), yields unbiased yet noisy gradient estimates with variance scaling cubically with the model dimensionality (Ren et al. 2022]. resulting in significant gap in model performance compared to backprop, that can only be partially mitigated when using
heuristics [Silver et al.) 2021, Ren et al. [2022] Fournier et al.) 2023 Chen et al., 2023]. On the other hand FF algorithms, as learning heuristics, suffer from a lack of theoretical guarantees which may impact the resulting model performance.
Limitations and future work. Since our recipe advocates EP-BP chaining by construction, is fair to say that ff-EBM learning partially inherits the pitfalls of BP. Fortunately, noth- ing prevents feedforward modules inside fT-EBMs from being trained by any BP alternative to mitigate specific issues. For instance: BP can be parameterized by feedback weights to ob- viate 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 |Erouult et al. 2022]; computed via implicit nerward-mode differentiation (Hiratani et al., 2022. Fournier et al. 2023, Malladi et al., 2023]; local layer-wise self-supervised or super-
vised
loss functions can be used to prevent "backward locking" [Belilovsky et al., 2019 Ren
et al.
2022 Hinton 2022 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
Tie
Chip
(DRNs) - developed by Scellier [024] and strongly inspired by
Kendall et al. [2020] are exact models of idealized analog cir- cuits, trainable by EP. promising fast simulation times. As such, using DRNs as EB blocks inside ff-EBMs is an exciting re- Fig. Still, further work in this direction presents new challenges especially given device non-idealities which may affect the inference pathway, such as analog-to- digital and digital-to-analie
Global buffer
I/0 nteefae
search direction
Tile
Anoing
Dightl processo
noise [Rasch et al. 2023 Lammic let al. 2024]. Finally, considerable work is needed prove m-EBM further at scale on more difficult tasks (e.g. standard ImageNet), considerably deeper architectures, and moving be- yond 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.
Aralog
Buffer
Analog Processor
RUNJDO Aqpaeda kd
Concluding remarks and broader impact. We show that ff-EBMs constitute a novel framework for deep-learning in het- erogeneous hardware settings. We hope that the proposed algo- rithm can help to overcome the typical division between digital versus analog or BP versus BP-free algorithms and that the greater energy-effie energy-efficiency afforded by this framework provides pragmatic, near-term blueprint mitigating the dramatic carbon footprint of AI training [Strubell et al. 2020). While we are still a long way from fully analog training accelerators at commercial maturity, we believe this work offers an incre- mental and sustainable roadmap to gradually integrate analog, energy-based based computational primitives as they are developed into existing digital accelerators.
Figure 6: ff-EBMs as hierarchi- cal systems implementing EP at chip scale (adapted from et al. 2023)) using energy-based analog processors made up of resistors (green edges), diodes (in blue), voltage sources (in purple). ADCs and DACs (adapted from (Scellier 20244), digital processors, memory buffers, all of these being connected by digital buses (red lines).
10

---

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, Gregory Kollmer and Mohammed Fouda for helpful feedback on the manuscript as well as Benjamin Scellier for useful discussions last year which led to an alternative derivation of our main result (Appendix A.2.2 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 computationnl 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 Al" 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 Al which commercializes technologies based on brain-ing inspired learning algorithms, as well Constance Castres Saint-Martin for her unwavering support at the maternity hospital where most of this manuscript was written.
Author contributions
TN was responsible for implementation, architecture design, coding all algorithmic details and rumning training experiments, well discovery of criteria for stable convergence. TN also participated in writing relevant portions of this manuscript. ME designed the study, derived all theoretical results, dehugged and refactored the initial codebase and wrote most of the manuscript.
References
A. Agarwala and S.S. S. Schoenholz. Deep equilibrium networks are sensitive initialization statistics. In International Conference on Machine Learning, pages 136-160. PMLR, 2022.
M. Akrout, C. Wilson, P. Humphreys, T. Lillicrap, and D. B. Tweed. Deep learning without weight transport. Advances in neural information processing systems, 32, 2019.
L. B. Almeida. A learning rule for asynchronous perceptrons with feedback in combinatorial environment. In Proceedings, 1st First International Conference on Neural Networks, volume 2, pages 609-618. IEEE, 1987.
Ambrogio, Narayanan, A. Okazaki, A. Fasoli, Mackin, Hosokawa, A. Nomura, Yasuda, A. Chen, A. Friz, et al. An analog-ai chip for energy-efficient speech recognition and transcription. Nature, 620(7975):768-7777 767-77, 2023.
S. Bai, J. 7. Kolter, and V. Koltun. Deep equilibrium models. Advances in neural information processing systems, 32, 2019.
E. Belilovsky, M. Eickenberg, and E. Oyallon. Greedy layerwise learning can scale to imagenet. In International conference on machine learning, pages 583-593. PMLR, 2019.
M. Blondel, Q. Berthet, M. Cuturi, R. Frostig, S. Hoyer, F. Llinares- F. Pedregosa, and J.-P. Vert. Efficient and modular implicit differentiation. Advances in neural information processing systems, 35:5230-5242. 2022.
A. Chen, Y. Zhang, J. Jia, J. Diffenderfer, J. Liu, K. Parasyris, Y. Zhang, Z. Zhang, B. Kailkhura, and S. Liu. Deepzero: Scaling up zeroth-orde optimization for deep model training. arXiv preprint arXiv:2310.020055, 2023.
Chrabaszez, I. Loshchilov, and F. Hutter. A downsampled variant of imagenet as an alternative to the cifar datasets. arXiv preprint 17XXi:1707008//7 2017.
S. Cosemans, B. Verhoef, J. Doevenspeck, Papistas, Catthoor, P. Debacker, A. Mallik, and D. Verkest. Towards 10000tops/w dnn inference with analog in-memory computing-a circuit blueprint, device options and requirements. In 2019 IEEE International Electron Devices Meeting (IEDM), pages 22-2. IEEE, 2019.

---

G. Dellaferrera and G. Kreiman. Error-driven input modulation: Solving the credit assignment problem without backward pass. In International Conference on Machine Learning. pages 4937-4955. PMLR, 2022.
S. Dillavou, B. Beyer, M. Stern, M. Miskin, A. Liu, and D. Durian. Transistor-based self-learning networks. In APS March Meeting Abstracts, volume 2023, pages D07-006, 2023.
A. L. Dontchev, R. T. Rockafellar, and R. T. Rockafellar. Implicit functions and solution mappings: A view from variational analysis, volume 616. Springer, 2009.
M. Emoult, J. Grollier, D. Querlioz, Y. Bengio, and B. Scellier. Updates of equilibrium prop match gradients of backprop through time in an rnn with static input. Advances in neural information processing systems, 32, 2019.
M. M. Ernoult, F. Normandin, A. Moudgil, S. Spinney, E. Belilovsky, I. Rish, B. Richards, and
Y. Bengio. Towards scaling difference target propagation by learning backprop targets. In International Conference on Machine Learning, pages 9968-5987. PMLR, 2022.
I. R. Fiete, M. S. Fee, and H. S. Seung. Model of birdsong learning based on gradient estimation by dynamic perturbation of neural conductances. Journal of neurophysiology. 98(4):2038-205) 2007.
L. Fournier, S. Rivaud, E. Belilovsky, M. Eickenberg, and E. Oyallon. Can forward gradient match backpropagation? In International Conference on Machine Learning, pages 10249-10264, PMLR, 2023.
W. Haensch, Gokmen, and Puri. The next generation of deep learning hardware: Analog computing. Proceedings of the IEEE, 107(1))188-122, 2018.
G. Hinton. The forward-forward algorithm: Some preliminary investigations. arXiv preprint arXiv:2212.13345, 2022.
N. Hiratani, Y. Mehta, T. Lillicrap, and P. E. Latham. On the stability and scalability of node perturbation learning. Advances in Neural Information Processing Systems, 35:31929-31911, 2022.
R. Høier, D. Staudt, and C. Zach. Dual propagation: Accelerating contrastive hebbian learning with dyadic neurons. In International Conference on Machine Learning, 2023. URL https
S. Jain, H. Tsai, C.-T. Chen, R. Muralidhar, I. Boybat, M. M. Frank, S. Wožnnak, M. Stanisavljevic, P. Adusumilli, P. Narayanan, al. A heterogeneous and programmable compute-in-memory accelerator architecture for analog-ai using dense 2-d mesh. IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 31(1): (1:::114-227 2022.
J. Kendall, R. Pantone, K. Manickavasagam, Y. Bengio, and B. Scellier. Training end-to-end analog neural networks with equilibrium propagation. arXiv preprint vrXiv:2006.0195 2020.
D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980) 2014.
A. Krizhevsky. Learning multiple layers of features from tiny images. Technical report, 2009.
A. Laborieux and F. Zenke. Holomorphic equilibrium propagation computes exact gradients through finite size oscillations. Advances Neural Information Processing Systems, 15:12950-12963, 2022.
A. Laborieux and F. Zenke. Improving equilibrium propagation without weight symmetry through jacobian homeostasis. arXiv preprint arXiv:2309.022: 02, 2023.
A. Laborieux, M. Emoull, B. Scellier, Y. Bengio, J. Grollier, and D. Querlioz. Scaling equilibrium propagation to deep convnets by drastically reducing its gradient estimator bias. Frontiers in neuroscience, 15:633674, 2021.
12

---

C. Lammie, F. Ponzina, Y. Wang, J. Klein, M. Zapater. I. Boybat, A. Sebastian, G. Ansaloni, and
D. Atienza. Lionheart: A layer-based mapping framework for heterogeneous systems with analog in-memory computing tiles. arXiv preprint arXiv:2401.0P42 .99400, 2024.
W. Li, M. Manley, J. Read, A. Kaul, M. S. Bakir, and S. Yu. H3datten: Heterogeneous 3-d integrated hybrid analog and digital compute-ing -n-putor- accelerator for vision transformer seff-attention. IEEE Transactions Very Large Scale Integration (VLSI) Systems, 2023.
R. Liao, Y. Xiong, E. Fetaya, L. Zhang, K. Yoon, X. Pitkow, R. Urtasun, and R. Zemel. Reviving and improving recurrent 3082-3091. PMLR, 2018.
In International Conference on Machine Learning, pages
Liu, C. Mu, H. Jiang. Wang. Zhang, Lin, K. Zhou, Q. Liu, and C. Chen. Hardsea: Hybrid analog-reran clustering and digital- in in-memory computing accelerator for dynamic sparse self-attention in transformer. IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 2023.
I. Loshchilov and F. Hutter. Sgdr: Stochastic gradient descent with warm restarts, In International Conference on Learning Representations, 2017.
S. Malladi, T. Gao, F. Nichani, A. Damian, J. D. Lee, D. Chen, and S. Arora. Fine-tuning language models with just forward passes. Advances in Neural Information Processing Systems, 36:53038- 53075, 2023.
A. Momeni, B. Rahmani, M. Malléjac, P. Del Hougne, and R. Fleury. Buckpropagationn-ii free training of deep physical neural networks. Science, 382(6676): 1297-1303, 2023.
A. Momeni, B. Rahmani, B. Scellier, L. G. Wright, L. McMahon, C. Wanjura, Y. Li, A. Skalli, N. G. Berloff, T. Onodera, et al. Training of physical neural networks. arXiv preprint arXiv:2406.03372, 2024.
S. Nandakumar, M. Le Gallo, C. Piveteau, V. Joshi, G. Mariani, I. Boybat, G. Karunaratne, R. Khaddam-Aljameh, U. Egger, A. Petropoulos, et al. Mixed-precision deep learning based on computational memory. Frontiers in neuroscience, 14:406, 2020.
1. Oguz, J. Ke, Q. Weng, F. Yang, M. Yildiri,, N. U. Dinc, J.-L. Hsieh, C. Moser, and D. Psaltis. Forward forward forward training of an optical neural network. Optics Letters, 48(20):5249-522 5252, 2023.
F.J. J. Pineda. Generalization of back-propagation recurrent neural networks. Physical review letters, 59(19)2222,, 1987.
M. J. Rasch, C. Mackin, M. Le Gallo, A. Chen, A. Fasoli, F. Odermatt, N. Li. S. Nandakumar, P. Narayanan, H. Tsai, et al. Hardware- training for large-scale and diverse deep learning inference workloads using in-memory computing-based accelerators. Nature communications, 14 (1):5282, 2023.
M. Ren, S. Korbblith, R. Liao, and G. Hinton. Scaling forward gradient with local losses. arXiv preprint arXiv:2210.03300 2022.
O. Russakovsky, J. Deng, H. Su, Krause, Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein, et al. Imagenet large scale visual recognition challenge. International journal of computer vision, 115:211-255 2015.
B. Scellier. A deep learning theory for neural networks grounded in physics. arXiv preprint arX1v:2103.09985. 2021.
B. Scellier. A fast algorithm to simulate nonlinear resistive networks. 11X4 2024.
arXiv preprint
B. Scellier and Y. Bengio. Equilibrium propagation: Bridging the gap between energy-based models and backpropagation. Frontiers computational neuroscience, 11:24, 2017.
B. Scellier and Y. Bengio. Equivalence of equilibrium propagation and recurrent backpropagation.
Neural computation, 1222):3
13

---

B. Scellier, S. Mishra, Y. Bengio, and Y. Ollivier. Agnostic Sehysics-driver deep learning. arXiv preprint arXiv:2205.15022 2022.
B. Scellier, M. Ernoult, Kendall, and S. Kumar. Energy-based learning algorithms for analog computing: a comparative study, Advances in Neural Information Processing Systems, 36, 2024.
A. Sebastian, M. Le Gallo, R. Khaddam-Aljameh and E. Eleftheriou. Memory devices and applica- tions for in-memory computing. Nature nanotechnology, 15(7):529-544 2020.
D. Silver, A. Goyal, I. Danihelka, M. Hessel, and H. van Hasselt. Learning by directional gradient descent. International Conference on Learning Representations, 2021.
C. Spall. Implementation of the simultaneous perturbation algorithm for stochastic optimization. IEEE Transactions on aerospace and electronic systems, 34(3):817-823 1998.
K. Spoon, H. Tsai, A. Chen, M. Rasch, S. Ambrogio, C. Mackin, A. Fasoli, A. M. Friz, P. Narayanan, M. Stanisavijevic, et al. Toward software-equivalent accuracy on transformer-based deep neural networks with analog memory devices. Frontiers in Computational Neuroscience, 15:675741, 2021.
M. Stern and A. Murugan. Learning without neurons in physical systems. Anmual Review of Condensed Matter Physics, 14:417-441, 2023.
M. Stern, S. Dillavou, Jayaraman, D. Durian, and A. J. Liu. Physical learning of power-efficient solutions. arXiv preprint arXiv:2310.10437, 2023.
E. Strubell, A. Ganesh, and A. McCallum. Energy and policy considerations for modern deep learning research. In Proceedings of the AAAI conference on artificial intelligence, volume 34, pages 13693-13696 2020.
N. C. Thompson, K. Greenewald. K. Lee, and G. F. Manso. The computational limits of deep learning. arXiv preprint arXiv:2007.05556 2020.
Z. Wang, P. S. Nalla, G. Krishnan, R. V. Joshi, N. C. Cady, D. Fan, J.-s. Seo, and Y. Cao. Digital- assisted analog in-memory computing with rram devices. In 2023 International VLSI Symposium on Technology, Systems and Applications (VLSI-TSAAVLSSI-DATT), pages 1-4. IEEE, 2023.
Z. Xue, T. Zhou, Z. Xu, S. Yu, Q. Dai, and L. Fang. Fully forward mode training for optical neural networks. Nature, 632(8024):280 286, 2024.
S.-i. Yi, J. D. Kendall, S. Williams, and S. Kumar. Activity-difference training of deep neural networks using memristor crossbars. Nature Electronics, 6(1):45-51, 2023..
C. Zach. Bilevel programs meet deep learning: A uniffying view on inference learning methods. arXiv preprint arXiv:2105.072:: 2021.
Z. Zhang and M. Brand. Convergent block coordinate descent for training tikhonov regularized deep neural networks. Advances in Neural Information Processing Systems, 30, 2017.
N. Zucchet and J. Sacramento. Beyond backpropagation: implicit gradients for bilevel optimization. arXiv preprint arXiv:2205.03076, 2022.
14

---

A
Appendix
Contents
A.1
Model details
16
A.1.1
Feedforward tied EBMs (ft-EBMs)
16
A. 1.2
Feedforward nets as a special case
16
A.1.3
Equilibrium computation
16
A.2
Main theoretical result
19
A.2.J
Proof of Theorem
19 ................................
A.2.2
An alternative proof of Theorem 3.1
23
A.3
Resulting algorithms
26
A.3.1
Explicit BP-EP chaining
26
A.3.2
Recovering backprop through feedforward nets as a special case
26
A.3.3
Detailed implementation of the implicit BP-EP chaining algorithm (Alg. 2
27
A.4
Statie gradient analysis
29
A.4.1
Algorithmic baselines
29
A.4.2
Proof of Theorem
31
A.4.3
Details about Fig.
33
A.5
Experimental Details
A.5.1
Datasets
34
A.5.2
Data preprocessing
34
A.5.3
Simulation details
34
15

---

A.1 Model details
A.1.1 Feedforward-tied EBMs (ft-EBMs)
We first formally define Feedforward-tied Energy-based Models (ff-EBMs) with precise assumptions on the energy-base based and feedforward blocks.
Definition A.1 (ff-EBMs). A Feedforward-tied Energy-based Model (ff-EBM) of size N comprises twice difserentiable feedforward mapping Ft,.... F FF and N - twice differentiable energy
functions E EN with respect all their variables. For a given the inference procedure reads as:
x = 1it (s! ,,,,,,
$ == arg min (s, =I......N-
(17)
O₂ = FNN (s ,
Finally, assume that vk = .....N....... 1, ( 0h ,w/) is invertible.
A.1.2 Feedforward nets as special case
We show that when energy-based blocks comprise single layer only, the ff-EBM becomes purely feedforward.
Lemma A.1. We consider fJ-BM per Def. (A.1) where the energy functions F have the form:
E* (s 0% ,0,,x*))=00 G
x + U* (U*,,00
(18)
We assume that U f = N...N.-- - VG(s) is invertible and we denote a == VG-1 T Then, the resulting model is feedfoowwrd model described by the following recursive equations:
x
lotk (s s. 0, == F. (5 (oNi-1 w
of =
-1
(19)
Proof of Lemma
A. Let k E [1,N N 1]. By definition of S and +::
V 1 gh " ;/,)x) =-== =0 VG* - 2 + V U* ph) -0
- a (s v (stt 0h))
(20)
Therefore Eq. 19 is immediately obtained from Eq. 20 with U
A.1.3
Equilibrium computation
For single block. form:
As mentioned in Section 3.1 the energy function of the k EB block has the
E* (s -
* Us*
(21)
where a is the output the preceding feedforward block. For a given choice of a continuously invertible activation function, G is defined as:
dimm("") 2
(s*) = E
0, 0 (u:)du; such that (s*), o, ( =i=
dim
(22)
16

---

To be more explicit and as we did previously, we we-write the augmented energy-function which encompasses both the k EB block and the feedforward module that precedes
Es 6t, st- ww) = E (s",0" 0" FR (s)
ww
(23)
Deep Hopfield Nets (DHNs) as EB blocks. In our experiments, we used weight matrices of the form:
0f
0
0
0
0
@h
0
02
(24)
whereby each layer e is only connected to its adjacent neighbors. Therefore, fully connected and convolutional DHNs with layers have an energy function of the form:
U 00)
1 KT
0h g ===- - 8+1 0f
1
Σ
(25)
2
(6 0h) =
2
1
(0k * s =
Σ
6+1
(0) * x
(26)
Synchronous fixed-point iteration. We showed that when G is chosen such that VG = o for some activation function a, then the steady state of the kth block reads:
: Vi
(27)
which justifies the following fixed fixed-point iteration scheme, when the block is influenced by some error signal 8s with nudging strength 8:
18.41 t - (s 8*) F B85*)
(28)
The dynamics prescribed by Eq. 28 are also used for the inference phase with B 0. further refine Eq. (28). let us re-write Eq. (28 with layer index é where E E [1, Lx] with L3 being the number of layers in the kth block, and replacing b² by its explicit expression:
v = 1 .... Ta L1 <- Ft (s--1, --k - V Uk (si F B55+)
(29)
As done in past EP works Ernoult al. 2019
Laborieux al. 2021
Laborieux and Zenke 2022
2023 Scellier et al. 2024 and for notational convenience, we introduce the primitive function of the
k block as:
** (s 0 st (w))) = st Fk (st-11 w - (**
t
(30)
such that Eq. 29 re-writes:
V = 1 Ta s1 t (si 0r 0+ st-1,, F 3355
(31)
Eq. 31
depicts synchromous scheme where all layers simultancously updated at each timestep.
17

---

Asynchronous fixed-point iteration. Another possible scheme, employed by Scellier et al. [2024).
instead prescribes asynchronously update odd and even layers and was shown speed up conver- gence in practice:
odd (1,.., Lx)
2,,3,44 + a
V V
st L 0h. of .o- wt) F 3305
(32)
even & E (1,-.. La}:
18,411 t (V 2
Oh st wu)) Bösh
We formally depict this procedure as the subroutine Asynchronous inside Alg. 3 In practice, we observe that was more practical use a fixed number of iterations rather than using a convergence criterion with fixed threshold.
Algorithm 3 Asynchronous (for all blocks until penultimate)
Input: T, 8%, wh, s-- B, 8s
k-1
&k Output: $
1: k (0
2: for 1...=
3: oodd& l fl, Lx} $
(V (st sh- .wa)
-555*)
4: Veven l (1, L&] &t
V
$. Ah
sh-1 ,w() Bösh
5: end for
Resulting ff-EBM inference algorithm. With the aforementioned details in hand, we re-write the inference algorithm Alg. presented in the main as a Forward subroutine.
Algorithm 4 Forward
Input: T, x. W {0] w/
Output: s s",... ... N-1 or ò depending on the context
1: S
2: for = l....N.
do s* +- Asynchronous (T,0*, 0h wh 1
3:
Alg
4: end for
5: ô t FN (s
18

---

A.2
Main theoretical result
A.2.1 Proof of Theorem 3.1
The proof of Theorem 3.1 is structured as follows:
We directly solve the multilevel problem optimization defined inside Eq. 9 using Lagrangian- approach (Lemma A.2). 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.3).
We demonstrate Theorem 3.1 by combining Lemma A.2 and Lemma
A.3
Finally, we highlight that when a ff-EBM feedforward net (Lemma A. then the proposed algorithm reduces to BP (Corollary
A.I. A.II.
Lemma A.2 (Lagrangian- based approach). Assuming a ff-EBM (Def.
A.7 we denote
sh .8N-1,0 Ô0 the states computed during the forward pass as prescribed by Eq. 117. Then, the gradients of the objective function C == ((0( NN11), ),,9) as defined in the multilevel optimization problem (Eq. ( where it is assumed that E is differentiable, read:
d 022 (sN- W -)e((o. v), do V4, (ot, Hs sh-1 wh) 2
vk = 1....N 1,
(33)
duhc = E*, gh sh .wt) 2
Wk 1... N-1
where No, .NN-1 satisfy the following conditions:
V.N Vxx-1 [(((s (6(s)_1)),20) ).8) + VIEN N (s 1.5-1,NN-1
2- w N-1) A 1N =0
Vk N ........ 1:
(34)
Vi. Ek+1 (st+1, gk++1 sh w+1) Aht1 ViE* (s 0 s wwt) 2
Proof of Lemma A.2
Denoting $ == (s',...,, ".......---"-" N the state variables of the energy-based
blocks, =
A 1 the Lagrangian multipliers associated with each of these variables. [ WT, ON 1, - the energy-based and feedforward parameters and N(1N == FNN w w- the logits, the Lagrangian of the multilevel optimization problem as defined in
W
Eq. reads:
N
C(s, = e +
Σ TT V (s gh, sh-1
sk-1, ,,a)), 50 ==
(35)
k=1
Writing the associated Karush- Tucker (KKT) conditions O S2,2L(s,, W) = 0 satisfied by optimal states and Lagrangian multipliers s, and X we get:
V Vi 0h k-1 , =0 VK Ik=1... 1
(36)
V Nー ((0(s N-1 )N + (1 NN N N2 N-1) X
NN-1 St w()) 2 =0 Vk N - -2 2,...
(37)
E k+1 (s gh+1 sti w+t) xh+1 ViEN (st th
(38)
Eq. 36) governs the bottom-up block-wise relaxation procedure (as depicted Alg D. while Eq. (37
and Eq. (38) governs error propagation in the last block and previous blocks respectively. Given
8, and X, by Eq. 366 Eq. 38 the total derivative of the loss function with respect the model
parameters read:
19

---

N-1
dw (w dw l (ó., +
Σ
VIET
El t 0t k-1 w
kml
(Eq.
(6)
dw Lox,, W)
ddss 0 01 de; W) +dw
02h(82, Ao, W)
+0gL(8x; 2x,W) W)
(Eq. 36
(Eq. 37 38
= A., W)
(39)
More precisely, applying Eq. 339 to the feedforward and energy-based block parameters yields:
devi(óo,, 3) F2p (s w W Vi((o.,
dor 4) E
Ak = 1....N-
dw V4 (8
wh X Ak -
Lemma A.3 (Computing Lagrangian multipliers by EP). Under the same hypothesis Lemma A.2
k. we define the nudged state of block ki denoted as Sy implicitly through V 3. 0$ with:
x 8s B)
( 8", x. ".. B) (s" " "0000""))) 35* 88
(40)
Defining ( (==)) I.....NN11 recursively as:
80 = dg (vzBatl
(8g ,6 ,f ((,.ar1))
k+1 ,0+1
4k = 1...
N-2,
3=0
(41)
then we have:
x dg (3 1==0
1
(42)
where (x))==1..N-1 are the Lagrangian multipliers associated to the multilevel optimization problem defined Eq. 49)
Proof of Lemma
A.3
We prove this result by backward induction on
Initialization (k
1). By definition, $ N-1 satisfies:
v x (ô., y) + V N- (s NN -N 1 N-2 (-1,
(43)
Differentiating Eq.
43 with respect and evaluating the resulting expression B obtain:
(b., + ViEN-1 (s NN1 N N-1) d3s N- - = 0
(44)
Substracting out Eq.
37 defining the Lagrangian multiplier AN-1 and Eq. 44pp. obtain:
ViEN (N N-1 NN N-2 2 w-1) (das -1 - X 1)
(45)
By invertibility of N- (1 -N . w-1, therefore have that:
N
d3s N 1 3
(46)
20

---

Backward induction step (k + + k). Let us assume that X 10=0. We want to prove that x dgs0la=o Again, st satisfies by definition:
3ös + V (5 0h sh w")) O, 55 = d3 (Volitt) (tt 10tt1,,00000 syt w)t)))
3=-
(47)
the one hand, proceeding as for the initialization step, differentiating Eq. with respect to and taking 3 0 yields:
85 + VIE st- dosøll=0 = O.
(48)
On the other hand, note that 5s rewrites
8s = d3 (vgEtti k+1 ghtt s w(ttt)))
1=0
Vit (sit1 gh- st wh1 ds dstt11 Vi.a (stt1 gh-1 st wht)) \ht1
3=0
(49)
where we used at the last step the recursion hypothesis. Therefore combining Eq. we get:
(48 and Eq.
449
Vi (
gh+1
st wttt
X VIE sh ,w,) dastla=0
(50)
Finally, we substract out Eq.
38 and Eq.
50 to obtain:
V
0h
sh wh) (das 21) = 0.
(51)
We conclude again by invertibility of VG (sh, 8h sh .w²) that X dasglamo
Theorem A.1 (Formal). Assuming a model of the form Eq.
(5) we denote s x1,
N-1 Ô2 the
states computed during the forward pass as prescribed by Alg. We define the nudged state of block & k, denoted as gg implicitly through V N 8% rkk 8t (3)
O with:
Fdh (s", 0hh wht 8 ,9) = (s 0h xx) 3s* 85
(52)
Denoting hs and & the error signals computed the input of the feedforward block Fk and of the energy- mergy-based block E respectively, Ig and 9 the gradients of the loss function:
Vk =k=1.... ,N
- 1 gor = Co 1k=1...N Dwc = Cu
(53)
then the following chain rule applies:
N = V.N (ow-i((0.,))) IWN = 2zF N (5.2-1 (3"--1,w-)) vie(ō.,
(54)
Vk 1.... N. 1:
/ Art = dg (Vz (s 0², 00, A=0² 9o- = dg (2 8h 0x, 11 D ph (s ( Art, Hut = O2 r (s+- w
3=0
(55)
T Art
Proof of Theorem A.7 Combining Lemma A.2 and Lemma A.3 the following chain rule computes loss gradients correctly:
21

---

58 N-1 = (( INN = =2FFN (sM-- ((")
T. Vi((ó.:y)
(56)
Vk- = ....N- -1
A = =d3 (V (5tg 0h st- - ( '=00 Jet 2i (s 0h
(4-1,wh) d3st 18==
O+ Ei (st, 0h, sh dgst
(57)
Therefore to conclude the proof, we need show that vk 1,... ,NN 1:
d3 (V2Ê* (s), 0r s ((a))) 1--0 = F₁ (s T. d3 (V3 0h ,000,2,))) 22)) (58)
Vi,2 (5 0 st ww/)
d3 (V2 (5 0t x 1
(59)
(s, 0h s wh) d=stو6ep = 2² ,wk) T_ - d3 (V3 0h 23)1
,,,,,))))
(60)
Let E [1, - 1). We prove Eq. 558 as:
dg Vvz ($ 0h 1 ) 3=0
= dg (V E* (5 0h F* (()), مدا
= F* (s+ T
dg (V3E*(s33 0h x 8=0
Eq. 59 can be obtained as:
0t, 5. ,w") d30s13=0 = dg (0), 0t ,5 ,
1=0
= da (V2F2 0² ,bti) 2=0
Finally and similarly, we have:
Vi, Ed (s ," -h ,w") dassle=o = d3 (V.E" 0h (._l_="))
3=0
dg (V F² 0F Fat (s- ,wn))) محا0
=dp (B2F (5. ,w/) V (3 @ ,x+)) 2=0
= (sf-2 w0)) T- dg (V gh. x o=0
22

---

A.2.2 An alternative proof of Theorem 3.1
An energy function for ft-EBMs? While it is clear that the energy function of ff-EBM is not E* (which would correspond in this case to the "single block" standard case), one may
F Σ wonder if:
ff-EBM inference (Alg. can still be described as the minimization of some energy
function?
Therefore, if Theorem 3.1 can be derived by directly applying EP to this energy function?
We show below that this is indeed the case. We follow Zach 2021 denoting 8 =
(s and W (W,,"* WN-1I. by picking the following energy function:
N-
F F x) = Σ
E (s": W 11-1)
+ Eg (shtl, wk+1 ,,, - EgEtt (s++1, wh+1 .t)]] (s (f)
+ -
N- (s WN-,, , 3) sē(s y. W
(61)
where we locally redefine x as the concatenation of all block inputs, i.e. (xT s. ..... NN )T and with S. Fi,(s,, v, 0) = O. In Lemma A.4 we show that the free steady state of the above energy function (s, obtained with by the ff-EBM inference scheme (Alg. D
(s+ ..... .sN-1 implicitly defined through
inside Eq. (6T) indeed corresponds to the states computed
Lemma AA. Let El.... EN be the bisck-wise energy functions of a ff-EBM defined per Def. A.1 Assume 8, implicitly defined through F(,, W, B = 0) = 0 where is defined by Eq. 161 Then:
Vk 1, 1: V F (s W
-
(62)
Proof of Lemma A.4 For k
N - 1, the stationarity condition FF-1F W, x., B) yields:
N- (s W N N - +0 0.
(63)
Then, for any
S N V F(..: 3) yields:
k W - +
+++1 (1 w Eb++ (st+1 Wh+1 sk)
(64)
Eq. 63 and Eq.
63 indeed correspond to ff-EBM inference depicted inside Alg. 1
The EP fundamental Lemma. result below inside Lemma A.3
For seff-containedness of this paper, we restate the fundamental EP
Lemma A.5 (Scellier 2021)). Let F(s, W. 3) be a twice differentiable function of the three variables W and B. For fixed W, and B, let 3 be point that satisfies the stationarity condition:
F W, x. f) 0
(65)
and suppose that ViF(sg. W. 3) is invertible. Then, in the neighborhood of this point, we can define a continuously differentiable function (x. W, W.B) sg such that Eq. (65) holds for any (a, W. B) in this neighborhood. Furthermore, we have the following identity:
dw (V F(p W. 3)) = d (V2 W, R))
(66)
In particular, Eq. (66) may be evaluated with F E + 38 at 3 = 0 to yield the EP learning rule,
denoting C : {(...)) v) (Scellier and Bengio 2017
dwC dg (V,F(ss, W.z, 9)) 1==0
(67)
23

---

Theorem 3.1 as direct application of EP. Now we know Eq. 61 defines a valid energy function for ff-EBMs and with Lemma A.5 in hand, we are ready to apply EP directly this energy function. We rewrite below the block-wise wise free energy functions at use inside Theorem 3.1 and used in practice
inside Alg. 2 to nudge block of energy E given some top-down error signal 8*:
(s k W 1 1,05 () == (s Wk. o. +Bs"" -884,
os = dæ Vz 5/. Wk+1, st st)) A=
if < N else V 1 y, WAY
(68)
In Theorem A.2 we show that the direct application of Lemma A.5 to as defined inside Eq. 61 yields the same gradient formula for each parameter W and the same nudged block states as those prescribed by Theorem 3.1 for sufficiently small B.
Theorem A.2 (Informal). Let F be defined in Eq. (67) satisfying the same assumptions in For fixed W, x and 3, let 8g satisfy the stationarity condition:
Lemma A.3
V1J W, x, 0.
Then, we have:
dwiC= = d3 (V2E" Wk: sk s4-11)
VIF* W 1k "05=, B) 0(52) B=0"
(69)
Proof of TThorrmm..... A.2
Lemma A.5 yielding:
dw:C = ds (VwiF(sp, w,,,
proving Eq. 669 amounts to show that:
dæ (V Wur W, (2,,,))) 3=0 = d3 (V2 (s Wh 1 ..4-11
(70)
0=0
V V, Br W k-1 1 = 0(33)
(71)
On the one hand, we have:
V F W,,,3) V2 (s) Se W
+ V22 Eth W
E W ( ..-1)
(35 - s s2-2)
(72)
For notational convenience, we define A(3)
V2
V E* - Wk s st-1))
W ,
Note that A(B 0) = 0. Differentiating Eq. 72
with respect and taking 8 0 yields:
d3 ((se B)) 1==0 = dg (v (s Wh , s4-11 l==0
+ dggA(3))
- 11 + A(O)
- s411)
=0
dg (V2E*(sh) W . st-12 l==0,
which proves Eq. 70.
On the other hand, the stationary condition V V,rF(ss, W, x, 3) on the last block (k = N - 1) yields:
V VIENI1 N- N-
,8, /V 1(( u W =0
V, NN E N- 3 N-2 + 3V 11,,- 0(82)
V,EN- -1..NNN V, N-2
2 - 033).
(73)
24

---

For previous blocks, i.c. lo < N 1. we have:
VetF () =0
V ( W* ) V3 Fo (sttt WWhh,, ... - Fz (s5 Whtl ,04)= = 0
V V. 1 d3 (z (s whtt (4), 0- 0
V W + 8s 3) (
(74)
Altogether, Eq.
73
and Eq. (4) finishes to prove Eq.
25

---

A.3 Resulting algorithms
A.3.1 Explicit BP-EP chaining
We presented in Alg. 2 a "pure" EP algorithm where the BP-EP gradient chaining is implicit. We show below, inside Alg. the use of BP through feedforward modules is explicit and which is the direct implementation of Theorem A.I We also show the resulting algorithm when the ff-EBM reduces feedforward net (Lemma A.I inside Alg. 7 highlight in blue the statements which differ from the general case
alternative implementation (equivalent in the limit 8 - 0) where
presented inside Alg.
Algorithm 5 Explicit BP-EP gradient chaining (Theorem 3.1
1: 88, N +- V ((., v). Www ((O..
Single backprop step
2: for N 1... do
3: $3 t Optim E* (s, ,,,,, ) + B8² 88]
b EP through E
4:
8 t Optim E (s, 0h x - =3st= ós]
5:
9ø - 23 (V2E* (8g, 0h
x - 2 -B 0h 241))
6:
Ar < 23 (VzEt (83, 1,, Hh x v3F (s-,, B, 0h .x1))
7:
INt t 0² (sh-1 1o4-1,w w
Ar
D Explicit BP through fob
8:
ōs + Oj pkk (st-1, Ar
9:
end for
A.3.2 Recovering backprop through feedforward nets as special case
Corollary A.1. Under the same hypothesis Theorem A.7 and Lemma A. then the following chain rule applies to compute error signals backward from the output layer:
N-1 = V V.N-il(d... v), Ar* ox @ 8s
IWN = VaN [
(75)
dsk-1= 0 Fh ( T, Art, Het O2 (s/-1 ,w+)T Ar
Proof of Corollary A.7
Let k E [1, - 1]. As we can directly apply Theorem A. here, proving the result simply boils down to showing that:
Ar = (x O 8s
(76)
First, we notice that when E of the form of Eq. (18). then Ar reads as:
Ar da (V 0h at)) 9=0 = de (5/) 3=D
(77)
of satisfies, by definition and when U 0:
o (55) x² + 38s
(z )
(78)
Combining Eq. and Eq. 78 yields Eq. 76 and therefore, along with Theorem A.I the
chain-rule Eq. 75).
We showcase in Alg. 6 and Alg. 7 the resulting algorithms implicit and explicit BP-EP chaining respectively, with lines in blue highlighting differences with the general algorithm Alg. 2
26

---

Algorithm 6 Implicit BP-EP gradient chaining with U
88, JNN t- x 1-x-1((0 f(0x, v),
(ww((on.))
Single backprop step
2: for k N- - -1. 1.... ...1. do
3: 83, B-B f (xf Bds*), (2 + 38s*)
P. EP through E
4: Iw t 28
(sp s sh1 ,w²) - V,E" 3 0º s--1 ,wa))
> i-BP through F
5: òs +
1
V3 13 0h, st w V3 (s 0h, s - .w,)
23
6: end for
Algorithm 7 Explicit BP-EP gradient chaining with U
1: 88, IN + V x- f(ō,, v), Van ((O,,))
D Single backprop step
2: for N 1. 1... do
3: Ar t- 1 (o [ 36s*) - -o (x+ + ))85
4: INt & O₂ Fk (st-1, w Ar
D Explicit BP through Fk
5: ōs 01 Fk
Ar
6: end for
A.3.3 Detailed implementation of the implicit BP-EP chaining algorithm (Alg.
Nudging the last block. From looking at the procedure prescribed by Theorem and algorithms thereof (Alg. 2 Alg. 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:
J- (s 100- -N x,, 1 = 1N -N ,
Bs + VAN re(0.. v), (79)
where d = PN around s.N-1 instead of the cost function itself. This is, however, in contrast with most EP implemen- tations where the nudging force acting upon the EB block usually elastic, Lee the nudging depends on the current state of the EB block. More precisely, instead of using Eq. (79). we instead use:
is constant. Up to constant, Eq. (80) uses the cost function linearized
JN J ON 2
.N ON-1,,NN-1 -NNN111 3)) E. -NN-1( ON11
3E(F (sN-1,,~")
(80)
This results in the following asynchronous fixed-point dynamics for the last block:
Yodd l (1,..., La}: 8,1 +- a
(V
(oe.. ( k-1 ww²) F BV (
Veven é E {1,... [b):
V
T sk ( 0h M , wa)) F BV ((s
2,14111
The resulting Asynchronous subroutine, applying for the last block, is depicted inside Alg.
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 readout in combination with the cross entropy loss function. Using our formalism, our readout is simply the last feedforward transformation used inside namely FN (-, N).
Detailed implicit EP-BP chaining algorithm. We provide detailed implementation of our algorithm presented the main (Alg. 2 in Alg. IT usually done for EP experiments, we always perform "free phase" to initialize the block states (Forward subroutine, Alg. 4). Then, two
27

---

Algorithm 8 Asynchronous (for last block)
Input: T, 1NN-1 w N-1 W sh-1 8, & (cost function), y
N-1 Output: B
1: s N- 0
2: for 1.... do
3:
A lodd l E { ... LiN:
4:
$1.9
N-1 to o
V i x
N-1 ,N
-N--2 3,--,,N-11 w- -38 op- eN-i (s _ v)
5: eevn l (1, .,. }:
6: $6 a
v x- $
BV 10 (s ,,0, v) N-1 aN-1 ,N sN w
7:
end for
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 (os" -2). Then, the same procedure recursed backward through blocks (Alg. TO]. until reaching first block.
Algorithm 9 BlockGradient (for last block)
Input: T, N-2 0V-1 N-1, W.,,N .. E. Output:ds
1: N- --1 Asynchronous
(T, AN w W.W ,B,,,)
Alg.
N-1 2: 8 -B
+- Asynchronous
(T, 0N w -1 W. -3,e,y)
3: Jw
V
C(FN
N 1 w
VoN-i 8(FN
N-1 -
4: -0 28 V2EN-IN (8) (SN-1 NN-1,
1, N-2 NN-2,wN-11 V,EN -B
AN .s N-2 W. N22)
5: N-N
V, (5 N-1 ON-
,84
N-2 w
N-1)
V. (1
N- N-1 NN NN- w2
28
6: N-2 K
VgEN-1 N-1 ON-
N
VgEN - N-1 ON N
N-2 ..
23
Algorithm 10 BlockGradient (for all blocks until penultimate)
Input: T, - 8k wh B, 8s
Output:ds k-
1: s Asynchronous 0h w* w (,,5s)
A
2: sk -B Asynchronous (T, th w -5,5s)
3: 9o
< 23
1
12
wh) - V Fo (s 0h 8 k-1 w.))
4: lank t 28
1
V v
w V V,E* (s
Qh she wh
5: &-h-1 --1 + 23
V 3ヤ
(
st w V Eu (s 0h sk w())))
Algorithm 11 Detailed implicit BP-EP gradient chaining
. NNN1 1: .
Forward (Tireer W)
Alg
2: 8s t -BlocGrrddient (Tuudge: sN N-2 22-2 -1 N-1 .3,,
b Alg
3: for k
2......1 do
4: ös + BlockGradient (Tmugges si wh 33,8s) 88)
D Alg
10
5: end for
28

---

A.4 Static gradient analysis
Important foreword. The whole subsection dedicated to an important tool when developing code for EP research. While EP is agnostic to how the steady states are obtained the EP theory only prescribes they are energy minimizers they can be obtained in practice (i.c. in simulations) through fixed-point iteration schemes (see Appendix I.IIII. The below formally defines the computational graph spanned by these schemes and abstract them away into a transition function K and defines three different techniques to compute gradients
this graph: Automatic Differentiation (AD, Prop. A.T). Implicit Differentiation (ID, Def. ( or Equilibrium Propagation (EP, Def. Aef. After defining each of these algorithms formally, will state and demonstrate equivalence between EP and ID (Theorem A.3)
which we test numerically and relied upon for the development of our codebase.
A.4.1 Algorithmic baselines
Definition of the computational graph being optimized. We abstract fixed exed-point iteration dynam- ics away into a kernel function K which, given some block state sf yields s $
Definition A.2 (Form of the computational graph through equilibrium).
vk 1, ... .N 1, Vt = 1, T x. n, of K(s W W x² 1)- C= C(F (s (.r C(FN(sN
= (= .)
(81)
Note that we emphasize, through the W = W* notation, that the parameters W are shared across all timesteps t = 1.... ,7 This will help define loss gradient with respect to W further below, i.e. how much 4 contributes time t to changing the loss C. The total contribution of W reads as the sum the elemental contributions of all W This intuition is more precisely illustrated further below. Given the computational graph defined in Def. A.2 can now formally define the Automatic Differentiation (AD) baseline.
Automatic Differentiation (AD).
Our goal is to compute:
AD (T)
with: (1)
Ow C
(82)
In plain words, W (t) denotes the loss gradient for parameter W truncated the tth step moving backward in time. We formally define below Automatic Differentiation (AD).
Proposition A.1 (Automatic Differentiation (AD)). The gradients AD (1) can be computed using the following recursive equations:
Vk N 1 1..
8sh 8x k+1 if N else Vrē(s,-1) .y)
8xt = O,
Iw (0) = 0
Vt = 1, .
{ 5 0 ((s W -² - -55/11
GAD (1) AD (1 1) 8g (ss
Wh 2 s 1-1 TTTT...
5st
(83)
dxt Sx 1 + 02K( ( WK xt s 65t
Proof of Prop. A.I This is straightforward application of the chain rule applied to Eq. (81)
Implicit Differentiation (ID). We define the steady state of block k, which we denote & as the fixed point of Eq. (8T). With this notation hand, we can define Implicit Differentiation (ID) this setting.
29

---

Definition A.3 (Implicit Differentiation (ID)). Denoting s the fixed point of Eq. 81
inside block k.
we define Implicit Differentiation (ID) through the following recursive equations:
4k = N 1:
5 6 ++1 if < N - else Vil(o,-
S 0, Ow* (0) 0
ID
1I 1, 1.. T
t 0 ( W r² = - -5s k-1
ID (t) = fw (t 1)
.ID
1) + 02 (s Wh at 8 -
(84)
dx, Sz,_1 0gK(s, W ,
I
We are ready to state simple algorithmic equivalence between ID and AD, which we built upon for our implementation Alg. 12
Corollary A.2 (Equivalence of ID and AD). Assuming that:
Vk = 1, N 1, Vt 1......
.
(85)
where s denotes the fixed-point of Eq. 81 then automatic differentiation (Prop. A.I and implicit differentiation (Def. AA3) are equivalent, namely:
Vk 1, N - 1, Vt
9 (t) 9 (t)
ID
(86)
Proof of Corollary A.2 This is a straightforward application of the definition of AD (Prop. A. along with the hypothesis made inside Corollary A.2
Resulting implementation of ID. We describe our implementation of ID inside Alg. 12 First, we relax all blocks sequentially equilibrium following Alg. 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 the previous step, we re-execute the same procedure (still with Alg. 44. this time tracking gradients and using Tuudge steps fixed-point iteration steps for each block. Then, we use automatic differentiation backpropagate through the last Tuudge steps for each block, namely backpropagating, backward in time, through equilibrium.
Algorithm 12 Our implementation of ID
1: Without tracking gradients:
e.g. with torch, no ho_grad (
2: ...... A. N +- Forward (Tinee I, W)
D Alg.
3: Initialize block states at s,.... N
4: ô, + Forwardd (Inudge, x, W)
> This time gradients are tracked
5: C "+E(
6: Backpropagate C backward through the last nudge steps for each block b e.g. C. backward ()
An important note about this implementation of ID. Note that this is not a standard implementa- tion of ID and it may be surprising at first glance to implement ID as AD, thereby loosing the constant (1) memory cost of with respect to the length the computational graph. Instead, the memory cost of Alg. 12 is O(N 1)r) However,
goal is not much to optimize for memory usage (as in the context of Deep Equilibrium Models [Bai et al.] 2019)) but to code an algorithmic baseline which we know to be equivalent to EP. Lastly, note that this implementation of ID is also known as Recurrent Backprop (RBP, [Almedda| 1987 Pineda 1987) ) or Von-Neumann RBP [Liao et al. 2018
and that ID generally comes in many more algorithmic flavors [Blondel al. 2022
We are not accounting for the spatial depth (L) of the computational graph in this cost. In this case, standard ID would have memory cost OLL) and our implementation inside Alg. OOL((N 1)7).
30

---

A.4.2 Proof of Theorem 4.1
In order to state formal equivalence between EP and ID, we first need to formally define EP in the context of the aforementioned computational graph defined in Def. [A..
Definition A.4 (Equilibrium Propagation (EP)). Denoting s the fixed point of Eq.
81
inside block l
and assuming that the transition kernel K has the form (s, W = V (( W x²),, we define
Equilibrium Propagation (EP) through the following recursive equations:
Vk -.... 1...1
8s Art k++ if N - else Vi((s,-11.1
Aay 0, 9 (0) 0, $3 5
Vt 1.... ;
9,1+1 = (50(5 W , rx = 1) 3ox*
8 W ", 14-11) Ary, - (5 Wk: xt 1 V32(s 1,14 Wk x 1)-1))
EP (3,1) = - 23
(2 (5 W x² = 8 ak-11 V2
11
Now that we have properly defined ID and EP, we are ready to state the main result of this section about the algorithmic equivalence between ID and EP which our coding work significantly built upon. Theorem A.3 (Extension of [Erooul et al. 2019 to ff-EBMs). Assuming that:
Vk =k=1.--- N 1, =t
=
(87)
where sh denotes the fixed-point of Eq. (81) and that the transition kernel K has the form (ss, W x*)) Vi Wh ,²*), then implicit differentiation (Def. A.3 and equilibrium prop- agation (Def. A.4 are equivalent in the limit /3 - O, namely:
vk =, N - =t 1..,
T
lim 8
F.P. 9
(3,t) = -9 (t)
ID
(88)
Proof of Theorem A.3 This proof follows the exact same methodology as that of Ernoult et al. [2019]
For self-containedness though and because of some subtle differences, carry out here the derivation.
We first define:
As; = des - dos
(89)
Note that since 2 -,,== 8x, d 1s-0 = 0 since does not depend on 0. Furthermore, note that by differentiating the equation satisfied by st.441 with respect to and evaluating the resulting expressions yields:
d 3, 1 1 1--0 O K(st, Wh = 14 dgss,rl8=0 55
(90)
In particular, evaluating Eq.
90) for = 0 yields:
As de d3st = - *
(91)
Therefore, substracting Eq.
90 across two timesteps yields altogether:
A 2 K(st, W rtk =k-1] A
VPT(stt W x. = =k--1 k-1) As
VP Wh at = k-k-T A
O Kss W t = A
(92)
31

---

Note that Oe (t) rewrites:
EP (3,t)
da (220( B,441 W , = - + (
- W II = 1
- W xt = A ( Wh ==== :dg + O(33)
2i Wh r² =
(3
22 wkt
82K( Wh 8 1TT (-Ai;) + + (3,1-1)+
(93)
Likewise, show that:
AZ B,t 0yK (s W x - 1T (-As;) + O(3))
(94)
Altogether, Eq.
Eq. 92 Eq. 933 and Eq. 94
yield, denoting J (t) lim 3-0 (8,1)
EP
and Ax = lim 3- Ar
vk N 1, 1:
- A $ = 8s
- ifk N else V V,E(s,--1 .)), IW (0) = 0. Art =0 (95)
11 ... T
A
Ks W , x = = (-As
I (t) A (sK Wk x st - (-Ast) + Ar
GEP 1) -g,((sh Wh x² sh sh-1) (-As;)
(96)
[1,] and N (8s;"1))eE(1) [[[] (Ax, (Ax and Starting from k 1, (- N-
)EE.. (WWN (t))ee[..1 and ( ~- (t) (t)) satisfy the same initial conditions and recursive equations, therefore there are all (pair-wise) equal for 1,..., Therefore in particular,
AxN-1 = No-N-1 such that (-As; N - ))E{1 and "os(">)ceeLL1) ")E[,.r from the previous (N - 2)th
block satisfy the same initial conditions, such that the reasonning applying to k N recurses for N 1. which yields Eq. 88).
32

---

A.4.3 Details about Fig. 3
Precise hyperparameters to reproduce Fig. 5 and Fig. 3 can be found inside our repository. Fig. precisely depict the architecture at use for these experiments.
Block
FastWon
Semdda
Block
Block 4
Comodds
Block 3
Cely
Block 2
Cemol28
Block
EastWorn
Figure 7: Architecture used for the static gradient analysis. The color code used to label layers matches that Fig. and Fig. 5 the context the static gradient analysis, "block" & is defined as all layers participating in A which therefore includes ph and E modules (rather than one of these taken alone).
33

---

A.5 Experimental Details
A.5.1 Datasets
Simulations were run on CIFAR-10, CIFAR-100 and Imagenet32 datasets, all consisting of color images of size 32 x 32 pixels. CIFAR-10 10 Krizhevsky 2009 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 (Krizhevsky] [209]] likewise comprises 60,000 and features 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 test set containing the remaining 10,000 images. The ImageNet32 dataset [Chrabasze7 et al., [2017) is downsampled version of the original ImageNet dataset Russakovsky et al. 2015 containing 1,000 classes with 1,281,167 training images, 50,000 validation images, 100,000 test images and 1000 classes.
A.5.2 Data preprocessing
All data were normalized according to statistics shown in and augmented with 50% random horizontal flips. Images were randomly cropped and padded with the last value along the edge of the image.
Table Data Normalization. Input images were normalized by conventional mean (p) and standard deviation (σ) values for each dataset. All images used are color (three channels).
Dataset
Mean (p)
Standard deviation (σ)
CIFAR-10/100
(0.4914, 0,4822, 0.4465) (0.485, 0,456, 0.406)
(0.2470, 0.2435, 0.2616)
Imagenet32
(0.3435, 0.336, 0.3375)
A.5.3 Simulation details
EP, similar to other machine learning paradigms reliant on fixed -oint iteration [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.
Weight initialization.
According to GOE, weights are intialized as:
[NO, K). ii [N N ifi=j
W
where NH,, o²) denotes Gaussian (normal) distribution with mean i and variance o². N was manually tuned for each architecture.
State initialization.
All layers 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 fixed- fixed-point procedure. the literature, activations (i.e. "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. [002]) depending on the loss function at use. For these experiments we used both the standard hard activation employed by Ernoult et al.) [019] and Scellier et al.) (2024), and the more conservative given in [Laborieux
2001 For the tolerance based and splitting experiments, we generalize the approach of Laborieux et al. 2021] by scaling values by a variable factor a instead
of fixed value 0.5. Details are given in Table
In practice, we find that the smaller scaling factors corresponding with the "laborieux" activation, in conjunction with GOE, and the omission of clamping at the output of each block significantly
34

---

Table Activation functions
Name
Description
Source
ernoult
(x))
max(min(r, 1),0)
Ernoult al. 2019
laborieux
(xx) = max(min(0.33 (0.5 x 1).0)
Laborieux et al. 2021
nest
0 (x) max(min(a 1).0)
This work
enhances gradient stability and speeds convergence in deep multi-block settings. In the interest of multi-scale uniformity and consistency with previous literature [Laborieux et al. 2021] Ernoult et al. (201)) apply clamping activations on all layers in our 6-layer architecture.
For the scaling experiments, we apply the "laborieux" activation 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 dynamics.
Cross-entropy loss and softmmax readout. Following [Lborriuxx 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: Wou of size dim (y) x dim(s), where is the state of the final layer. For each time-step 9 = softmax Doftmax(Wousss (oon))).
The cross- eoos--yntoppy cost function associated with the softmax readout is then:
l(s, y, Wou) Ye log(softmax. s)).
Convention to count layers. It is important to note that by convention we refer to architectures throughout this text the exclusion of the softmax readout, which technically an additional layer, though not involved in the relaxation process.
Architecture. All convolutional layers used experiments are of kernel size 3 and stride and padding 1. Max-pooling was applied with a window of 2 x 2 and stride of 2. For the 6-layer model used in Table [. batchnorm was applied after the first layer convolution and pooling operation. All other models in both experiments use batch-normalization convolution and pooling (where applied). These details exclude the linear softmax readout of size classes.
the first layer of each block after
Hyperparameters. Detailed hyperparameters for to reproduce Table I and Table 2 are given inside the configuration files of our repository. Note that all architectural details for the 12-layer models are identical across splitting and scaling experiments. Additionally, identical hyperparameters used for CIFAR100 and Imagenet experiments of Table 2 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 Table 2 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 algorithms. While in principle any root-finding algorithm may be used for the two relaxation phases of our EP implementation (for inference and gradient computation). our implementation utilized a simple fixed- exed-point iteration procedure, 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 Alg- As indicated in Section 4.3 two variants of the convergence procedure were employed, one in which the average value of current state 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 1c-4. This was known
35

---

as the tolerance- (TOL) procedure. Notably, tolerance- beasedce-based convergence criteria were applied on the free phase only, with nudging computed with a fixed value of iterations. This was to ensure consistency between ID) and EP, though in practice tolerance can be applied equally to the nudging phase.
Algorithm 13 Asynchronous with Tolerance (for all blocks until penultimate)
Input: T, 0%, wk, s4-1, B, 8s
Output: $B
&k
1: s +-0 0
2: c c 00
3: for tort=11 1... do do
4: V Wodd l E {1, Lb} ) B,ttmmp + (V & 1 $. 0h s5-2 wa)
5: Veven le (1,... [k}: &
Btttppp V T
(0 OK sh-1 :wo B55
6:
i 2 then
7:
VLE (1, ,,] ct +
8:
if mean (c*) S Tol then
9:
BREAK;
10: end if
11: end
12: $i. f- sk St,ttomp
13: end for
Supplementary results with a fixed number of iterations. In addition the TOL-based pro-
cedure, we obtained results for for using the more conventional approach of (Scellier and Bengio the first and second relaxation phases (see 1). This approach was also the default used for our scaling experiments
2019
Laborieux al. 2021 [Emoult al. 2019], applying fixed number of iterations
4.4
Importantly, with the TOL procedure described above Alg 3 becomes Alg 13] Results using
fixed iteration root-finding scheme are shown in [S]
Table Validation accuracy and Wall Clock Time (WCT) obtained on CIFAR-10 by EP (Alg. and ID on models with different number of layers (L) and block sizes ("bs"). 3 seeds are used.
EP
ID
Top-1 (%)
WCT
Top-1 1 (%)
WCT
L=
bs=6
88,8
8:06
87.3 10.6
8:05
bs=3
89.5
8:01
89.2
7:40
bs=2
90.1 1002
7:47
90.0
10.2
7:18
L=12
bs=4
91.6+0.11
7:49
91.6+0.1
7:08
bs=3
92.2
6:06
92.2
5:59
bs=2
91.7
6:10
91.8
6:08
Other details. All experiments were run using Adam optimizer (Kingma and Ba 2014 Cosine Annealing scheduler[Loshchilo and Hutter] 2017]) specifying some minimum learning rates and setting maximum T equal to epochs (i.e. warm restarts). Code was implemented in Pytorch 2.0 and simulations were run on NVIDIA A 100 SXM4 40GB GPUs.
36

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: we have a dedicated paragraph in the "Discussion" section of the paper which explicitly mentions limitations and future work.
3. 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: [Yes]
Justification: All proofs are included in the appendix.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: detailed information to reproduce experiments along with configuration files are provided on our github repository.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The URL to our github repository is provided in the main.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: all these details are provided in the appendix.
7. 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:
[Yes]
Justification: perform each of our training simulations on 3 different seeds and reported mean and standard deviation of the resulting performance.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: this information is also inside appendix.
9. Code Of Ethics
37

---

Question: Does the research conducted in the paper conform, in every respect. with the
NeurIPS Code of Ethics httigittere...accr....c///////it
Answer: [Yes]
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: we wrote dedicated paragraph inside our "Discussion" section.
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: Our work poses such risk present as it only provides proof-of-coneepts for systems which do not yet exist.
12. Licenses for existing assets
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?
Question:
Answer: [NA]
Justification: this work does not use existing assets.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: this paper does not release new assets.
14. 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: [NA]
Justification: paper does not involve crowdsourcing research with human subjects. 15. 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: [NA]
Justification: our paper does not involve crowdsourcing nor research with human subjects.
38


Reviews:
[{"review_id": "VEhJKvQPCG", "reviewer": "Reviewer_Jt5z", "summary": "Analog in-memory computing is gaining traction as an energy-efficient platform for deep learning. However, fully analog-based accelerators are challenging to construct, necessitating a training solution for digital-analog hybrid accelerators. This paper introduces Feedforward-tied Energy-based Models (ff-EBMs), a hybrid model that integrates feedforward components, typical in digital systems, with energy-based blocks suitable for analog circuits. An algorithm is derived to compute gradients end-to-end for training ff-EBMs. Experimental results show that ff-EBMs achieve superior accuracy on ImageNet32 classification compared to conventional equilibrium propagation literature.", "questions": "As I understand it, analog in-memory computing naturally finds the node voltage that minimizes the energy function through Kirchhoff's laws, classifying it as an EP problem. Therefore, if we disregard the natural minimization of the energy function, the training process with analog in-memory computing appears very similar to the conventional backpropagation-based training procedure. Hence, the derived algorithm for calculating end-to-end gradients of the ff-EBM seems trivial. In this context, I don't see the difference between this research and previous studies that have conducted training on analog in-memory systems using EP or backpropagation. Could you explain the differences between the previous approaches and the proposed approach in more detail? Providing examples of situations where ff-EBM is needed, as discussed in the Weakness section of this review, would be helpful.", "limitations": "Please check the weakness.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 1, "contribution": 3, "strengths": "1. This paper developed training algorithms specifically for digital-analog hybrid computing platforms, addressing a realistic computing scenario that leverages analog in-memory computing.", "weaknesses": "1. The paper lacks examples of situations where ff-EBM is needed. It does not adequately explain which parts of the overall system are digital and which are analog when using ff-EBM. Additionally, it does not clarify whether the existing training techniques can be used in these scenarios or what advantages ff-EBM offers compared to traditional methods.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer Jt5z for their honest feedback. We consider these clarifications important to highlight our contribution for readers.\n\ni) *Difference between this research and previous studies that have conducted training on analog in-memory systems using EP or backpropagation?*\n\nWe detail two research trends and contrast them with our work:\n\na) **Training analog systems using BP**. In this context, analog hardware sustains **feedforward models**, with two main techniques:\n\n- *On-chip BP training*: BP is directly executed on the analog hardware itself alongside inference. In general, forward and backward passes are necessarily executed on two **distinct** dedicated circuits. However, weights and activations must be *transported* from the inference circuit to the gradient computation circuit. Since weights/activations are encoded as noisy analog quantities (upon read and write operations), transporting them yields a *device mismatch* between these two circuits affecting in turn the resulting training performance [1].\n\n- *Off-chip BP training*: a digital proxy of the analog system is trained by BP on GPUs and the resulting weights are mapped onto the analog system, yielding in turn an analog *inference-only engine* [2]. This approach circumvents the difficulty of on-chip training of analog systems.\n\nb) **Training analog systems using EP**. When considering **energy-based** (EB) instead of feedforward models, the *same circuit* can be used for both the forward and backward passes. Therefore, when mapping an EB model onto a *single* analog circuit, the resulting system no longer suffers from device mismatch [1], with all the quantities required to compute gradients locally available at the same location, unlike BP. See global rebuttal and PDF for greater details.\n\nc) **This work**. While it is clear how to map models *small enough to fit single analog chip* and with standard weight stationary operations such as fully connected layers onto energy-based, analog compute engines trainable by EP [3, 4], models spanning *multiple analog chips* alongside many feedforward operations (e.g. batchnorm) to be maintained in digital with high precision result in systems which may not be readily trainable by EP – see Fig 1 of our attached PDF. **No existing work has tackled the extension of EP to this setting** where EP may only be applied within some subparts of the model, with BP applied elsewhere, in a principled fashion.\n\n**We accordingly propose to add an extra related work section related to this literature**.\n\nii) *Examples of situations where ff-EBMs are needed?*\n\nSee our global rebuttal and attached PDF. Additionally, the architecture considered for our experiments (L.213) is a very concrete example of this where it is unclear how to map batchnorm onto an energy-based, analog piece of hardware, and which may therefore be maintained in digital. Another example (mentioned L.304) would be transformer architectures – see Reviewer’s kDhg rebuttal and Fig. 2 of the PDF.\n\niii) *Which parts of the overall system are digital and which are analog when using ff-EBM?*\n\nAs explained in L.71-72, we “model digital and analog parts as feedforward and EB modules respectively”. This is further illustrated on Fig. 2, with “yellow and pink blocks denoting EB and feedforward transformations” and using the exact same notations, Eq. 12 explicitly defines the energy function of EB blocks and the feedforward transformations used for our experiments.\n\n**We propose to extend the color code used inside Fig. 2 to algorithms and equations to better emphasize which parts of the model and of the algorithm would be sustained in digital and analog**\n\niv) *Can existing training techniques can be used in these scenarios and what advantages ff-EBM offer compared to traditional methods?*\n\n- As explained in L. 229, ff-EBMs can be trained by traditional automatic differentiation (AD) through the root finding algorithm used to compute the equilibrium state of the ff-EBM (Eq. 13), which here reduces to Implicit Differentation (ID) – its implementation is detailed in Algorithm 12 below L. 605 inside the Appendix.\n\n- About the advantages of **EP-based training of ff-EBMs** compared to ID or AD: we want to emphasize that **AD is only used here for simulation benchmarking purposes** to demonstrate that our algorithm achieves the best possible performance and would not be practical to deploy onto analog hardware for all reasons mentioned above. Namely, exactly as for BP as described previously, ID would require two separated circuits to compute *the steady state of the ff-EBM* on the one hand, and *the associated Lagrangian multipliers* on the other hand.\n\n- Another advantage of **ff-EBMs alone**: they are *faster to simulate than standard EBM counterparts* (see Table 1 of the PDF and details inside Reviewer's).\n\n**We propose to add a more detailed mathematical definition of AD and ID when applied to ff-EBMs to better emphasize the benefits of EP-based training of ff-EBMs**\n\nv) *Disregarding energy minimization, the proposed algorithm appears very similar to BP and therefore seems trivial*\n\nWe assume that Reviewer Jt5z refers to the edge case where each EB block consists of a single layer such that the resulting ff-EBM is purely feedforward (“Recovering a feedforward net”, L.160), with our proposed algorithm reducing BP in this case (“Recovering backprop”, L.198, corollary A.5.1, Alg. 3 and Alg. 5 in appendix A.2). If so, we agree with Reviewer Jt5z **on this particular edge case**, and this is why **we discarded it upfront for our experiments** (L.201). However, in any other situation where EB blocks don’t reduce to a single layer, the proposed algorithm, deeply rooted in a Lagrangian-based approach thoroughly described in Appendix A.2, **doesn’t trivially reduce to any known algorithm in the literature**.\n\n[1] Yi et al, Nature Electronics, 2022\n[2] Wright et al, Nature, 2022\n[3] Kendall et al, ArXiV, 2020\n[4] Scellier, ICML 2024", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "xGvIVuSkXH", "reviewer": "Reviewer_k7Rs", "summary": "This paper presents Feedforward-tied Energy-based Models (ff-EBMs), a hybrid model that integrates feedforward and energy-based components, accounting for both digital and analog circuits. A novel algorithm is proposed to compute gradients end-to-end in ff-EBMs by backpropagating and \"eq-propagating\" through feedforward and energy-based sections, respectively, allowing EP to be applied to more flexible and realistic architectures. It has been shown that ff-EBMs can be trained on ImageNet32, achieving new state-of-the-art performance in the EP literature with a top-1 accuracy of 46%.", "questions": "See weaknesses.", "limitations": "NA", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "-- The proposed ff-EBMs as high-level models of mixed precision systems, where the inference pathway is composed of feedforward and EB modules, is interesting and novel. Specially, gradients computations as an end-to-end backpropagation through feedforward blocks and “eq-propagating” through EB blocks.\n\n-- The results are also encouraging specially on CIFAR datasets.\n\n-- The paper is easy to read and understand.", "weaknesses": "-- The primary limitation of this work is its accuracy performance on the ImageNet dataset. Although the paper has set a new state-of-the-art accuracy, it still lags behind the results achieved by Transformers and CNNs. There is a lack of compelling reasons to use this method given its comparatively lower performance.\n\n-- Additionally, it is crucial to measure the energy consumption of this training method and compare it with traditional methods. How much energy savings does your method offer?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer k7Rs for their comments and we are pleased they've appreciated our work.\n\ni) *Accuracy performance on the ImageNet dataset*\n\nAs  acknowledged inside the “Limitations and future work” paragraph (L.302-305), ff-EBM training by EP remains to be proved at scale on deeper models and more complex tasks (see Fig. 3 of the PDF about our research roadmap). Therefore, we  fully agree that the resulting performance of ff-EBMs **trained by EP** on ImageNet “lags behind” SOTA performance achieved by modern deep learning architectures, such as transformers or deeper CNN models, **trained by BP** on the same dataset.\n\nHowever, when comparing *different* models (ff-EBMs and transformers) trained by *different* training algorithms (EP and BP respectively), it is hard to disentangle improvements that come from the model *or* the training algorithm taken separately. Our results displayed inside Tables 1 and 2 of our original submission clearly show that **EP always performs as well as the automatic differentiation baseline on equal models**. This suggests (we believe) that the scalability of our approach is not limited by the proposed algorithm *itself*, but by the **ability to efficiently simulate large ff-EBMs**. To further support this, note that (as seen in Table 2)  training a ff-EBM of 15 layers takes at least 40 hours, *regardless of the training algorithm used*. These simulation times come from the lengthy fixed-point iteration used to compute equilibrium states inside EB blocks.\n\nTherefore, we envision two paths towards achieving better performance on ImageNet (see Reviewer kDhg's rebuttal and Fig. 3 of our PDF for further details):\n- using more GPU resources for longer times on deeper ff-EBMs, whilst preserving the core architectural components of the ff-EBM used in the present paper;\n- using EB blocks which converge faster to equilibrium. As mentioned in L. 297, using Deep Resistive Networks [1] instead of Deep Hopfield Networks as EB blocks, or even employing techniques such as Anderson Acceleration for root-finding, would yield a significant speed-up which would allow to train much deeper ff-EBMs.\n\n**We propose to make it clearer, in the “Limitations and future work” paragraph, that our simulations are bottlenecked by the ability to simulate large-scale ff-EBMs rather than the proposed algorithm itself, and will accordingly mention the two research paths mentioned above**\n\nii) *There is a lack of compelling reasons to use this method given its comparatively lower performance*\n\n- *Clarification*. We would like to clarify that the goal of this work is only to bring a **proof-of-concept** to guide the design of mixed-precision training accelerators with analog and digital parts. As ff-EBMs remain inefficient to simulate on GPUs, we do not advocate any practical use of ff-EBM over standard modern feedforward architectures, nor training them by EP instead of automatic differentiation, to achieve the best possible performance **when using GPUs**. For lack of widely accessible *real* mixed precision architectures to test our training algorithm against, ff-EBMs are high-level abstractions thereof which we *simulate* using GPUs.\n\n- ***ff-EBMs are much faster to simulate than their EBM counterparts***. With regards to **simulating** hardware systems as mentioned above, a compelling reason which did not appear clearly enough in our original submission is that **ff-EBMs are much faster to simulate than their single-block EBM counterpart of same depth**. To demonstrate this, we re-ran experiments of Table 1 (the \"splitting experiment\", section 4.3) with a *slightly different method to compute equilibria inside EB blocks*. We contrast these two settings below:\n  + *Former setting of the original submission* (Table 1 of the paper). Our initial goal was to demonstrate that starting from a standard (\"single-blocked\") EBM that ff-EBMs of *same depth* retained the same expressive power across all possible even splits, per the resulting performance on CIFAR-10. In this case, a fixed-point iteration scheme was applied inside each EB block, with the same *fixed* (and possibly large) number of iterations across all splits and throughout training. However, EB blocks of smaller size are expected to converge faster and therefore may require fewer fixed-point iterations [2-4].\n  + *New setting* (Table 1 of the PDF attached to the rebuttal). In order to take advantage of the faster convergence of smaller EB blocks and ensure the fairest comparison of wall-clock times across all splits *for a given depth* ($L=$6 or 12), we employ a *relative convergence criterion* to compute equilibria inside EB blocks. Namely, denoting $s(t)$ the state of a given EB block at the t-th fixed-point iteration, the fixed-point dynamics (Eq. 13 of our submission) are executed until $|(s(t + 1) - s(t))/ s(t)| < \\epsilon$ where $\\epsilon$ denotes some threshold. In spite of its simplicity, this trick was never employed in past EP works. Using this criterion, we can observe from Table 1 of the PDF attached to the rebuttal that ff-EBM training **can be up to 4 times faster that with equivalent EBM, while maintaining and often improving performance**.\n\nTo summarize: instead of the **superlinear** scaling of the convergence time of DHNs (with respect to the number of layers) empirically observed in past EP works [2--6], the convergence time ff-EBMs made up of DHNs are EB blocks are guaranteed, *by construction*, to scale **linearly** with the number of blocks, the convergence time of a single block decreasing with its size.\n\niii) *How much energy savings does your method offer*?\n\nWe kindly invite the Reviewer to read our global rebuttal for this question.\n\n[1] Scellier, ICML 2024\n[2] Ernoult et al, NeurIPS 2019\n[3] Laborieux et al, Frontiers in Neuro., 2021\n[4] Laborieux & Zenke, NeurIPS 2022\n[5] Scellier et al, NeurIPS 2023\n[6] Laborieux & Zenke, ICLR 2024", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "3xXS2E6YMe", "reviewer": "Reviewer_9JF8", "summary": "This paper proposes a new building model block with analog forward circuits and energy-based blocks built on a digital-analog hybrid setup.  A novel algorithm is further proposed to train the new block model. Experiments show the SOTA accuracy in the EP literature.", "questions": "N/A", "limitations": "N/A", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "* This paper achieves SOTA accuracy compared to recent literature.\n* A solid optimization method is proposed for the hybrid block.", "weaknesses": "* The paper's motivation is vague to me, as I am new to this topic. Why must we combine an analog forward circuit with a digital one for EP-based training? For the analog accelerator, we can use several other methods to do training, e.g., forward-forward only and zeroth-order optimization. What are the benefits of incorporating an energy-based model?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer 9JF8 for their comments and we are happy that they've appreciated our work.\n\ni) *Why must we combine an analog forward circuit with a digital one for EP-based training?*\n\nWe first want to clarify that per our modelling choice, analog and digital parts are respectively modeled as energy-based and feedforward models. Therefore, we do *not* assume that feedforward models are mapped onto analog. Modelling analog parts as energy-based models is a plausible assumption as it was shown that arrays of analog resistive devices minimize an energy function at equilibrium called the “pseudo-power” function [1, 2]. We kindly invite Reviewer 9JF8 to read details about these aspects inside the global rebuttal.\n\nAdditionally, EP-based training does not  *by itself* require the combination of analog (i.e. energy-based) and digital (i.e. feedforward) circuits. EP only applies, by default, to models which are *fully* energy-based, and that would map to a *fully* analog system. Rather, this combination is a realistic hardware design constraint, as analog hardware can currently only easily support vector-matrix multiplications and certain non-linear functions, with all other operations (e.g. max pooling, batchnorm) carried out in digital (see Fig. 1 of the PDF). Therefore, as such a system may not be *fully* energy-based but only *by parts*, it was unclear, before this work **how EP training would extend to such systems**.\n\nLastly, there is also an algorithmic advantage to use ff-EBMs as abstraction of such systems over EBMs: they are **faster to simulate**. Please read Reviewer k7Rs's rebuttal for further details, along with new results inside Table 1 of the PDF attached to our global rebuttal.\n\nii) *What are the benefits of incorporating an energy-based model, compared to using forward-forward only or zeroth-order optimization?*\n\nThis is an excellent question. Assuming Reviewer 9JF8 refers to the “Forward-Forward” (FF) algorithm [3] when mentioning “forward-forward only” optimization, both FF and zeroth-order (ZO) techniques apply to *any* model, be it energy-based or not, and estimate gradients through multiple forward passes through the *same* circuit, which is the holy grail for analog substrates. So why considering energy-based models at all?\n\nFF and ZO algorithms, while being mechanistically appealing for analog hardware [4], do not match the performance of automatic differentiation (i.e. backprop for feedforward models) on *equal* models, of around the same size as those studied in our work.\n\n- **ZO**. When using random weight perturbation [5] (WP, the best-known variant of ZO optimization, also known as “SPSA” [6]), resulting gradient estimates, although unbiased, have variance which scale cubically with the model dimensionality [7], yielding a significant gap in resulting model performance compared to backprop. This issue can be mitigated by architectural changes allowing the use of layer-wise, or even patch-wise losses [7,8], effectively reducing dimensionality and therefore the variance of ZO gradients or using small auxiliary networks to compute good “guesses” in the weight space (instead of random ones) [8]. These heuristics, while being effective at nearing backprop performance when training feedforward models from scratch, still don’t entirely close this performance gap on equal models. For instance, training a ResNet-18 by weight perturbation yields in the best case 45.8% top-1 test accuracy on ImageNet32 (using in this case auxiliary networks) while BP achieves 53.7% on the same architecture [8]. In contrast with ZO optimization, the use of EP on EB models or ff-EBMs **always match the performance of automatic differentiation on equal models**, in a principled fashion and without any heuristic, as observed in this work and past EP works [9, 10, 11]. Indeed, the energy-based requirement guarantees that EP gradients match automatic differentiation gradient in the limit of small nudging [12, 13].\nThis being said, we are aware that ZO techniques may however be sufficient to *fine-tune* pre-trained models [14], as pre-trained models somehow behave as small models (i.e. their Hessian has small rank). Therefore, neither first-order optimization techniques as EP, nor energy-based models, may be necessarily needed in the fine-tuning context, depending on the difficulty of the downstream task.\n\n- **FF**. The FF algorithm is endowed with the same high-level features of WP or ZO order algorithms, namely with the use of multiple forward passes with the same circuit and of local losses. However, it is a heuristic algorithm with no theoretical guarantees: the weight update does not estimate the gradient of the loss of interest, nor is it guaranteed to decrease the loss. Instead, the algorithm is designed to increase (resp. decrease) the “goodness” of each layer on positive (resp. negative) samples, with the goodness function and negative samples being heuristically motivated [3]. As a result, the resulting model accuracy trained by FF on CIFAR-10 does not exceed 55% on CIFAR-10 [15].\n\nIn summary, both ZO and FF algorithms, while possessing certain qualities and advantages similar to EP, are as yet not known to scale comparably to BP, whereas our algorithm (in principle) does.\n\n\n**We propose to add this detailed discussion around FF and ZO in the related work section**\n\n[1] *Johnson, 2010*\n[2] *Kendall et al, 2020*\n[3] *Hinton, 2022*\n[4] *Oguz et al, 2023*\n[5] *Fiete et al, 2010*\n[6] *Spall, 1992*\n[7] *Ren et al, 2022*\n[8] *Fournier et al, 2023*\n[9] *Laborieux et al, 2021*\n[10] *Laborieux & Zenke 2022*\n[11] *Scellier et al 2023*\n[12] *Ernoult et al 2019*\n[13] *Scellier & Bengio 2019*\n[14] *Malladi et al, 2023*\n[15] *Aminifar, 2024*", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "NAUN9OEG3V", "reviewer": "Reviewer_kDhg", "summary": "This paper discusses a novel approach to improving the power efficiency of AI training by integrating analog and digital hardware. The authors introduce a hybrid model called Feedforward-tied Energy-based Models (ff-EBMs) that combines feedforward neural networks with energy-based models (EBMs). This model aims to leverage the benefits of both analog and digital systems to create a more efficient computational framework for neural network training.\n\nThe paper introduces a new model, ff-EBMs, which integrates feedforward and EB modules. These models are designed to perform end-to-end gradient computation, combining backpropagation through feedforward blocks and \"eq-propagation\" through EB blocks. And bilevel and multi-level optimization are concepts used to frame the learning process of the Feedforward-tied Energy-based Models (ff-EBMs).\n\nThe paper promises to demonstrate the effectiveness of the proposed approach experimentally, using DHNs as EB blocks and achieving new state-of-the-art performance on the ImageNet32 dataset.", "questions": "In the proposed implicit BP-EP chaining algorithm, is it necessary to ensure the satisfaction of the first-order optimality conditions derived from the implicit theory, or can the algorithm be robust and effective without strictly guaranteeing these conditions?", "limitations": "A notable limitation of our study is the constrained scope of network architectures and datasets utilized for validation. While the implicit BP-EP chaining algorithm has shown success with Deep Hopfield Networks on ImageNet32, its performance on the standard ImageNet and transformer-based models, which are more complex and widely used, is yet to be established.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The paper presents a robust framework that leverages the strengths of bilevel and multi-level optimization to address the challenges of training neural networks in a hybrid digital-analog setting. One of the significant strengths lies in the rigorous mathematical proofs provided for the bilevel optimization inherent to energy-based models (EBMs). The authors extend this foundation to a multi-level optimization approach, which is not only theoretically sound but also practically applicable to the training of Feedforward-tied Energy-based Models (ff-EBMs).\n\n- The bilevel optimization is solidly established with the inner problem accurately reflecting the equilibrium state search, crucial for EBMs, while the outer problem adeptly encapsulates the parameter adjustment to minimize the cost function. This nested structure is then expanded into a multi-level optimization problem that inherently captures the complexity of training ff-EBMs, where each level represents an optimization challenge within the model's architecture.\n\n- The theoretical robustness is complemented by a well-thought-out algorithm that intertwines backpropagation and equilibrium propagation, offering a practical solution to the end-to-end training of these hybrid models. This algorithm is not just a theoretical construct; it is underpinned by a series of proofs that validate its effectiveness. The experimental results are particularly compelling, as they demonstrate the effectiveness of the proposed approach on ff-EBMs, achieving state-of-the-art performance on the ImageNet32 dataset.", "weaknesses": "While the paper introduces an innovative approach to hybrid digital-analog neural network training, there are a few areas where it shows some limitations. \n- The paper could benefit from clearer articulation regarding the energy efficiency claims, providing more detailed comparisons with existing systems to substantiate these claims. \n- Additionally, the reliance on simple datasets, ImageNet32, for experimental validation, while common, might limit the generalizability of the findings. A broader range of datasets with varying characteristics could strengthen the paper's conclusions. The choice of ImageNet32 should be justified with respect to its relevance to the research goals and its limitations for testing the model's robustness.\n- Lastly, while the authors acknowledge the need for further research, a more explicit discussion on the current limitations and a roadmap for future work would provide a more comprehensive view of the research's trajectory.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer kDhg for their detailed review and we are happy that they've appreciated our work.\n\ni) *Energy efficiency claims, detailed comparisons with existing systems to substantiate these claims*. Please see our global rebuttal.\n\nii) *Why only ImageNet32?*\n\nImageNet32 was used for two reasons:\na) we wanted to align our work with the most recent EP literature evaluating algorithms on vision tasks only with dimensionality at most 32 [1, 2, 3];\nb) The reason we restricted ourselves to such a low-dimensional dataset is because ff-EBM training simulations are already *relatively slow* on simple datasets – see below for a more detailed discussion.\n\niii) *Explicit discussion on the current limitations and a roadmap for future work*.\n\nWe herein detail the roadmap for future work which addresses three major limitations of our  paper: 1) Scalability to complex tasks; 2) Generalizability to other datasets and architectures; 3)  Hardware realistic training of ff-EBMs.\n\n**Please find the associated Fig. 3 inside the PDF attached to our global rebuttal**.\n\n1) **Scalability to complex tasks**, i.e. *higher input dimensionality* (e.g. going from ImageNet32 to ImageNet224) and *deeper and larger architectures*. As we wrote to Reviewer k7RS, we insist that what bottlenecks the scalability of our work is not the proposed algorithm itself but the **ability to efficiently simulate ff-EBM inference**, because of the lengthy root finding algorithms needed to compute equilibria inside the EB blocks. As highlighted by Table 2, training a ff-EBM of 15 layers on ImageNet32 takes at least 40 hours, *regardless of the training algorithm used*. Therefore, “scaling” either requires:\n - a) *maintaining the current architecture*, i.e. with Deep Hopfield Networks (DHNs) as EB blocks, on ImageNet224, *with more GPUs and more compute time*, or:\n - b) *changing the EB block*, i.e. use EB blocks with better convergence properties than Deep Hopfield Networks (DHNs). A compelling class of models is called Deep Resistive Nets (DRNs) [4] which have the double advantage to be more hardware plausible and faster to simulate owing to their convex energy function.\n\n2) **Generalizability to other datasets and architectures**. As pointed out by Reviewer 9JF8, our approach has yet to be proved on *transformer architectures* and associated datasets, which subsumes:\n\n- a) *Defining a ff-EBM counterpart of standard transformers*. To map a transformer into a ff-EBM architecture, we need to break it up into feedforward / digital and energy-based / analog blocks. *As illustrated in Fig. 2 inside the PDF attached to our rebuttal*, the ff-EBM counterpart of a transformer encoder block would comprise: i) a feedforward block with normalization and attention layers, ii) an energy-based block with the two fully connected layers. We denote below this ff-EBM counterpart of a transformer architecture “TF-ff-EBM”.\n\n- b) *Testing a TF-ff-EBM with DHN blocks*. As we did in this study, we would use *DHNs as EB blocks* in the proposed TF-ff-EBM architecture and test it on a relatively simple dataset to start with (e.g. the PTB dataset).\n\n- c) *Testing a TF-ff-EBM with DRN blocks*. See 3a) below.\n\n- d) *Scaling hardware-realistic training of TF-ff-EBM to more complex tasks*. We would scale the previous approach to deeper TF-ff-EBMs on the Wikitext datasets, with several additions to the architecture depicted in c) to make it more hardware realistic (see 3c) below).\n\n3) **Hardware realistic training of ff-EBMs**. Moving towards more realistic proof-of-concepts of mixed-precision training of ff-EBMs, the following milestones need to be hit applying to both vision models considered in this paper as well as language models to be investigated:\n\n- a) *Using DRNs as EB blocks*.  As mentioned earlier, DRNs are *exact* models of nonlinear resistive analog circuits, made up of resistors, diodes, operational amplifiers with actual voltage nodes as neurons.\n\n- b) *Using quantized activations and gradients within feedforward blocks*. A realistic assumption for feedforward transformation would be to employ ***quantized** activations, weights and **gradients***. Techniques such as Quantization-aware Scaling (QAS) could be employed to compute meaningful quantized gradients inside feedforward blocks [5].\n\n- c) *Taking into account analog non-idealities*. As pointed out by Reviewer pevw, we need to consider a “detailed analog computing device model”, for instance taking “device variation” into account or any other device imperfection which may affect gradient computation. Taking these into account is crucial to assess how much extra digital circuitry it would take to mitigate these as it entails energy consumption.\n\n- d) *A comprehensive energy consumption analysis of hardware realistic ff-EBM training*. Finally, *quantitatively* and *precisely* assessing the energy efficiency gains of ff-EBM training with the more hardware-realistic assumptions introduced above is crucial.\n\n**We propose to add this detailed research roadmap inside the supplementary material of the camera-ready version of our paper in case of acceptance**\n\niv) *Should the first-order optimality conditions be strictly guaranteed for the BP-EP chaining to work?*\n\nIn practice, meeting optimality conditions can be quantitatively assessed with the static gradient analysis presented inside Fig. 3 and Fig. 4 (since $\\lambda_\\star$ can be typically be computed by automatic differentiation) and we observe that **the “more optimal” $s_{\\star}^{\\pm \\beta}$ are** (in the sense of satisfying the KKT conditions) **,  the better the resulting training performance by EP**. For a further mathematical analysis of approximate first-order optimality in implicit models, we refer to [6, 7].\n\n[1] *Laborieux & Zenke, NeurIPS 2022*\n[2] *Laborieux & Zenke, ICLR 2024*\n[3] *Scellier et al, NeurIPS 2023*\n[4] *Scellier, ICML 2024*\n[5] *Lin et al, NeurIPS 2022*\n[6] *Zucchet & Sacramento, 2022*\n[7] *Meulemans et al, NeurIPS 2022*", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "SOETgUyKIm", "reviewer": "Reviewer_pevw", "summary": "State-of-the-art (SOTA) analog hardware accelerators consist of both analog and digital components supporting major and auxiliary operations. Moreover, they typically suffer from device imperfections on their analog parts. In this paper, the authors propose feedforward-tied energy-based models (ff-EBMs) for digital and analog circuits. ff-EBMs compute gradients by both backpropagating on feedforward parts and eq-propagating on energy-based parts. In this paper, ff-EBMs use a Deep Hopfield Network (DHN) as one example, which can be arbitrarily partitioned into uniform size. ff-EBMs achieves a SOTA performance on ImageNet32.", "questions": "Please comment on the points in the weakness section.", "limitations": "no limitations.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper works on an important and interesting problem.\n2. The paper flows well.", "weaknesses": "1. The paper does NOT consider a detailed analog computing device model. It will be interesting to see how significant variations on the analog devices the ff-EBMs can tolerate. And what is the relationship between the convergence speed and the device variations?\n\n2. The paper does not include energy-related data. What is the energy saving for the energy-based models? What is the energy saving of skipping the accurate gradient computations?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer pevw for their detailed review and we are happy that they've appreciated our work.\n\ni) *The paper does NOT consider a detailed analog computing device model*.\n\nWe fully agree with Reviewer pevw that this is a current limitation of our work, as explicitly acknowledged inside the “Limitations and future work” paragraph (L.294-301). This is why moving towards **hardware realistic ff-EBM training is one of the main thrusts of our research trajectory**  - see section 3 of the proposed roadmap to Reviewer kDhg as well as Fig. 3 inside the pdf attached with the global rebuttal).\n\nThis being rightfully pointed out, there are still good reasons not to start the investigation of ff-EBM training by EP with realistic analog computing device model but instead Deep Hopfield Networks (DHNs). DHNs, while being remotely inspired by physics [1] and without “analog-realistic” noise, serve as a useful model baseline to start with in the context of EP for at least three reasons:\n\n- In all past studies investigating the scalability of EP training and motivated by analog computing, DHNs were used as EB models [2, 3, 4, 5].\n\n- DHNs are *already difficult to train by EP*. Compared to feedforward models trained by backprop, training DHNs by EP comes with extra critical and subtle hyperparameter tuning, most importantly: a) the root find algorithm at use to compute equilibria [6, 7, 8]; b) the number of iterations needed to converge in the first, second and third phases of EP; c) the value of the nudging parameter $\\beta$; d) the way nudging is applied in the output layer [9]; e) weight initialization statistics, which implicit learning more generally is highly sensitive to [6, 7]. For instance, our study reveals that initializing the weights of the EB blocks as *Gaussian Orthogonal Ensembles* (abbreviated as “GOE weight initialization” in our paper) was instrumental in having ff-EBM training by EP working properly.\n\n- Finally, simulating energy-based analog resistive networks comes with even more hyperparameters compared to DHNs [8, 10], making in turn EP training harder to debug.\n\nii)\t*How much can ff-EBMs tolerate significant variations on the analog devices?*\n\nWe would like to clarify the scope of the question by defining two broad classes of analog device imperfections which Reviewer pevw may have in mind:\n\n- *Device imperfections affecting the **weight update** given a gradient value*. In analog systems, weights are typically encoded as conductance values of resistive devices [10]. *Given a gradient estimate*, programming a resistive device to update its conductance by this gradient estimate *precisely* is generally very hard, if not impossible because conductance updates suffer from device-to-device (“static”) **variability**, cycle-to-cycle (“dynamic”) **variability**, conductance drift, non-linearity with respect to the applied programming voltage, asymmetry between device potentiation and depression, etc. [11]. However, these imperfections **apply to any gradient computation algorithm** since they only pertain to the conductance update *given a gradient estimate*. As highlighted in the most recent study of “analog SGD” [14], the gradient estimate is *given*. Therefore, all conclusions regarding “analog parameter optimization” *given a gradient value* hold in particular for ff-EBM training by EP.\n\n- *Device imperfections affecting the **gradient estimation***. Once an analog nonlinear resistive circuit is instantiated with some *fixed* conductance values, it may still be subject to *static variability* (i.e. two analog devices do not share the exact same input-to-output characteristic) and *dynamic variability* (i.e. a single analog device may not react exactly the same when exposed twice to the same input). However, as highlighted in the global rebuttal and that addressed to Reviewer Jt5z, **EP-based training of EB blocks is inherently tolerant to this static variability** since this very same circuit is used for *both* inference *and* gradient computation [12]. Therefore, **ff-EBM training by EP may only be impacted by dynamic variability**. Investigating quantitatively the tolerance of ff-EBM training to dynamic variability and coming up with strategies to mitigate its potential negative effects is a project of its own that naturally falls in our research agenda (see task 3c) inside rebuttal to Reviewer kDhg and in Fig. 3 inside the pdf attached to the global rebuttal.\n\niii)\t*What is the relationship between the convergence speed and the device variations?*\nWe assume that  Reviewer pevw has in mind the convergence of *the learning dynamics of ff-EBMs **in the parameter space*** when trained by the proposed algorithm. The “device variations” may refer in this case to those involved in the *conductance updates* which, as mentioned above, are *agnostic to how gradients are computed*. It is known that the learning dynamics of \"analog SGD\" [14] are severely impacted by such non idealities even on simple tasks as MNIST [13]. Therefore, any such observation, or theoretical prediction [14] as well as any heuristic to mitigate non idealities of the conductance updates, e.g. Tiki-Taka [15], *are also directly applicable to ff-EBM training by EP*.\n\niv)\t*What is the energy saving for the energy-based models? What is the energy saving of skipping the accurate gradient computations?*\nPlease, see our global rebuttal.\n\n\n[1] *Hopfield, Science, 1986*\n[2] *Laborieux et al, Frontiers in Neuroscience, 2021*\n[3] *Scellier et al, 2022*\n[4] *Laborieux & Zenke, NeurIPS 2022*\n[5] *Laborieux & Zenke, ICLR 2024*\n[6] *Bai et al, NeurIPS 2019*\n[7] *Agarwala and Schoenholz, ICML 2022*\n[8] *Scellier, ICML 2024*\n[9] *Scellier et al, NeurIPS 2023*\n[10] *Kendall et al, 2020*\n[11] *Burr et al, Advance in Physics X, 2017*\n[12] *Yi et al, Nature Electronics, 2022*\n[13] *Burr et al, IEEE Transactions on Electron Devices, 2015*\n[14] *Wu et al, ArXiV, 2024*\n[15] *Gokmen & Haensch, Frontiers in Neurosciences, 2020*", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank reviewers for their time and highly valuable comments. In light of these, we propose several clarifications which we hope render the value of our work more clear to readers and can quell some of the concerns expressed.

**I-Proposed additions to our paper**
- A  discussion about the **relevance of ff-EBMs** in analog computing and **energy efficiency gains of the proposed method** (see *inside this global rebuttal*), including a sketch of a hybrid chip with analog and digital parts with EP training at chip scale (Fig. 1 inside PDF).
- A **detailed roadmap for future research** (Reviewer kDhg’s rebuttal and Fig. 3 inside the PDF). This includes a  discussion about **hardware realism of ff-EBM training** (Reviewer pevw’s rebuttal).
- Explanation of **how transformer architectures can be mapped onto ff-EBM architectures** and hence trained by the proposed algorithm (Reviewer kDhg’s rebuttal and Fig. 2 of PDF).
- A  related work section about **training analog in-memory systems using EP or BP** (Reviewer Jt5z’s rebuttal).
- A related work section with **detailed comparison of energy-based learning with 0-th order optimization (ZO) and the forward-forward algorithm (FF)** (Reviewer 9JF8’s rebuttal).
- Detailed mathematical definition of automatic differentiation (AD) and implicit differentiation (ID) when applied to ff-EBMs **to better emphasize the benefits of EP-based training of ff-EBMs** (Reviewer Jt5z’s rebuttal).

**II- New experiments**

To further justify why ff-EBMs are useful in practice (Reviewers Jt5z, 9JF8), we highlight in Table 1 of the PDF that **they are much faster to simulate than EBM counterparts of equal depth** (see details below and inside Reviewer k7Rs's rebuttal).

**III-The relevance of ff-EBM modelling**

For some readers, “the paper’s motivation [can look] vague” (Reviewer 9JF8) and it may be unclear “when ff-EBMs [are] needed” (Reviewer Jt5z). Therefore, we clarify some fundamentals underpinning our work, alongside new results (Fig. 1 and Table 1 inside the PDF) of benefit to all Reviewers. Importantly, our answer below relies on an **end-to-end experimental realization of EP training on analog resistive networks** [6].

i) *Why analog computing?* Analog in-memory computing (AIMC) accelerators achieves greater energy efficiency by: a) encoding weights as conductance values of resistive elements in an analog crossbar array *which can be read and written at very low cost* [5]; b) leveraging analog physics, i.e. continuous physical observables, to perform *vector-matrix multiplications (VMMs) at very low cost* [1 – 5].

ii) *Why modelling AIMC systems as energy-based models?* Nonlinear resistive networks as AIMC systems are *energy-based* models [1, 4], as Kirchhoff laws, which govern these systems, obey a variational (**energy minimization**) principle. Energy-based models are convenient as they can compute loss gradients using multiple relaxations to equilibrium using a *single circuit*, as prescribed by EP. This is why, in line with these works, we also model AIMC cores as energy-based models.

iii) *Why are AIMC systems not sufficient alone?* As highlighted in the scalability study of an experimental realization of EP [6], **AIMC energy-based cores may only constitute the smallest compute unit of a larger system comprising many digital parts**. We adapt from [6] a simplified outline of an “architecture to implement [EP training] on a chip scale” in *Fig. 1 of the attached PDF*. Quoting from the authors of [6], “[this] architecture consists of **hybrid (analogue-digital) operations** [and is] hierarchical”: as seen from our Fig. 1, several analog processors form a tile, and several tiles form the chip, with digital buses connecting analog cores along with digital coprocessors to support analytical operations (e.g. maxpool, as in [6]). Therefore a **ff-EBM**, as inherently hierarchical and hybrid (Eq. 8 of our paper), **is a coarse-grained abstraction of such a modular architecture**, which EP is extended to in our work.

iv) **ff-EBMs are faster to simulate**. We reproduced the splitting experiments with a *uniform convergence criterion* to compute equilibria inside EB blocks across all splits (instead of tuning the number of iterations for each of these). **See Reviewer k7Rs's rebuttal for details and Table 1 of the PDF**.

**IV-Energy efficiency of ff-EBM training**

Reviewers k7RS, kDhg and pevw asked about the “energy consumption of [our] training method” and “detailed comparisons with existing systems”. As mentioned earlier, **energy savings fundamentally would come from lesser costs for VMMs and read/write operations on resistive devices inside EB blocks**. Providing an accurate answer to this question requires research beyond the scope of this paper (see Fig. 3 for a detailed outline of our research roadmap) and is highly system dependent. Based on [6], it is possible to convey this complexity and provide an estimation of the energy gains *compared to a NVIDIA V100 device* (for inference + gradient computation + weight update) depending on a variety of  factors. Most importantly:

- *System’s size*: if the model fits a single AIMC core (a 64x64 array of RRAM devices), energy gains are around $10^5$. If instead considering a model spanning multiple cores within the hybrid architecture previously described from [6], **which ff-EBMs most resemble**, energy gains reduce to $10^4$ because of the “architectural overheads, dominated by the peripheral circuitry” and “analog-to-digital conversions” [6]
- *Resistive devices*: the above numbers assume RRAM resistive devices. Considering Flash memory, writing is around 100x more energy consuming than on RRAM, yielding lesser energy savings if batch size is small. These numbers also change for transistor-based [2, 3] or magnetic-based [6] synapses.

[1] *Kendall et al, 2020*
[2] *Dillavou et al, 2022*
[3] *Dillavou et al, 2023*
[4] *Scellier, 2024*
[5] *Burr et al, 2017*
[6] *Yi et al, Nature Electronics, 2022*

Meta-review: The authors present feedforward tied Energy Based Models (ff-EBMs) a new class of models that combine traditional feedforward layers with EBM layers. The paper also describes an end-to-end training methodology for ff-EBMs that chains the gradient estimations via Backprop (for ff layers) and Eq Prop (for EBM layers). This training methodology is supported by a detailed proof that confirms its correctness and shows that the ff-EBM training methodology works for any arbitrary split of a network into ff/EBM layers.

In the experiments section, the authors apply ff-EBM models for CIFAR10/CIFAR100/ImageNet32 datasets and show competitive accuracy. They also validate the equivalence of the gradients estimated using the proposed training methodology against traditional backprop based gradients and the validity of their training methodology for a range of ff/EBM splits for target model. 

While the paper presents a strong theoretical explanation of ff-EBMs training and inference, many reviewers raised interesting questions about why ff-EBMs should be of interest and what are the benefits of Energy based models. The Authors have provided a detailed response to every question and provided useful figures and explanations of ff-EBMs usefulness and the HW system they have in mind to deploy these models. 

A couple of points that stood out to me in terms of this work’s importance are as follows. Digital AIML HW implementations show high accuracy but suffer from high power cost, Analog HW implementation promises to counter some of this high-power cost. However, the analog HW suffers from reliability concerns and slow weight updates, they aren’t suitable for backprop based training and non-weight stationery operations (e.g. Batch Normalization). While EBMs offer a solution that enables reliable training on analog HW, it still depends on Digital HW for implementing non-weight stationery operations. The proposed ff-EBMs are a pioneering solution to combine the best of both worlds, by using analog HW wherever it enables low energy operations and outsourcing the rest to digital HW. Thanks to the training methodology demonstrated in this paper, any given model can be arbitrarily split onto analog/digital HW. 

Further, EBMs currently suffer from long simulation times due to the delays involved in the root-finding step for each layer of EBMs. By combining ff-EBMs the authors speed up the training steps significantly (~4x as shown in Table 1 of rebuttal pdf). This in turn enables deploying networks with deeper EBM layers interspersed with ff layers to achieve superior accuracy compared to existing EBMs. 

The authors also add a detailed timeline (Figure 3 rebuttal pdf) of their future work indicating extension to other datasets (ImageNet244) and model architectures (Transformer). In conclusion, ff-EBMs and their proposed training algorithm is a pioneering step in the direction of hybrid (digital/analog) HW systems implementing efficient AIML algorithms. This paper should be accepted for NeurIPS.

Edit: Updated my recommendation to spotlight (prev Oral) post discussion with SAC.
Decision: Accept (spotlight)