CONTINUAL WEIGHT UPDATES
Year: 2019

Abstract:
Equilibrium Propagation (EP) is a learning algorithm that bridges Machine Learning and Neuroscience, by computing gradients closely matching those of Backpropagation Through Time (BPTT), but with a learning rule local in space. Given an input x and associated target y, EP proceeds in two phases: in the first phase neurons evolve freely towards a first steady state; in the second phase output neurons are nudged towards y until they reach a second steady state. However, in existing implementations of EP, the learning rule is not local in time: the weight update is performed after the dynamics of the second phase have converged and requires information of the first phase that is no longer available physically. In this work, we propose a version of EP named Continual Equilibrium Propagation (C-EP) where neuron and synapse dynamics occur simultaneously throughout the second phase, so that the weight update becomes local in time. Such a learning rule local both in space and time opens the possibility of an extremely energy efficient hardware implementation of EP. We prove theoretically that, provided the learning rates are sufficiently small, at each time step of the second phase the dynamics of neurons and synapses follow the gradients of the loss given by BPTT (Theorem 1). We demonstrate training with C-EP on MNIST and generalize C-EP to neural networks where neurons are connected by asymmetric connections. We show through experiments that the more the network updates follows the gradients of BPTT, the best it performs in terms of training. These results bring EP a step closer to biology by better complying with hardware constraints while maintaining its intimate link with backpropagation.

heading: 1 INTRODUCTION
text: A motivation for deep learning is that a few simple principles may explain animal intelligence and allow us to build intelligent machines, and learning paradigms must be at the heart of such principles, creating a synergy between neuroscience and Artificial Intelligence (AI) research. In the deep learning approach to AI (LeCun et al., 2015), backpropagation thrives as the most powerful algorithm for training artificial neural networks. Unfortunately, its implementation on conventional computer or dedicated hardware consumes more energy than the brain by several orders of magnitude (Strubell et al., 2019). One path towards reducing the gap between brains and machines in terms of power consumption is by investigating alternative learning paradigms relying on locally available information, which would allow radically different hardware implementations: such local learning rules could be used for the development of extremely energy efficient learning-capable hardware. Investigating such bioplausible learning schemes with real-world applicability is therefore of interest not only for neuroscience, but also for developing neuromorphic computing hardware that takes inspiration from information-encoding, dynamics and topology of the brain to reach fast and energy efficient AI (Ambrogio et al., 2018; Romera et al., 2018). In these regards, Equilibrium Propagation (EP) is an alternative style of computation for estimating error gradients that presents significant advantages (Scellier and Bengio, 2017).
EP belongs to the family of contrastive Hebbian learning (CHL) algorithms (Ackley et al., 1985; Movellan, 1991; Hinton, 2002) and therefore benefits from an important feature of these algorithms: neural dynamics and synaptic updates depend solely on information that is locally available. As a CHL algorithm, EP applies to convergent RNNs, i.e. RNNs that are fed by a static input and converge to a steady state. Training such a convergent RNN consists in adjusting the weights so that the
steady state corresponding to an input x produces output values close to associated targets y. CHL algorithms proceed in two phases: in the first phase, neurons evolve freely without external influence and settle to a (first) steady state; in the second phase, the values of output neurons are influenced by the target y and the neurons settle to a second steady state. CHL weight updates consist in a Hebbian rule strengthening the connections between co-activated neurons at the first steady state, and an anti-Hebbian rule with opposite effect at the second steady state. A difference between Equilibrium Propagation and standard CHL algorithms is that output neurons are not clamped in the second phase but elastically pulled towards the target y.
A second key property of EP is that, unlike CHL and other related algorithms, it is intimately linked to backpropagation. It has been shown that synaptic updates in EP follow gradients of recurrent backpropagation (RBP) (Scellier and Bengio, 2019) and backpropagation through time (BPTT) (Ernoult et al., 2019). This makes it especially attractive to bridge the gap between neural networks developed by neuroscientists, neuromorphic researchers and deep learning researchers.
Nevertheless, the bioplausibility of EP still undergoes two major limitations. First, although EP is local in space, it is non-local in time. In all existing implementations of EP the weight update is performed after the dynamics of the second phase have converged, when the first steady state is no longer physically available. Thus the first steady state has to be artificially stored. Second, the network dynamics have to derive from a primitive function, which is equivalent to the requirement of symmetric weights in the Hopfield model. These two requirements are biologically unrealistic and also hinder the development of efficient EP computing hardware.
In this work, we propose an alternative implementation of EP (called C-EP) which features temporal locality, by enabling synaptic dynamics to occur throughout the second phase, simultaneously with neural dynamics. We then address the second issue by adapting C-EP to systems having asymmetric synaptic connections, taking inspiration from Scellier et al. (2018) ; we call this modified version C-VF.
More specifically, the contributions of the current paper are the following:
• We introduce Continual Equilibrium Propagation (C-EP, Section 3.1-3.2), a new version of EP with continual weight updates: the weights of the network are adjusted continually in the second phase of training using local information in space and time. Neuron steady states do not need to be stored after the first phase, in contrast with standard EP where a global weight update is performed at the end of the second phase. Like standard EP, the C-EP algorithm applies to networks whose synaptic connections between neurons are assumed to be symmetric and tied. • We show mathematically that, provided that the changes in synaptic strengths are sufficiently
slow (i.e. the learning rates are sufficiently small), at each time step of the second phase the dynamics of neurons and synapses follow the gradients of the loss obtained with BPTT (Theorem 1 and Fig. 2, Section 3.3). We call this property the Gradient Descending Dynamics (GDD) property, for consistency with the terminology used in Ernoult et al. (2019). • We demonstrate training with C-EP on MNIST, with accuracy approaching the one obtained
with standard EP (Section 4.2). • Finally, we adapt our C-EP algorithm to the more bio-realistic situation of a neural network
with asymmetric connections between neurons. We call this modified version C-VF as it is inspired by the Vector Field method proposed in Scellier et al. (2018). We demonstrate this approach on MNIST, and show numerically that the training performance is correlated with the satisfaction of Gradient Descending Dynamics (Section 4.3).
For completeness, we also show how the Recurrent Backpropagation (RBP) algorithm of Almeida (1987); Pineda (1987) relates to C-EP, EP and BPTT. We illustrate the equivalence of these four algorithms on a simple analytical model (Fig. 3) and we develop their relationship in Appendix A.

heading: 2 BACKGROUND: CONVERGENT RNNS AND EQUILIBRIUM PROPAGATION
text: Convergent RNNs With Static Input. We consider the supervised setting, where we want to predict a target y given an input x. The model is a recurrent neural network (RNN) parametrized by θ
and evolving according to the dynamics:
st+1 = F (x, st, θ) . (1)
F is the transition function of the system. Assuming convergence of the dynamics before time step T , we have sT = s∗ where s∗ is the steady state of the network characterized by
s∗ = F (x, s∗, θ) . (2)
The number of timesteps T is a hyperparameter that we choose large enough so that sT = s∗ for the current value of θ. The goal is to optimize the parameter θ in order to minimize a loss:
L∗ = ` (s∗, y) . (3)
Algorithms that optimize the loss L∗ for RNNs include Backpropagation Through Time (BPTT) and the Recurrent Backpropagation (RBP) algorithm of Almeida (1987); Pineda (1987), presented in Appendix B.
Equilibrium Propagation (EP). EP (Scellier and Bengio, 2017) is a learning algorithm that computes the gradient of L∗ in the particular case where the transition function F derives from a scalar function Φ, i.e. with F of the form F (x, s, θ) = ∂Φ∂s (x, s, θ). The algorithm consists in two phases (see Alg. 1 of Fig. 1). During the first phase, the network follows a sequence of states s1, s2, s3 . . . and converges to a steady state denoted s∗. In the second phase, an extra term β ∂`∂s pertubs the dynamics of the neurons (where β > 0 is a scalar hyperparameter): starting from the steady state sβ0 = s∗, the network follows a second sequence of states s β 1 , s β 2 , s β 3 . . . and converges to a new steady state denoted sβ∗ . Scellier and Bengio (2017) have shown that the gradient of the loss L∗ can be estimated based on the two steady states s∗ and sβ∗ . Specifically, in the limit β → 0,
1
β
( ∂Φ
∂θ
( x, sβ∗ , θ ) − ∂Φ ∂θ (x, s∗, θ) ) → −∂L ∗ ∂θ . (4)

heading: 3 EQUILIBRIUM PROPAGATION WITH CONTINUAL WEIGHT UPDATES (C-EP)
text: This section presents the main theoretical contributions of this paper. We introduce a new algorithm to optimize L∗ (Eq. 3): a new version of EP with continual parameter updates that we call C-EP. Unlike typical machine learning algorithms (such as BPTT, RBP and EP) in which the weight updates occur after all the other computations in the system are performed, our algorithm offers a mechanism in which the weights are updated continuously as the states of the neurons change.

heading: 3.1 FROM EP TO C-EP: AN INTUITION BEHIND CONTINUAL WEIGHT UPDATES
text: The key idea to understand how to go from EP to C-EP is that the gradient of EP appearing in Eq. (4) reads as the following telescopic sum:
1
β
( ∂Φ
∂θ
( x, sβ∗ , θ ) − ∂Φ ∂θ (x, s∗, θ) ) ︸ ︷︷ ︸
global parameter gradient in EP
= ∞∑ t=1 1 β ( ∂Φ ∂θ ( x, sβt , θ ) − ∂Φ ∂θ ( x, sβt−1, θ )) ︸ ︷︷ ︸
parameter gradient at time t in C-EP
. (5)
In Eq. (5) we have used that sβ0 = s∗ and s β t → s β ∗ as t→∞. Here lies the very intuition of continual updates motivating this work; instead of keeping the weights fixed throughout the second phase and updating them at the end of the second phase based on the steady states s∗ and s β ∗ , as in EP (Alg. 1 of Fig. 1), the idea of the C-EP algorithm is to update the weights at each time t of the second phase between two consecutive states sβt−1 and s β t (Alg. 2 of Fig. 1). One key difference in C-EP compared to EP though, is that, in the second phase, the weight update at time step t influences the neural states at time step t+ 1 in a nontrivial way, as illustrated in the computational graph of Fig. 2. In the next subsection we define C-EP using notations that explicitly show this dependency.
Algorithm 1 EP Input: x, y, θ, β, η. Output: θ.
1: s0 ← 0 . First Phase 2: repeat 3: st+1 ← ∂Φ∂s (x, st, θ) 4: until st = s∗ 5: Store s∗ 6: sβ0 ← s∗ . Second Phase 7: repeat 8: sβt+1 ← ∂Φ∂s ( x, sβt , θ ) − β ∂`∂s ( sβt , y
) 9: until sβt = s β ∗
10: . Global Parameter Update 11: θ ← θ + ηβ ( ∂Φ ∂θ ( sβ∗ , θ ) − ∂Φ∂θ (s∗, θ) )
Algorithm 2 C-EP (with simplified notations) Input: x, y, θ, β, η. Output: θ.
1: s0 ← 0 . First Phase 2: repeat 3: st+1 ← ∂Φ∂s (x, st, θ) 4: until st = s∗ 5: sβ0 ← s∗ . Second Phase 6: repeat 7: sβt+1 ← ∂Φ∂s ( x, sβt , θ ) − β ∂`∂s ( sβt , y
) 8: . Parameter Update at Time t 9: θ ← θ + ηβ ( ∂Φ ∂θ ( sβt+1 ) − ∂Φ∂θ ( sβt
)) 10: until sβt and θ are converged.
Figure 1: Left. Pseudo-code of EP. This is the version of EP for discrete-time dynamics introduced in Ernoult et al. (2019). Right. Pseudo-code of C-EP with simplified notations (see section 3.2 for a formal definition of C-EP). Difference between EP and C-EP. In EP, one global parameter update is performed at the end of the second phase ; in C-EP, parameter updates are performed throughout the second phase. Eq. 5 shows that the continual updates of C-EP add up to the global update of EP.

heading: 3.2 DESCRIPTION OF THE C-EP ALGORITHM
text: The first phase of C-EP is the same as that of EP (see Fig. 1). In the second phase of C-EP the parameter variable is regarded as another dynamic variable θt that evolves with time t along with st. The dynamics of st and θt in the second phase of C-EP depend on the values of the two hyperparameters β (the hyperparameter of influence) and η (the learning rate), therefore we write sβ,ηt and θβ,ηt to show explicitly this dependence. With now both the neurons and the synapses evolving in the second phase, the dynamic variables sβ,ηt and θ β,η t start from s β,η 0 = s∗ and θ β,η 0 = θ and follow:
∀t ≥ 0 :
 sβ,ηt+1 = ∂Φ ∂s ( x, sβ,ηt , θ β,η t ) − β ∂` ∂s ( sβ,ηt , y ) ,
θβ,ηt+1 = θ β,η t +
η
β
( ∂Φ
∂θ
( x, sβ,ηt+1, θ β,η t ) − ∂Φ ∂θ ( x, sβ,ηt , θ β,η t )) .
(6)
The difference in C-EP compared to EP is that the value of the parameter used to update sβ,ηt+1 in Eq. (6) is the current θβ,ηt , not θ. Provided the learning rate η is small enough, i.e. the synapses are slow compared to the neurons, this effect is weak. Intuitively, in the limit η → 0, the parameter changes are negligible so that θβ,ηt can be approximated by its initial value θ β,η 0 = θ. Under this approximation, the dynamics of sβ,ηt in C-EP and the dynamics of s β t in EP are the same. See Fig. 3 for a simple example, and Appendix A.3 for a proof in the general case.

heading: 3.3 GRADIENT-DESCENDING DYNAMICS (GDD)
text: Now we prove that, provided the hyperparameter β and the learning rate η are small enough, the dynamics of the neurons and the weights given by Eq. (6) follow the gradients of BPTT (Theorem 1 and Fig. 2). For a formal statement of this property, we define the normalized (continual) updates of C-EP, as well as the gradients of the loss L = ` (sT , y) after T time steps, computed with BPTT:
 ∆C−EPs (β, η, t) = 1 β ( sβ,ηt+1 − s β,η t ) ,
∆C−EPθ (β, η, t) = 1
η
( θβ,ηt+1 − θ β,η t ) ,
 ∇BPTTs (t) = ∂L ∂sT−t ,
∇BPTTθ (t) = ∂L
∂θT−t .
(7)
More details about L and the gradients∇BPTTs (t) and∇BPTTθ (t) are provided in Appendix B. Note that injecting Eq. (6) in Eq. (7), the normalized updates of C-EP consequently read:
∆C−EPθ (β, η, t) = 1
β
( ∂Φ
∂θ
( x, sβ,ηt+1, θ β,η t ) − ∂Φ ∂θ ( x, sβ,ηt , θ β,η t )) , (8)
which corresponds to the parameter gradient at time t, defined informally in Eq. (5). The following result makes this statement more formal. Theorem 1 (GDD Property). Let s0, s1, . . . , sT be the convergent sequence of states and denote s∗ = sT the steady state. Further assume that there exists some step K where 0 < K ≤ T such that s∗ = sT = sT−1 = . . . sT−K . Then, in the limit η → 0 and β → 0, the first K normalized updates in the second phase of C-EP are equal to the negatives of the first K gradients of BPTT, i.e.
∀t = 0, 1, . . . ,K :  limβ→0 limη→0 ∆ C−EP s (β, η, t) = −∇BPTTs (t),
lim β→0 lim η→0
∆C−EPθ (β, η, t) = −∇ BPTT θ (t).
(9)
Theorem 1 rewrites sβ,ηt+1 ≈ s β,η t − β ∂L∂sT−t and θ β,η t+1 ≈ θ β,η t − η ∂L∂θT−t , showing that in the second pĥase of C-EP, neurons and synapses descend the gradients of the loss L obtained with BPTT, with the hyperparameters β and η playing the role of learning rates for sβ,ηt and θ β,η t , respectively. Fig. 3 illustrates Theorem 1 with a simple dynamical system for which the normalized updates ∆C−EP and the gradients∇BPTT are analytically tractable - see Appendix C for derivation details.

heading: 4 NUMERICAL EXPERIMENTS
text: In this section, we validate our continual version of Equilibrium Propagation against training on the MNIST data set with two models. The first model is a vanilla RNN with tied and symmetric weights: the dynamics of this model approximately derive from a primitive function, which allows training with C-EP. The second model is a Discrete-Time RNN with untied and asymmetric weights, which is therefore closer to biology. We train this second model with a modified version of C-EP which we call C-VF (Continual Vector Field) as it is inspired from the algorithm with Vector-Field dynamics of Scellier et al. (2018). Ernoult et al. (2019) showed with simulations the intuitive result that, if a model is such that the normalized updates of EP ‘match’ the gradients of BPTT (i.e. if they are approximately equal), then the model trained with EP performs as well as the model trained with BPTT. Along the same lines, we show in this work that the more the EP normalized updates follow the gradients of BPTT before training, the best is the resulting training performance. We choose to implement C-EP and C-VF on vanilla RNNs to accelerate simulations (Ernoult et al., 2019).

heading: 4.1 MODELS: VANILLA RNNS WITH SYMMETRIC AND ASYMMETRIC WEIGHTS
text: Vanilla RNN with symmetric weights trained by C-EP. The first phase dynamics is defined as:
st+1 = σ (W · st +Wx · x) , (10) where σ is an activation function, W is a symmetric weight matrix connecting the layers s and Wx is a matrix connecting the input x to the layers s. Although the dynamics are not directly defined in terms of a primitive function, note that st+1 ≈ ∂Φ∂s (st,W ) with Φ(s,W ) = 1 2s > ·W · s if we ignore the activation function σ. Following Eq. (6) and Eq. (7), we define the normalized updates of this model as:
∆C−EPs (β, η, t) = 1
β
( sβ,ηt+1 − s β,η t ) , ∆C−EPW (β, η, t) = 1
β
( sβ,η >
t+1 · s β,η t+1 − s
β,η> t · s β,η t ) .
(11)
Note that this model applies to any topology as long as existing connections have symmetric values: this includes deep networks with any number of layers – see Appendix E for detailed descriptions of the models used in the experiments. More explicitly, for a network whose layers of neurons are s0, s1, ..., sN , with Wn,n+1 connecting the layers sn+1 and sn in both directions, the corresponding primitive function is Φ = ∑ n(s n)> ·Wn,n+1 · sn+1 + sN > ·Wx · x - see Appendix E.1 for more details.
Vanilla RNN with asymmetric weights trained by C-VF. In this model, the dynamics in the first phase is the same as Eq. (10) but now the weight matrix W is no longer assumed to be symmetric, i.e. the reciprocal connections between neurons are not constrained. In this setting the weight dynamics in the second phase is replaced by a version for asymmetric weights: W β,ηt+1 = W β,ηt + η β s β,η> t · ( sβ,ηt+1 − s β,η t ) , so that the normalized updates are equal to:
∆C−VFs (β, η, t) = 1
β
( sβ,ηt+1 − s β,η t ) , ∆C−VFW (β, η, t) = 1
β sβ,η
> t · ( sβ,ηt+1 − s β,η t ) . (12)
Like the previous model, the vanilla RNN with asymmetric weights also applies to deep networks with any number of layers. Although in C-VF the dynamics of the weights is not one of the form of Eq. (6) that derives from a primitive function, the (bioplausible) normalized weight updates of Eq. (12) can approximately follow the gradients of BPTT, provided that the values of reciprocal connections are not too dissimilar: this is illustrated in Fig. 5 (as well as in Fig. 12 and Fig. 13 of Appendix E.6) and proved in Appendix D.2. This property motivates the following training experiments.

heading: 4.2 C-EP TRAINING EXPERIMENTS
text: Experiments are performed with multi-layered vanilla RNNs (with symmetric weights) on MNIST. The table of Fig. 4.1 presents the results obtained with C-EP training benchmarked against standard EP training (Ernoult et al., 2019) - see Appendix E for model details and Appendix F.1 for training conditions. Although the test error of C-EP approaches that of EP, we observe a degradation in accuracy. This is because although Theorem 1 guarantees Gradient Descending Dynamics (GDD) in the limit of infinitely small learning rates, in practice we have to strike a balance between having a learning rate that is small enough to ensure this condition but not too small to observe convergence within a reasonable number of epochs. As seen on Fig. 5 (b), the finite learning rate η of continual updates leads to ∆C−EP(β, η, t) curves splitting apart from the −∇BPTT(t) curves. As seen per Fig. 5 (a), this effect is emphasized with the depth: before training, angles between the normalized updates of C-EP and the gradients of BPTT reach 50 degrees for two hidden layers. The deeper the network, the more difficult it is for the C-EP dynamics to follow the gradients provided by BPTT. As an evidence, we show in Appendix F.2 that when we use extremely small learning rates throughout the second phase (θ ← θ + ηtiny∆C−EPθ ) and rescale up the resulting total weight update (θ ← θ −∆θtot + ηηtiny ∆θtot), we recover standard EP results.

heading: 4.3 CONTINUAL VECTOR FIELD (C-VF) TRAINING EXPERIMENTS
text: Depending on whether the updates occur continuously during the second phase and the system obey general dynamics with untied forward and backward weights, we can span a large range of deviations from the ideal conditions of Theorem 1. Fig. 5 (b) qualitatively depicts these deviations with a model for which the normalized updates of EP match the gradients of BPTT (EP) ; with continual weight
updates, the normalized updates and gradients start splitting apart (C-EP), and even more so if the weights are untied (C-VF).
Protocol. In order to create these deviations from Theorem 1 and study the consequences in terms of training, we proceed as follows. For each C-VF simulations, we tune the initial angle between forward weights (θf ) and backward weights (θb) between 0 and 180◦. We denote this angle Ψ(θf , θb) - see Appendix F.1 for the angle definition and the angle tuning technique employed. For each of these weight initialization, we compute the angle between the total normalized update provided by C-VF, i.e. ∆C−VF(β, η, tot) = ∑K−1 t=0 ∆
C−VF(β, η, t) and the total gradient provided by BPTT, i.e. ∇BPTT(tot) = ∑K−1 t=0 ∇BPTT(t) on random mini-batches before training. We denote this
angle Ψ ( ∆C−VF(tot),−∇BPTT(tot) ) . Finally for each weight initialization, we perform training in the discrete-time setting of Ernoult et al. (2019) - see Appendix E.4 for model details. We proceed in the same way for EP and C-EP simulations, computing Ψ ( ∆EP(tot),−∇BPTT(tot) ) and
Ψ ( ∆C−EP(tot),−∇BPTT(tot) ) before training. We use the generic notation ∆(tot) to denote the
total normalized update. This procedure yields (x, y) data points with x = Ψ ( ∆(tot),−∇BPTT(tot) ) and y = test error, which are reported on Fig. 5 (a) - see Appendix F.1 for the full table of results.
Results. Fig. 5 (a) shows the test error achieved on MNIST by EP, C-EP on the vanilla RNN with symmetric weights and C-VF on the vanilla RNN with asymmetric weights for different number of hidden layers as a function of the angle Ψ ( ∆(tot),−∇BPTT(tot) ) before training. This graphical representation spreads the algorithms between EP which best satisfies the GDD property (leftmost point in green at∼ 20◦) to C-VF which satisfies the less the GDD property (rightmost points in red and orange at ∼ 100◦). As expected, high angles between gradients of C-VF and BPTT lead to high error rates that can reach 90% for Ψ ( ∆C−VF(tot),−∇BPTT(tot) ) over 100◦. More precisely, the inset of Fig. 5 shows the same data but focusing only on results generated by initial weight angles lying below 90◦, i.e. Ψ(θf , θb) = {0◦, 22.5◦, 45◦, 67.5◦, 90◦}. From standard EP with one hidden layer to C-VF with two hidden layers, the test error increases monotonically with Ψ ( ∆(tot),−∇BPTT(tot) ) but
does not exceed 5.05% on average. This result confirms the importance of proper weight initialization when weights are untied, also discussed in other context (Lillicrap et al., 2016). When the initial weight angle is of 0◦, the impact of untying the weights on classification accuracy remains constrained, as shown in table of Fig. 4.1. Upon untying the forward and backward weights, the test error increases by ∼ 0.2% with one hidden layer and by ∼ 0.5% with two hidden layers compared to standard C-EP.

heading: 5 DISCUSSION
text: Equilibrium Propagation is an algorithm that leverages the dynamical nature of neurons to compute weight gradients through the physics of the neural network. C-EP embraces simultaneous synapse and neuron dynamics, resolving the initial need of artificial memory units for storing the neuron values between different phases. The C-EP framework preserves the equivalence with Backpropagation Through Time: in the limit of sufficiently slow synaptic dynamics (i.e. small learning rates), the system satisfies Gradient Descending Dynamics (Theorem 1). Our experimental results confirm this theorem. When training our vanilla RNN with symmetric weights with C-EP while ensuring convergence in 100 epochs, a modest reduction in MNIST accuracy is seen with regards to standard EP. This accuracy reduction can be eliminated by using smaller learning rates and rescaling up the total weight update at the end of the second phase (Appendix F.2). On top of extending the theory of Ernoult et al. (2019), Theorem 1 also appears to provide a statistically robust tool for C-EP based learning. Our experimental results show as in Ernoult et al. (2019) that, for a given network with specified neuron and synapse dynamics, the more the updates of Equilibrium Propagation follow the gradients provided by Backpropagation Through Time before training (in terms of angle in this work), the better this network can learn.
Our C-EP and C-VF algorithms exhibit features reminiscent of biology. C-VF extends C-EP training to RNNs with asymmetric weights between neurons, as is the case in biology. Its learning rule, local in space and time, is furthermore closely acquainted to Spike Timing Dependent Plasticity (STDP), a learning rule widely studied in neuroscience, inferred in vitro and in vivo from neural recordings in the hippocampus (Dan and Poo, 2004). In STDP, the synaptic strength is modulated by the relative timings of pre and post synaptic spikes within a precise time window (Bi and Poo, 1998; 2001).
θ (t), ∆
C−VF
θ (t)) and the gradients −∇BPTTθ (t).
Each randomly selected synapse corresponds to one color. While dashed and continuous lines coincide for standard EP, they split apart upon untying the weights and using continual updates.
Strikingly, the same rule that we use for C-VF learning can approximate STDP correlations in a rate-based formulation, as shown through numerical experiments by Bengio et al. (2015). From this viewpoint our work brings EP a step closer to biology. However, C-EP and C-VF do not aim at being models of biological learning per se, in that it would account for how the brain works or how animals learn, for which Reinforcement Learning might be a more suited learning paradigm. The core motivation of this work is to propose a fully local implementation of EP, in particular to foster its hardware implementation.
When computed on a standard computer, due to the use of small learning rates to mimic analog dynamics within a finite number of epochs, training our models with C-EP and C-VF entail long simulation times. With a Titan RTX GPU, training a fully connected architecture on MNIST takes 2 hours 39 mins with 1 hidden layer and 10 hours 49 mins with 2 hidden layers. On the other hand, C-EP and C-VF might be particularly efficient in terms of speed and energy consumption when operated on neuromorphic hardware that employs analog device physics (Ambrogio et al., 2018; Romera et al., 2018). To this purpose, our work can provide an engineering guidance to map our algorithm onto a neuromorphic system. Fig. 5 (a) shows that hyperparameters should be tuned so that before training, C-EP updates stay within 90◦ of the gradients provided by BPTT. More concretely in practice, it amounts to tune the degree of symmetry of the dynamics, for instance the angle between forward and backward weights - see Fig. 4.1. Our work is one step towards bridging Equilibrium Propagation with neuromorphic computing and thereby energy efficient implementations of gradient-based learning algorithms.

heading: Appendix

heading: A PROOF OF THEOREM 1
text: In this appendix, we prove Theorem 1, which we recall here. Theorem 1 (GDD Property). Let s0, s1, . . . , sT be the convergent sequence of states and denote s∗ = sT the steady state. Further assume that there exists some step K where 0 < K ≤ T such that s∗ = sT = sT−1 = . . . sT−K . Then, in the limit η → 0 and β → 0, the first K normalized updates in the second phase of C-EP are equal to the negatives of the first K gradients of BPTT, i.e.
∀t = 0, 1, . . . ,K :  limβ→0 limη→0 ∆ C−EP s (β, η, t) = −∇BPTTs (t),
lim β→0 lim η→0
∆C−EPθ (β, η, t) = −∇ BPTT θ (t).
(9)

heading: A.1 A SPECTRUM OF FOUR COMPUTATIONALLY EQUIVALENT LEARNING ALGORITHMS
text: Proving Theorem 1 amounts to prove the equivalence of C-EP and BPTT. In fact we can prove the equivalence of four algorithms, which all compute the gradient of the loss:
1. Backpropagation Through Time (BPTT), presented in Section B.2, 2. Recurrent Backpropagation (RBP), presented in Section B.3, 3. Equilibrium Propagation (EP), presented in Section 2, 4. Equilibrium Propagation with Continual Weight Updates (C-EP), introduced in Section 3.
In this spectrum of algorithms, BPTT is the most practical algorithm to date from the point of view of machine learning, but also the less biologically realistic. In contrast, C-EP is the most realistic in terms of implementation in biological systems, while it is to date the least practical and least efficient for conventional machine learning (computations on standard Von-Neumann hardware are considerably slower due to repeated parameter updates, requiring memory access at each time-step of the second phase).

heading: A.2 SKETCH OF THE PROOF
text: Theorem 1 can be proved in three phases, using the following three lemmas. Lemma 2 (Equivalence of C-EP and EP). In the limit of small learning rate, i.e. η → 0, the (normalized) updates of C-EP are equal to those of EP:
∀t ≥ 0 :  lim η→0 (η>0) ∆C−EPs (β, η, t) = ∆ EP s (β, t),
lim η→0 (η>0)
∆C−EPθ (β, η, t) = ∆ EP θ (β, t).
(13)
Lemma 3 (Equivalence of EP and RBP). Assume that the transition function derives from a primitive function, i.e. that F is of the form F (x, s, θ) = ∂Φ∂s (x, s, θ). Then, in the limit of small hyperparameter β, the normalized updates of EP are equal to the gradients of RBP:
∀t ≥ 0 :  lim β→0 (β>0) ∆EPs (β, t) = −∇RBPs (t),
lim β→0 (β>0)
∆EPθ (β, t) = −∇RBPθ (t). (14)
Lemma 4 (Equivalence of BPTT and RBP). In the setting with static input x, suppose that the network has reached the steady state s∗ after T −K steps, i.e. sT−K = sT−K+1 = · · · = sT−1 = sT = s∗. Then the first K gradients of BPTT are equal to the first K gradient of RBP, i.e.
∀t = 0, 1, . . . ,K : { ∇BPTTs (t) = ∇RBPs (t), ∇BPTTθ (t) = ∇RBPθ (t).
(15)
Proofs of the Lemmas can be found in the following places:
• The link between BPTT and RBP (Lemma 2) is known since the late 1980s and can be found e.g. in Hertz (2018). We also prove it here in Appendix B. • Lemma 3 was proved in Scellier and Bengio (2019) in the setting of real-time dynamics. • Lemma 4 is the new ingredient contributed here, and we prove it in Appendix A.3.
Also a direct proof of the equivalence of EP and BPTT was derived in Ernoult et al. (2019).

heading: A.3 EQUIVALENCE OF C-EP AND EP
text: First, recall the dynamics of C-EP in the second phase: starting from sβ,η0 = s∗ and θ β,η 0 = θ we have
∀t ≥ 0 :  sβ,ηt+1 = ∂Φ ∂s ( x, sβ,ηt , θ β,η t ) − β ∂` ∂s ( sβ,ηt , y ) ,
θβ,ηt+1 = θ β,η t +
η
β
( ∂Φ
∂θ
( x, sβ,ηt+1, θ β,η t ) − ∂Φ ∂θ ( x, sβ,ηt , θ β,η t )) .
(16)
We have also defined the normalized updates of C-EP:
∀t ≥ 0 :  ∆C−EPs (β, η, t) = 1 β ( sβ,ηt+1 − s β,η t ) ,
∆C−EPθ (β, η, t) = 1
η
( θβ,ηt+1 − θ β,η t ) .
(17)
We also recall the dynamics of EP in the second phase:
sβ0 = s∗ and s β t+1 =
∂Φ
∂s
( x, sβt , θ ) − β ∂`
∂s
( sβt , y ) , (18)
as well as the normalized updates of EP, as defined in Ernoult et al. (2019):
∀t ≥ 0 :  ∆EPs (β, t) = 1 β ( sβt+1 − s β t ) ,
∆EPθ (β, t) = 1
β
( ∂Φ
∂θ
( x, sβt+1, θ ) − ∂Φ ∂θ ( x, sβt , θ )) .
(19)
Lemma 2 (Equivalence of C-EP and EP). In the limit of small learning rate, i.e. η → 0, the (normalized) updates of C-EP are equal to those of EP:
∀t ≥ 0 :  lim η→0 (η>0) ∆C−EPs (β, η, t) = ∆ EP s (β, t),
lim η→0 (η>0)
∆C−EPθ (β, η, t) = ∆ EP θ (β, t).
(13)
Proof of Lemma 2. We want to compute the limits of ∆C−EPs (β, η, t) and ∆ C−EP θ (β, η, t) as η → 0 with η > 0. First of all, note that under mild assumptions – which we made here – of regularity on the functions Φ and ` (e.g. continuous differentiability), for fixed t and β, the quantities sβ,ηt and θ β,η t are continuous as functions of η ; this is straightforward from the form of Eq. (16). As a consequence, ∆C−EPs (β, η, t) is a continuous function of η, which implies in particular that
lim η→0 (η>0)
∆C−EPs (β, η, t) = ∆ C−EP s (β, 0, t). (20)
Now, taking η = 0 in the bottom equation of Eq. (16) yields the recurrence relation θβ,0t+1 = θ β,0 t , so that θβ,0t = θ β,0 0 = θ for every t. Injecting θ β,0 t = θ in the top equation of Eq. (16) yields for s β,0 t the same recurrence relation as that of sβt (Eq. 18), so that s β,0 t = s β t for every t. Therefore, for η = 0, we have ∆C−EPs (β, 0, t) = 1
β
( sβ,0t+1 − s β,0 t ) = 1
β
( sβt+1 − s β t ) = ∆EPs (β, t). (21)
It follows from Eq. (20) and Eq. (21) that
lim η→0 (η>0)
∆C−EPs (β, η, t) = ∆ EP s (β, t). (22)
Now let us compute limη→0 (η>0) ∆C−EPs (β, η, t). Using Eq. (16), we have
∆C−EPθ (β, η, t) = 1
η
( θβ,ηt+1 − θ β,η t ) (23)
= 1
β
( ∂Φ
∂θ
( x, sβ,ηt+1, θ β,η t ) − ∂Φ ∂θ ( x, sβ,ηt , θ β,η t )) . (24)
Similarly as before, for fixed t, ∂Φ∂θ ( x, sβ,ηt , θ β,η t ) is a continuous function of η. Therefore
lim η→0 (η>0)
∆C−EPθ (β, η, t) = 1
β
( ∂Φ
∂θ
( x, sβ,0t+1, θ β,0 t ) − ∂Φ ∂θ ( x, sβ,0t , θ β,0 t )) (25)
= 1
β
( ∂Φ
∂θ
( x, sβt+1, θ ) − ∂Φ ∂θ ( x, sβt , θ )) = ∆EPθ (β, t). (26)
A consequence of Lemma 2 is that the total update of C-EP matches the total update of EP in the limit of small η, so that we retrieve the standard EP learning rule of Eq. (4). More explicitly, after K steps in the second phase and starting from θβ,η0 = θ0:
θβ,ηK − θ0 = K−1∑ t=0 θβ,ηt+1 − θ β,η t
= K−1∑ t=0 η∆C−EPθ (β, η, t) (definition of ∆ C−EP θ , Eq. (7))
≈ η&0 K−1∑ t=0 η∆EPθ (β, t) (Lemma 2)
= K−1∑ t=0 η 1 β ( ∂Φ ∂θ (x, st+1, θ0)− ∂Φ ∂θ (x, st+1, θ0) ) (definition of ∆EPθ , Eq. (19))
= η
β
( ∂Φ
∂θ (x, sβK , θ0)−
∂Φ ∂θ (x, s∗, θ0)
)

heading: B EQUIVALENCE OF BPTT AND RBP
text: In this section, we recall Backprop Through Time (BPTT) and the Almeida-Pineda Recurrent Backprop (RBP) algorithm, which can both be used to optimize the loss L∗ of Eq. 3. Historically, BPTT and RBP were invented separately around the same time. RBP was introduced at a time when convergent RNNs (such as the one studied in this paper) were popular. Nowadays, convergent RNNs are less popular ; in the field of deep learning, RNNs are almost exclusively used for tasks that deal with sequential data and BPTT is the algorithm of choice to train such RNNs. Here, we present RBP in a way that it can be seen as a particular case of BPTT.

heading: B.1 SKETCH OF THE PROOF OF LEMMA 4
text: Lemma 4, which we recall here, is a consequence of Proposition 5 and Definition 6 below. Lemma 4 (Equivalence of BPTT and RBP). In the setting with static input x, suppose that the network has reached the steady state s∗ after T −K steps, i.e. sT−K = sT−K+1 = · · · = sT−1 = sT = s∗. Then the first K gradients of BPTT are equal to the first K gradient of RBP, i.e.
∀t = 0, 1, . . . ,K : { ∇BPTTs (t) = ∇RBPs (t), ∇BPTTθ (t) = ∇RBPθ (t).
(15)
However, in general, the gradients∇BPTT(t) of BPTT and the gradients∇RBP(t) of RBP are not equal for t > K. This is because BPTT and RBP compute the gradients of different loss functions:
• BPTT computes the gradient of the loss after T time steps, i.e. L = ` (sT , y),
• RBP computes the gradients of the loss at the steady state, i.e. L∗ = ` (s∗, y).

heading: B.2 BACKPROPAGATION THROUGH TIME (BPTT)
text: Backpropagation Through Time (BPTT) is the standard method to train RNNs and can also be used to train the kind of convergent RNNs that we study in this paper. To this end, we consider the cost of the state sT after T time steps, denoted L = ` (sT , y), and we substitute the loss after T time steps L as a proxy for the loss at the steady state L∗ = ` (s∗, y). The gradients of L can then be computed with BPTT.
To do this, we recall some of the inner working mechanisms of BPTT. Eq. (1) rewrites in the form st+1 = F (x, st, θt+1 = θ), where θt denotes the parameter of the model at time step t, the value θ being shared across all time steps. This way of rewriting Eq. (1) enables us to define the partial derivative ∂L∂θt as the sensitivity of the loss L with respect to θt when θ1, . . . θt−1, θt+1, . . . θT remain fixed (set to the value θ). With these notations, the gradient ∂L∂θ reads as the sum:
∂L ∂θ = ∂L ∂θ1 + ∂L ∂θ2 + · · ·+ ∂L ∂θT . (27)
BPTT computes the ’full’ gradient ∂L∂θ by first computing the partial derivatives ∂L ∂st and ∂L∂θt iteratively, backward in time, using the chain rule of differentiation. In this work, we denote the gradients that BPTT computes:
∀t ∈ [0, T − 1] :  ∇BPTTs (t) = ∂L ∂sT−t ,
∇BPTTθ (t) = ∂L
∂θT−t .
(28)
Proposition 5 (Gradients of BPTT). The gradients ∇BPTTs (t) and ∇BPTTθ (t) satisfy the recurrence relationship
∇BPTTs (0) = ∂`
∂s (sT , y) , (29)
∀t = 1, 2, . . . , T, ∇BPTTs (t) = ∂F
∂s (x, sT−t, θ)
> · ∇BPTTs (t− 1), (30)
∀t = 1, 2, . . . , T, ∇BPTTθ (t) = ∂F
∂θ (x, sT−t, θ)
> · ∇BPTTs (t− 1). (31)

heading: B.3 FROM BACKPROP THROUGH TIME (BPTT) TO RECURRENT BACKPROP (RBP)
text: In general, to apply BPTT, it is necessary to store in memory the history of past hidden states s1, s2, . . . , sT in order to compute the gradients∇BPTTs (t) and∇BPTTθ (t) as in Eq. 30-31. However, in our specific setting with static input x, if the network has reached the steady state s∗ after T −K steps, i.e. if sT−K = sT−K+1 = · · · = sT−1 = sT = s∗, then we see that, in order to compute the first K gradients of BPTT, all one needs to know is ∂F∂s (x, s∗, θ) and ∂F ∂θ (x, s∗, θ). To this end, all one needs to keep in memory is the steady state s∗. In this particular setting, it is not necessary to store the past hidden states sT , sT−1, . . . , sT−K since they are all equal to s∗.
The Almeida-Pineda algorithm (a.k.a. Recurrent Backpropagation, or RBP for short), which was invented independently by Almeida (1987) and Pineda (1987), relies on this property to compute the gradients of the loss L∗ using only the steady state s∗. Similarly to BPTT, it computes quantities ∇RBPs (t) and∇RBPθ (t), which we call ‘gradients of RBP’, iteratively for t = 0, 1, 2, . . .
Definition 6 (Gradients of RBP). The gradients ∇RBPs (t) and ∇RBPθ (t) are defined and computed iteratively as follows:
∇RBPs (0) = ∂`
∂s (s∗, y) , (32)
∀t ≥ 0, ∇RBPs (t+ 1) = ∂F
∂s (x, s∗, θ)
> · ∇RBPs (t), (33)
∀t ≥ 0, ∇RBPθ (t+ 1) = ∂F
∂θ (x, s∗, θ)
> · ∇RBPs (t). (34)
Unlike in BPTT where keeping the history of past hidden states is necessary to compute (or ‘backpropagate’) the gradients, in RBP Eq. 33-34 show that it is sufficient to keep in memory the steady state s∗ only in order to iterate the computation of the gradients. RBP is more memory efficient than BPTT.
Algorithm 3 BPTT Input: x, y, θ. Output: θ.
1: s0 ← 0 2: for t = 0 to T − 1 do 3: st+1 ← F (x, st, θ) 4: end for 5: ∇BPTTs (0)← ∂`∂s (sT , y) 6: for t = 1 to T do 7: ∇s(t)← ∂F∂s (x, sT−t, θ)
> · ∇s(t− 1) 8: ∇θ(t)← ∂F∂θ (x, sT−t, θ)
> · ∇s(t− 1) 9: end for
10: ∇BPTTθ (tot)← ∑T−1 t=0 ∇BPTTθ (t)
Algorithm 4 RBP Input: x, y, θ. Output: θ.
1: s0 ← 0 2: repeat 3: st+1 ← F (x, st, θ) 4: until st = s∗ 5: ∇RBPs (0)← ∂`∂s (s∗, y) 6: repeat 7: ∇RBPs (t)← ∂F∂s (x, s∗, θ)
> ·∇RBPs (t−1) 8: ∇RBPθ (t)← ∂F∂θ (x, s∗, θ)
> ·∇RBPs (t−1) 9: until∇RBPθ (t) = 0.
10: ∇RBPθ (tot)← ∑∞ t=0∇RBPθ (t)
Figure 6: Left. Pseudo-code of BPTT. The gradients ∇(t) denote the gradients ∇BPTT(t) of BPTT. Right. Pseudo-code of RBP. Difference between BPTT and RBP. In BPTT, the state sT−t is required to compute ∂F∂s (x, sT−t, θ) and ∂F ∂θ (x, sT−t, θ) ; thus it is necessary to store in memory the sequence of states s1, s2, . . . , sT . In contrast, in RBP, only the steady state s∗ is required to compute ∂F ∂s (x, s∗, θ) and ∂F ∂θ (x, s∗, θ) ; it is not necessary to store the past states of the network.

heading: B.4 WHAT ‘GRADIENTS’ ARE THE GRADIENTS OF RBP?
text: In this subsection we motivate the name of ‘gradients’ for the quantities∇RBPs (t) and∇RBPθ (t) by proving that they are the gradients of L∗ in the sense of Proposition 7 below. They are also the gradients of what we call the ‘projected cost function’ (Proposition 8), using the terminology of Scellier and Bengio (2019). Proposition 7 (RBP Optimizes L∗). The total gradient computed by the RBP algorithm is the gradient of the loss L∗ = ` (s∗, y), i.e.
∞∑ t=1 ∇RBPθ (t) = ∂L∗ ∂θ . (35)
∇RBPs (t) and ∇RBPθ (t) can also be expressed as gradients of Lt = ` (st, y), the cost after t time steps. In the terminology of Scellier and Bengio (2019), Lt was named the projected cost. For t = 0, L0 is simply the cost of the initial state s0. For t > 0, Lt is the cost of the state projected a duration t in the future. Proposition 8 (Gradients of RBP are Gradients of the Projected Cost). The ‘RBP gradients’∇RBPs (t) and ∇RBPθ (t) can be expressed as gradients of the projected cost:
∀t ≥ 0, ∇RBPs (t) = ∂Lt ∂s0 ∣∣∣∣ s0=s∗ , ∇RBPθ (t) = ∂Lt ∂θ0 ∣∣∣∣ s0=s∗
(36)
where the initial state s0 is the steady state s∗.
Proof of Proposition 7. First of all, by Definition 6 (Eq. 32-34) it is straightforward to see that
∀t ≥ 0, ∇RBPs (t) = ( ∂F
∂s (x, s∗, θ) > )t · ∂` ∂s (s∗, y) , (37)
∀t ≥ 1, ∇RBPθ (t) = ∂F
∂θ (x, s∗, θ)
> · ( ∂F
∂s (x, s∗, θ)
> )t−1
· ∂` ∂s (s∗, y) . (38)
Second, recall that the loss L∗ is L∗ = ` (s∗, y) , (39)
where s∗ = F (x, s∗, θ) . (40) By the chain rule of differentiation, the gradient of L∗ (Eq. 39) is ∂L∗
∂θ = ∂` ∂s (s∗, y) · ∂s∗ ∂θ . (41)
In order to compute ∂s∗∂θ , we differentiate the steady state condition (Eq. 40) with respect to θ, which yields
∂s∗ ∂θ = ∂F ∂s (x, s∗, θ) · ∂s∗ ∂θ + ∂F ∂θ (x, s∗, θ) . (42)
Rearranging the terms, and using the Taylor expansion (Id−A)−1 = ∑∞ t=0A
t with A = ∂F ∂s (x, s∗, θ), we get
∂s∗ ∂θ =
( Id− ∂F
∂s (x, s∗, θ) )−1 · ∂F ∂θ (x, s∗, θ) (43)
= ∞∑ t=0 ( ∂F ∂s (x, s∗, θ) )t · ∂F ∂θ (x, s∗, θ) . (44)
Therefore ∂L∗
∂θ = ∂` ∂s (s∗, y) · ∂s∗ ∂θ
(45)
= ∞∑ t=0 ∂` ∂s (s∗, y) · ( ∂F ∂s (x, s∗, θ) )t · ∂F ∂θ (x, s∗, θ) (46)
= ∞∑ t=0 ∇RBPθ (t). (47)
Proof of Proposition 8. By the chain rule of differentiation we have ∂Lt+1 ∂s0 = ∂F ∂s (x, s0, θ) > · ∂Lt+1 ∂s1 . (48)
Evaluation this expression for s0 = s∗ we get
∂Lt+1 ∂s0 ∣∣∣∣ s0=s∗ = ∂F ∂s (x, s∗, θ) > · ∂Lt+1 ∂s1 ∣∣∣∣ s0=s∗ . (49)
Finally note that ∂Lt+1 ∂s1 ∣∣∣∣ s0=s∗ = ∂Lt+1 ∂s1 ∣∣∣∣ s1=s∗ = ∂Lt ∂s0 ∣∣∣∣ s0=s∗
(50)
Therefore ∂Lt∂s0 ∣∣∣ s0=s∗ and∇RBPs (t) satisfy the same recurrence relation, thus they are equal. Proving the equality of ∂Lt∂θ0 ∣∣∣ s0=s∗ and∇RBPθ (t) is analogous.
C ILLUSTRATING THE EQUIVALENCE OF THE FOUR ALGORITHMS ON AN ANALYTICALLY TRACTABLE MODEL
Model. To illustrate the equivalence of the four algorithms (BPTT, RBP, EP and CEP), we study a simple model with scalar variable s and scalar parameter θ:
s0 = 0, st+1 = 1
2 (st + θ) , L∗ =
1 2 s2∗, (51)
where s∗ is the steady state of the dynamics (it is easy to see that the solution is s∗ = θ). The dynamics rewrites st+1 = F (st, θ) with the transition function F (s, θ) = 12 (s+ θ), and the loss rewrites L∗ = ` (s∗) with the cost function `(s) = 12s
2. Furthermore, a primitive function of the system 1 is Φ(s, θ) = 14 (s + θ)
2. This model has no practical application ; it is only meant for pedagogical purpose.
Backpropagation Through Time (BPTT). With BPTT, an important point is that we approximate the steady state s∗ by the state after T time steps sT , and we approximate L∗ (the loss at the steady state) by the loss after T time steps L = ` (sT ). In order to compute (i.e. ‘backpropagate’) the gradients of BPTT, Proposition 5 tells us that we need to compute ∂`∂s (sT ) = sT , ∂F ∂s (st, θ) = 1 2 and ∂F ∂θ (st, θ) = 1 2 . We get
∀t = 0, 1, . . . , T − 1, ∇BPTTs (t) = sT 2t , ∇BPTTθ (t) = sT 2t+1 . (52)
Recurrent Backpropagation (RBP). Similarly, to compute the gradients of RBP, Definition 6 tells us that we need to compute ∂`∂s (s∗) = s∗, ∂F ∂s (s∗, θ) = 1 2 and ∂F ∂θ (s∗, θ) = 1 2 . We have
∀t ≥ 0, ∇RBPs (t) = s∗ 2t , ∇RBPθ (t) = s∗ 2t+1 . (53)
The state after T time steps in BPTT converges to the steady state s∗ as T → ∞, therefore the gradients of BPTT converge to the gradients of RBP. Also notice that the steady state of the dynamics is s∗ = θ.
Equilibrium Propagation (EP). Following the equations governing the second phase of EP (Fig. 1), we have:
sβ0 = θ, s β t+1 =
( 1
2 − β
) sβt + 1
2 θ. (54)
This linear dynamical system can be solved analytically:
∀t ≥ 0, sβt = θ
1 + 2β
( 1 + 2β ( 1
2 − β
)t) . (55)
Notice that sβt → θ as β → 0 ; for small values of the hyperparameter β, the trajectory in the second phase is close to the steady state s∗ = θ.
Using Eq. 19, it follows that the normalized updates of EP are
∀t ≥ 0, ∆EPs (β, t) = − θ
2t (1− 2β)t , ∆EPθ (β, t) = −
θ
2t+1 (1− 2β)t . (56)
Notice again that the normalized updates of EP converge to the gradients of RBP as β → 0.
Continual Equilibrium Propagation (C-EP). The system of equations governing the system is:
{ sβ,η0 = s∗,
θβ,η0 = θ, ∀t ≥ 0 :  sβ,ηt+1 = ( 1 2 − β ) sβ,ηt + 1 2 θβ,ηt ,
θβ,ηt+1 = θ β,η t +
η
2β
( sβ,ηt+1 − s β,η t ) .
(57)
1The primitive function Φ is determined up to a constant.
First, rearranging the terms in the second equation, we get
1
η
( θβ,ηt+1 − θ β,η t ) = 1
2β
( sβ,ηt+1 − s β,η t ) . (58)
It follows that ∆C−EPθ (β, η, t) = 1
2 ∆C−EPs (β, η, t). (59)
Therefore, all we need to do is to compute ∆C−EPs (β, η, t). Second, by iterating the second equation over all indices from t = 0 to t− 1 we get
θβ,ηt = θ + η
2β
( sβ,ηt − s∗ ) . (60)
Using s∗ = θ and plugging this into the first equation we get
sβ,ηt+1 =
( 1
2 − β + η 4β
) sβ,ηt + ( 1
2 − η 4β
) θ. (61)
Solving this linear dynamical system, and using the initial condition sβ,η0 = θ we get
sβ,ηt = θ
1− η2β + 2β
[ 1− η
2β + 2β
( 1
2
)t( 1− 2β + η
2β
)t] (62)
Finally:
∆C−EPs (β, η, t) = − θ
2t
( 1− 2β + η
2β
)t (63)

heading: D COMPLEMENT ON GRADIENT-DESCENDING DYNAMICS (GDD)
text: Step-by-step equivalence of the dynamics of EP and gradient computation in BPTT was shown in Ernoult et al. (2019) and was refered to as the Gradient-Descending Updates (GDU) property. In this appendix, we first explain the connection between the GDD property of this work and the GDU property of Ernoult et al. (2019). Then we prove another version of the GDD property (Theorem 9 below), more general than Theorem 1.

heading: D.1 GRADIENT-DESCENDING UPDATES (GDU) OF ERNOULT ET AL. (2019)
text: The GDU property of Ernoult et al. (2019) states that the (normalized) updates of EP are equal to the gradients of BPTT. Similarly, the Gradient-Descending Dynamics (GDD) property of this work states that the normalized updates of C-EP are equal to the gradients of BPTT. The difference between the GDU property and the GDD property is that the term ‘update’ has slightly different meanings in the contexts of EP and C-EP. In C-EP, the ‘updates’ are the effective updates by which the neuron and synapses are being dynamically updated throughout the second phase. In contrast in EP, the ‘updates’ are effectively performed at the end of the second phase.

heading: D.2 A GENERALISATION OF THE GDD PROPERTY
text: The Gradient Descending Dynamics property (GDD, Theorem 1) states that, when the system dynamics derive from a primitive function, i.e. when the transition function F is of the form F = ∂Φ∂s , then the normalized updates of C-EP match the gradients provided by BPTT. Remarkably, even in the case of the C-VF dynamics that do not derive from a primitive function Φ, Fig. 5 shows that the biologically realistic update rule of C-VF follows well the gradients of BPTT. More illustrations of this property are shown on Fig. 12 and Fig. 13. In this section we give a theoretical justification for this fact by proving a more general result than Theorem 1.
First, recall the dynamics of the C-VF model. In the first phase:
st+1 = σ (W · st) , (64)
where σ is an activation function and W is a square weight matrix. In the second phase, starting from sβ,η0 = s∗ and W β,η 0 = W , the dynamics read:
∀t ≥ 0 :
 sβ,ηt+1 = σ ( W β,ηt · s β,η t ) − β ∂` ∂s ( sβ,ηt ) ,
W β,ηt+1 = W β,η t +
η β sβ,η >
t · ( sβ,ηt+1 − s β,η t ) .
(65)
Now let us define the transition function F (s,W ) = σ(W · s), so that the dynamics of the first phase rewrites
st+1 = F (st,W ) . (66)
As for the second phase, notice that ∂F∂W (s,W ) = σ ′(W · s) · s, so that if we ignore the factor σ′(W · s), Eq. (65) rewrites
∀t ≥ 0 :
 sβ,ηt+1 = F ( sβ,ηt ,W β,η t ) − β ∂` ∂s ( sβ,ηt ) ,
W β,ηt+1 = W β,η t +
η
β
∂F
∂W
( sβ,ηt ,W β,η t )> · ( sβ,ηt+1 − s β,η t ) .
(67)
Now, recall the definition of the normalized updates of C-VF, as well as the gradients of the loss L = ` (sT , y) after T time steps, computed with BPTT: ∆C−VFs (β, η, t) = 1 β ( sβ,ηt+1 − s β,η t ) ,
∆C−VFW (β, η, t) = 1
η
( W β,ηt+1 −W β,η t ) ,
 ∇BPTTs (t) = ∂L ∂sT−t ,
∇BPTTW (t) = ∂L
∂WT−t .
(68)
The loss L and the gradients∇BPTTs (t) and ∇BPTTθ (t) are defined formally in Appendix B.2.
Theorem 9 (Generalisation of the GDD Property). Let s0, s1, . . . , sT be the convergent sequence of states and denote s∗ = sT the steady state. Further assume that there exists some step K where 0 < K ≤ T such that s∗ = sT = sT−1 = . . . sT−K . Finally, assume that the Jacobian of the transition function at the steady state is symmetric, i.e. ∂F∂s (s∗,W ) = ∂F ∂s (s∗,W )
>. Then, in the limit η → 0 and β → 0, the first K normalized updates of C-VF follow the the first K gradients of BPTT, i.e.
∀t = 0, 1, . . . ,K :  limβ→0 limη→0 ∆ C−VF s (β, η, t) = −∇BPTTs (t),
lim β→0 lim η→0
∆C−VFW (β, η, t) = −∇ BPTT W (t).
(69)
A few remarks need to be made:
1. Observe that ∂F
∂s (s,W ) = σ′(W · s) ·W>. (70)
Ignoring the factor σ′(W · s), we see that if W is symmetric then the Jacobian of F is also symmetric, in which case the conditions of Theorem 9 are met.
2. Theorem 1 is a special case of Theorem 9. To see why, notice that if the transition function F is of the form F (s,W ) = ∂Φ∂s (s,W ), then
∂F ∂s (s,W ) =
∂2Φ ∂s2 (s,W ) = ∂F ∂s (s,W )> (71)
In this case the extra assumption in Theorem 9 is automatically satisfied.

heading: D.3 PROOF OF THEOREM 9
text: Theorem 9 is a consequence of Proposition 5 (Appendix B.2), which we recall here, and Lemma 10 below. Proposition 5 (Gradients of BPTT). The gradients ∇BPTTs (t) and ∇BPTTθ (t) satisfy the recurrence relationship
∇BPTTs (0) = ∂`
∂s (sT , y) , (29)
∀t = 1, 2, . . . , T, ∇BPTTs (t) = ∂F
∂s (x, sT−t, θ)
> · ∇BPTTs (t− 1), (30)
∀t = 1, 2, . . . , T, ∇BPTTθ (t) = ∂F
∂θ (x, sT−t, θ)
> · ∇BPTTs (t− 1). (31)
Lemma 10 (Updates of C-VF). Define the (normalized) neural and weight updates of C-VF in the limit η → 0 and β → 0:
∀t = 0, 1, . . . ,K :  ∆ C−VF s (t) = lim β→0 lim η→0 ∆C−VFs (β, η, t),
∆C−VFθ (t) = lim β→0 lim η→0
∆C−VFθ (β, η, t). (72)
They satisfy the recurrence relationship
∆C−VFs (0) = − ∂`
∂s (s∗, y) , (73)
∀t ≥ 0, ∆C−VFs (t+ 1) = ∂F
∂s (x, s∗, θ) ·∆C−VFs (t), (74)
∀t ≥ 0, ∆C−VFθ (t+ 1) = ∂F
∂θ (x, s∗, θ)
> ·∆C−VFs (t). (75)
The proof of Lemma 10 is similar to the one provided in Ernoult et al. (2019).

heading: E MODELS
text: In this section, we describe the C-EP and C-VF algorithms when implemented on multi-layered models, with tied weights and untied weights respectively. In the fully connected layered architecture model, the neurons are only connected between two consecutive layers (no skip-layer connections and no lateral connections within a layer). We denote neurons of the n-th layer as sn with n ∈ [0, N − 1], where N is the number of hidden layers. Layers are labelled in a backward fashion: n = 0 labels the output layer, n = 1 the first hidden layer starting from the output layer, and n = N − 1 the last hidden layer (before the input layer). Thus, there are N hidden layers in total. Fig. 7 shows this architecture with N = 2. Each model are presented here in a "real-time" and "discrete-time" settings For each model we lay out the equations of the neuron and synapse dynamics, we demonstrate the GDD property and we specify in which part of the main text they are used.
We present in this order:
1. Discrete-Time RNN with symmetric weights trained with EP,
2. Discrete-Time RNN with symmetric weights trained with C-EP,
3. Real-Time RNN with symmetric weights trained with C-EP,
4. Discrete-Time RNN with asymmetric weights trained with C-VF,
5. Real-Time RNN with asymmetric weights trained with C-VF.
Our ‘Discrete-Time RNN model’ is also commonly called ‘vanilla RNN’, and it is refered to as the ‘prototypical model’ in Ernoult et al. (2019). Our ‘Real-Time RNN model with symmetric weights’ is also commonly called ‘continuous Hopfield model’ and is refered to as the ‘energy-based model’ in Ernoult et al. (2019).
Demonstrating the Gradient Descending Dynamics (GDD) property (Theorem 1) on MNIST. For this experiment, we consider the 784-512-. . . -512-10 network architecture, with 784 input neurons, 10 ouput neurons, and 512 neurons per hidden layer. The activation function used is σ(x) = tanh(x). The experiment consists of the following: we take a random MNIST sample (of size 1× 784) and its associated target (of size 1× 10). For a given value of the time-discretization parameter , we perform the first phase for T steps. Then, we perform on the one hand BPTT over K steps (to compute the gradients ∇BPTT), on the other hand C-EP (or C-VF) over K steps for given values of β and η (to compute the normalized updates ∆C−EP or ∆C−VF) and compare the gradients and normalized updates provided by the two algorithms. Precise values of the hyperparameters , T , K, β and η are given in Tab. E.6.

heading: E.1 DISCRETE-TIME RNN WITH SYMMETRIC WEIGHTS TRAINED WITH EP
text: Context of use. This model is used for training experiments in Section 4.2 and Table 4.1.
Equations with N = 2. We consider the layered architecture of Fig. 7, where s0 denotes the output layer, and the feedback connections are constrained to be the transpose of the feedforward connections, i.e. Wnn−1 = W>n−1n. In the discrete-time setting of EP, the dynamics of the first phase
are defined as:
∀t ∈ [0, T ] :  s 0 t+1 = σ ( W01 · s1t ) , s1t+1 = σ ( W12 · s2t +W>01 · s0t ) ,
s2t+1 = σ ( W2x · x+W>12 · s1t ) .
(76)
In the second phase the dynamics reads:
∀t ∈ [0,K] :  s0,βt+1 = σ ( W01 · s1,βt ) + β ( y − s0,βt ) , s1,βt+1 = σ ( W12 · s2,βt +W>01 · s 0,β t ) , s2,βt+1 = σ ( W2x · x+W>12 · s 1,β t ) .
(77)
As usual, y denotes the target. Consider the function:
Φ ( x, s2, s1, s0 ) = (s0)> ·W01 · s1 + (s1)> ·W12 · s2 + (s2)> ·W2x · x. (78)
We can compute, for example:
∂Φ ∂s1 = W12 · s2 +W>01 · s0. (79)
Comparing Eq. (76) and Eq. (79), and ignoring the activation function σ, we can see that
s1t ≈ ∂Φ ∂s1 ( x, s2t−1, s 1 t−1, s 0 t−1 ) . (80)
And similarly for the layers s0 and s2.
According to the definition of ∆C−EPθ in Eq. (19), for every layer and every t ∈ [0,K]: ∆EPW2x(β, t) = 1 β ( s2,βt+1 · x> − s 2,β t · x> ) , ∆EPW12(β, t) = 1 β ( s1,βt+1 · s 2,β> t+1 − s 1,β t · s 2,β> t ) ,
∆EPW01(β, t) = 1 β ( s0,βt+1 · s 1,β> t+1 − s 0,β t · s 1,β> t ) .
(81)
Simplifying the equations with N = 2. To go from our multi-layered architecture to the more general model presented in section 4.1. we define the state s of the network as the concatenation of all the layers’ states, i.e. s = (s2, s1, s0)> and we define the weight matrices W and Wx as:
W =  0 W>12 0W12 0 W>01 0 W01 0  , Wx = (W2x0 0 ) . (82)
Note that Eq. (76) and Eq. (78) can be vectorized into:
st+1 = σ(W · st +Wx · x), (83)
Φ = 1 2 sT ·W · s+ 1 2 sT ·Wx · x. (84)
Generalizing the equations for any N . For a general architecture with a given N , the dynamics of the first phase are defined as:
∀t ∈ [0, T ] :  s0t+1 = σ ( W01 · s1t ) snt+1 = σ ( Wnn+1 · sn+1t +W>n−1n · sn−1t ) ∀n ∈ [1, N − 1]
sNt+1 = σ ( WN,x · x+W>N−1N · s N−1 t ) ,
(85)
and those of the second phase as:
∀t ∈ [0,K] :  s0,βt+1 = σ ( W01 · s1,βt ) + β(y − s0,β) sn,βt+1 = σ ( Wnn+1 · sn+1,βt +W>n−1n · s n−1,β t ) ∀n ∈ [1, N − 1] sN,βt+1 = σ ( WN,x · x+W>N−1N · s N−1,β t ) ,
(86)
where y denotes the target. Defining:
Φ(x, sN , . . . , s0) = N−1∑ n=0 sn > ·Wnn+1 · sn+1 + sN ·WN,x · x, (87)
ignoring the activation function σ, Eq. (85) rewrites:
snt+1 ≈ ∂Φ
∂sn (x, sN , . . . , s0) ∀n ∈ [1, N − 1] (88)
According to the definition of ∆C−EPθ in Eq. (19), for every layer Wnn+1 and every t ∈ [0,K]: ∆ EP WN,x (β, t) = 1β ( sN,βt+1 · x> − s N,β t · x> ) ,
∆EPWnn+1(β, t) = 1 β ( sn,βt+1 · s n+1,β> t+1 − s n,β t · s n+1,β> t ) ∀n ∈ [0, N − 1]
(89)
Defining s = (sN , sN−1, . . . , s0)> and:
W =  0 W>N−1N 0 0 0 WN−1N 0 W > N−2N−1 0 0 0 WN−2N−1 0 . . . 0
0 0 . . . 0 W>01 0 0 0 W01 0
 , Wx =  WN,x 0 ... 0  , (90)
Eq. (85) and Eq. (87) can also be vectorized into:
st+1 = σ(W · st +Wx · x) (91)
Φ(x, s,W,Wx) = 1 2 sT ·W · s+ 1 2 sT ·Wx · x. (92)
Thereafter we introduce the other models in this general case.

heading: E.2 DISCRETE-TIME RNN WITH SYMMETRIC WEIGHTS TRAINED WITH C-EP
text: Context of use. This model is used for training experiments in Section 4.2 and Table 4.1.
Equations. Recall that we consider the layered architecture of Fig. 7, where s0 denotes the output layer. Just like in the discrete-time setting of EP, the dynamics of the first phase are defined as:
∀t ∈ [0, T ] :  s0t+1 = σ ( W01 · s1t ) snt+1 = σ ( Wnn+1 · sn+1t +W>n−1n · sn−1t ) ∀n ∈ [1, N − 1]
sNt+1 = σ ( WN,x · x+W>N−1N · s N−1 t ) (93)
Again, as in EP, the feedback connections are constrained to be the transpose of the feedforward connections, i.e. Wnn−1 = W>n−1n. In the second phase the dynamics reads:
∀t ∈ [0,K] :  s0,β,ηt+1 = σ(W01 · s 1,β,η t ) + β ( y − s0,β,ηt ) sn,β,ηt+1 = σ(Wnn+1 · s n+1,β,η t +W > n−1n · s n−1,β,η t ) ∀n ∈ [1, N − 1], sN,β,ηt+1 = σ ( WN,x · x+W>N−1N · s N−1,β,η t ) θβ,ηt+1 = θ β,η t + η∆ C−EP θ (β, η, t) ∀θ ∈ {Wnn+1}
(94) As usual, y denotes the target. Since Eq. (93) and Eq. (85) are the same, the equations describing the C-EP model can also be written in a vectorized block-wise fashion, as in Eq. (91) and Eq. (92). We can consequently define the C-EP model in Section 4.1 per Eq. (10).
According to the definitions of Eq. (6) and Eq. (7), for every layer Wnn+1 and every t ∈ [0,K]: ∆ C−EP WN,x (β, t) = 1β ( sN,β,ηt+1 · x> − s N,β,η t · x> ) ,
∆C−EPWnn+1(β, t) = 1 β ( sn,β,ηt+1 · s n+1,β,η> t+1 − s n,β,η t · s n+1,β,η> t ) ∀n ∈ [0, N − 1]
(95)

heading: E.3 REAL-TIME RNN WITH SYMMETRIC WEIGHTS TRAINED WITH C-EP
text: Context of use. This model has not been used in this work. We only introduce it for completeness with respect to Ernoult et al. (2019).
Equations. For this model, the primitive function is defined as:
Φ ( s0, s1, . . . , sN−1 ) = 1
2 (1− ) ( N∑ n=0 ||sn||2 ) + N−1∑ n=0 σ(sn)·Wnn+1·σ(sn+1)+σ(sN )·WN,x·σ(x)
(96)
so that the equations of motion read:
∀t ∈ [0, T ] :  s0t+1 = (1− )s0t + W01 · σ(s1t ) snt+1 = (1− )snt + (Wnn+1 · σ ( sn+1t ) +W>n−1n · σ(sn−1t )) ∀n ∈ [1, N − 1]
sNt+1 = (1− )sNt + (WN,x · σ (x) +W>N−1N · σ(s N−1 t ))
In the second phase:
∀t ∈ [0,K] :  s0,β,ηt+1 = (1− )s 0,β,η t + W01 · σ(s 1,β,η t ) + β (y − s0,β,η(t)) sn,β,ηt+1 = (1− )s n,β,η t + (Wnn+1 · σ(s n+1,β,η t ) +W > n−1n · σ(s n−1,β,η t )) ∀n ∈ [1, N − 1], sN,β,ηt+1 = (1− )s N,β,η t + (WN,x · σ (x) +W>N−1N · σ(s N−1,β,η t ))
θβ,ηt+1 = θ β,η t + η∆ C−EP θ (β, η, t) ∀θ ∈ {Wnn+1}
(97) where is a time-discretization parameter and y denotes the target.
According the definition of the C-EP dynamics (Eq. (6)), the definition of ∆C−EPθ (Eq. (7)) and the explicit form of Φ (Eq. 96), for all time step t ∈ [0,K], we have: ∆C−EPWN,x (β, η, t) = 1 β ( σ ( sN,β,ηt+1 ) · σ (x)> − σ ( sN,β,ηt ) · σ (x)> ) ∆C−EPWnn+1(β, η, t) = 1 β ( σ ( sn,β,ηt+1 ) · σ ( sn+1,β,ηt+1 )> − σ ( sn,β,ηt ) · σ ( sn+1,β,ηt )>) ∀n ∈ [0, N − 1]

heading: E.4 DISCRETE-TIME RNN WITH ASYMMETRIC WEIGHTS TRAINED WITH C-VF
text: Context of use. This model is used for training experiments in Section 4.3 and Table 4.1.
Equations. Recall that we consider the layered architecture of Fig. 7, where s0 denotes the output layer. The dynamics of the first phase in C-VF are defined as:
∀t ∈ [0, T ] :  s0t+1 = σ(W01 · s1t ) snt+1 = σ(Wnn+1 · sn+1t +Wnn−1 · sn−1t ) ∀n ∈ [1, N − 1] sNt+1 = σ ( WN,x · x+WNN−1 · sN−1t
) (98) Here, note the difference with EP and C-EP: the feedforward and feedback connections are unconstrained. In the second phase of C-VF:
∀t ∈ [0,K] :  s0,β,ηt+1 = σ(W01 · s 1,β,η t ) + β (y − s 0,β,η t ) sn,β,ηt+1 = σ(Wnn+1 · s n+1,β,η t +Wnn−1 · s n−1,β,η t ) ∀n ∈ [1, N − 1], sN,β,ηt+1 = σ ( WN,x · x+WNN−1 · sN−1,β,ηt
) θβ,ηt+1 = θ β,η t + η∆ C−EP θ (β, η, t) ∀θ ∈ {Wnn+1,Wn+1n}
(99) As usual y denotes the target. Note that Eq. (98) can also be in a vectorized block-wise fashion as Eq. (91) with s = (s0, s1, . . . , sN−1)> and provided that we define W and Wx as:
W =  0 WNN−1 0 0 0 WN−1N 0 WN−1N−2 0 0 0 WN−2N−1 0 . . . 0
0 0 . . . 0 W10 0 0 0 W01 0
 , Wx =  WN,x
0 ... 0  , (100) For all layers Wnn+1 and Wn+1n, and every t ∈ [0,K], we define: ∆C−VFWN,x (β, η, t) = 1 β (s N,β,η t+1 − s N,β,η t ) · x ∆C−VFWnn+1(β, η, t) = 1 β (s n,β,η t+1 − s n,β,η t ) · s n+1,β,η> t
∆C−VFWn+1n(β, η, t) = 1 β (s n+1,β,η t+1 − s n+1,β,η t ) · s
n,β,η> t

heading: E.5 REAL-TIME RNN WITH ASYMMETRIC WEIGHTS TRAINED WITH C-VF
text: Context of use. This model is used to generate Fig. 5 - see Table E.6 for precise hyperparameters.
Equations. For this model, the dynamics of the first phase are defined as:
∀t ∈ [0, T ] :  s0t+1 = (1− )s0t + W01 · σ ( s1t ) snt+1 = (1− )snt + ( Wnn+1 · σ ( sn+1t ) +Wnn−1 · σ ( sn−1t )) ∀n ∈ [1, N − 1]
sNt+1 = (1− )sNt + ( WN,x · σ (x) +WNN−1 · σ ( sN−1t )) where is the time-discretization parameter. Again, as in the discre-time version of C-VF, the feedforward and feedback connections Wnn−1 and Wn−1n are unconstrained. In the second phase, the dynamics reads:
∀t ∈ [0,K] :  s0,β,ηt+1 = (1− )s 0,β,η t + W01 · σ ( s1,β,ηt ) + β ( y − s0,β,ηt ) sn,β,ηt+1 = (1− )s n,β,η t + ( Wnn+1 · σ ( sn+1,β,ηt ) +Wnn−1 · σ ( sn−1,β,ηt )) ∀n ∈ [1, N − 1] sN,β,ηt+1 = (1− )s N,β,η t + ( WN,x · σ (x) +WNN−1 · σ ( sN−1,β,ηt )) θβ,ηt+1 = θ β,η t + η ∆ C−VF θ (β, η, t) ∀θ ∈ {Wnn+1,Wn+1n} (101) where y denotes the target, as usual. For every feedforward connection matrix Wnn+1 and every feedback connection matrix Wn+1n, and for every time step t ∈ [0,K] in the second phase, we define ∆C−VFWN,x (β, η, t) = 1 β (s N,β,η t+1 − s N,β,η t ) · σ(x) ∆C−VFWnn+1(β, η, t) = 1 β ( sn,β,ηt+1 − s n,β,η t ) · σ ( sn+1,β,ηt )> ∆C−VFWn+1n(β, η, t) = 1
β
( sn+1,β,ηt+1 − s n+1,β,η t ) · σ ( sn,β,ηt )>

heading: E.6 FIGURES FOR THE GDD EXPERIMENTS
text: In the following figures, we show the effect of using continual updates with a finite learning rate in terms of the ∆C−EP and −∇BPTT processes on different models introduced above. These figures have been realized either in the discrete-time or continuous-time setting with the fully connected layered architecture with one hidden layer on MNIST. Dashed an continuous lines respectively represent the normalized updates ∆ and the gradients ∇BPTT. Each randomly selected synapse or neuron correspond to one color. We add an s or θ index to specify whether we analyse neuron or synapse updates and gradients. Each C-VF simulation has been realized with an angle between forward and backward weights of 0 degrees (i.e. Ψ(θf , θb) = 0◦). For each figure, left panels demonstrate the GDD property with C-EP with η = 0 and the right panels show that, upon using η > 0, dashed and continuous lines start to split appart.
θ (t) normalized updates (η > 0) and −∇BPTTθ (t)
gradients.
θ (t) normalized updates (η > 0) and −∇BPTTθ (t) gradients.
θ (t) normalized updates (η > 0) and −∇BPTTθ (t)
gradients.
θ (t) normalized updates (η > 0) and −∇BPTTθ (t)
gradients.

heading: F EXPERIMENTAL DETAILS

heading: F.1 TRAINING EXPERIMENTS (TABLE 4.1)
text: Simulation framework. Simulations have been carried out in Pytorch. The code has been attached to the supplementary materials upon submitting this work on OpenReview. We have also attached a readme.txt with a specification of all dependencies, packages, descriptions of the python files as well as the commands to reproduce all the results presented in this paper.
Data set. Training experiments were carried out on the MNIST data set. Training set and test set include 60000 and 10000 samples respectively.
Optimization. Optimization was performed using stochastic gradient descent with mini-batches of size 20. For each simulation, weights were Glorot-initialized. No regularization technique was used and we did not use the persistent trick of caching and reusing converged states for each data sample between epochs as in Scellier and Bengio (2017).
Activation function. For training, we used the activation function
σ(x) = 1
1 + exp(−4(x− 1/2)) . (102)
Although it is a shifted and rescaled sigmoid function, we shall refer to this activation function as ‘sigmoid’.
Use of a randomized β. The option ’Random β’ appearing in the detailed table of results (Table 3) refers to the following procedure. During training, instead of using the same β accross mini-batches, we only keep the same absolute value of β and sample its sign from a Bernoulli distribution of probability 12 at each mini-batch iteration. This procedure was hinted at by Scellier and Bengio (2017) to improve test error, and is used in our context to improve the model convergence for Continual Equilibrium Propagation - appearing as C-EP and C-VF in Table 4.1 - training simulations.
Tuning the angle between forward and backward weights. In Table 4.1, we investigate C-VF initialized with different angles between the forward and backward weights - denoted as Ψ in Table 4.1. Denoting them respectively θf and θb, the angle κ between them is defined here as:
κ(θf , θb) = cos −1  Tr(θf · θ>b )√ Tr(θf · θ>f ) √ Tr(θb · θ>b )  , where Tr denotes the trace, i.e. Tr(A) = ∑ iAii for any squared matrix A. To tune arbitrarily well enough κ(θf , θb), the procedure is the following: starting from θb = θf , i.e. κ(θf , θb) = 0, we can gradually increase the angle between θf and θb by flipping the sign of an arbitrary proportion of components of θb. The more components have their sign flipped, the larger is the angle. More formally, we write θb in the form θb = M(p) θf and we define:
Ψ(p) = κ(θf ,M(p) θf), (103)
where M(p) is a mask of binary random values {+1, -1} of the same dimension of θf : M(p) = −1 with probability p and M(p) = +1 with probability 1− p. Taking the cosine and the expectation of Eq. (103), we obtain:
〈cos(Ψ(p))〉 = p×−Tr(θf · θ > f ) Tr(θf · θ>f ) + (1− p)× Tr(θf · θ > f )
Tr(θf · θ>f ) = 1− 2p
Thus, the angle Ψ between θf and θf M(p) can be tuned by the choice of p through:
p(Ψ) = 1
2 (1− 〈cos(Ψ)〉) (104)
Hyperparameter search for EP. We distinguish between two kinds of hyperparameters: the recurrent hyperparameters - i.e. T , K and β - and the learning rates. A first guess of the recurrent hyperparameters T and β is found by plotting the ∆C−EP and ∇BPTT processes associated to synapses and neurons to see qualitatively whether the theorem is approximately satisfied, and by conjointly computing the proportions of synapses whose ∆C−EPW processes have the same sign as its ∇BPTTW processes. K can also be found out of the plots as the number of steps which are required for the gradients to converge. Morever, plotting these processes reveal that gradients are vanishing when going away from the output layer, i.e. they lose up to 10−1 in magnitude when going from a layer to the previous (i.e. upstream) layer. We subsequently initialized the learning rates with increasing values going from the output layer to upstreams layers. The typical range of learning rates is [10−3, 10−1], [10, 1000] for T , [2, 100] for K and [0.01, 1] for β. Hyperparameters where adjusted until having a train error the closest to zero. Finally, in order to obtain minimal recurrent hyperparameters - i.e. smallest T and K possible - we progressively decreased T and K until the train error increases again.
Table 3: Training results on MNIST with EP, C-EP and C-VF. "#h" stands for the number of hidden layers. We indicate over five trials the mean and standard deviation for the test error, the mean error in parenthesis for the train error. T (resp. K) is the number of iterations in the first (resp. second) phase.
Full table of results. Since Table 4.1 does not show C-VF simulation results for all initial weight angles, we provide below the full table of results, including those which were used to plot Fig. 5.
Initial Ψ(θf , θb) (◦) Error (%) T K Random β Epochs
Test Train
EP-1h − 2.00± 0.13 (0.20) 30 10 No 30 EP-2h − 1.95± 0.10 (0.14) 100 20 No 50 C-EP-1h − 2.85± 0.18 (0.83) 40 15 No 100 C-EP-1h − 2.28± 0.16 (0.41) 40 15 Yes 100 C-EP-2h − 2.44± 0.14 (0.31) 100 20 No 150 C-VF-1h 0 2.43± 0.08 (0.77) 40 15 Yes 100
22.5 2.38± 0.15 (0.74) 40 15 Yes 100 45 2.37± 0.06 (0.78) 40 15 Yes 100 67.5 2.48± 0.15 (0.81) 40 15 Yes 100 90 2.46± 0.18 (0.78) 40 15 Yes 100
112.5 4.51± 3.96 (2.92) 40 15 Yes 100 135 86.61± 4.27 (88.51) 40 15 Yes 100 157.5 91.08± 0.01 (90.98) 40 15 Yes 100 180 92.82± 3.47 (92.71) 40 15 Yes 100
C-VF-2h 0 2.97± 0.19 (1.58) 100 20 Yes 150 22.5 3.54± 0.75 (2.70) 100 20 Yes 150 45 3.78± 0.78 (2.86) 100 20 Yes 150
67.5 4.59± 0.92 (4.68) 100 20 Yes 150 90 5.05± 1.17 (4.81) 100 20 Yes 150
112.5 20.33± 13.03 (20.30) 100 20 Yes 150 135 59.04± 17.97 (60.53) 100 20 Yes 150 157.5 77.90± 13.49 (78.04) 100 20 Yes 150 180 74.17± 12.76 (74.05) 100 20 Yes 150

heading: F.2 WHY C-EP DOES NOT PERFORM AS WELL AS STANDARD EP ?
text: We provide here further ground for the training performance degradation observed on the MNIST task when implementing C-EP compared to standard EP. In practice, when training with C-EP, we have to make a trade-off between:
1. having a learning rate that is small enough so that C-EP normalized updates are subsequently close enough to the gradients of BPTT (Theorem 1),
2. having a learning rate that is large enough to ensure convergence within a reasonable number of epochs.
In other words, the degradation of accuracy observed in the table of Fig. 4.1 is due to using a learning rate that is too large to observe convergence within 100 epochs. To demonstrate this, we implement Alg. 5 which simply consists in using a very small learning rate throughout the second phase (denoted as ηtiny), and artificially rescaling the resulting weight update by a bigger learning rate (denoted as η). Applying Alg. 5 to a fully connected layered architecture with one hidden layer, T = 30, K = 10, β = 0.1, yields 2.06± 0.13% test error and 0.18± 0.01% train error over 5 trials, where we indicate mean and standard deviation. Similarly, applying Alg. 5 to a fully connected layered architecture with two hidden layers, T = 100, K = 20, β = 0.5, yields 1.89 ± 0.22% test error and 0.02 ± 0.02% train error. These results are exactly the same as the one provided by standard EP - see Table 3.
Algorithm 5 Debugging procedure of C-EP Input: x, y, θ, β, η, ηtiny = 10−5η. Output: θ.
1: s0 ← 0 . First Phase 2: ∆θ ← 0 . Temporary variable accumulating parameter updates 3: repeat 4: st+1 ← ∂Φ∂s (x, st, θ) 5: until st = s∗ 6: sβ0 ← s∗ . Second Phase 7: repeat 8: sβt+1 ← ∂Φ∂s ( x, sβt , θ ) − β ∂`∂s ( sβt , y
) 9:
10: θ ← θ + ηtinyβ ( ∂Φ ∂θ ( sβt+1 ) − ∂Φ∂θ ( sβt )) 11: ∆θ ← ∆θ + ηtinyβ ( ∂Φ ∂θ ( sβt+1 ) − ∂Φ∂θ ( sβt
)) 12: until sβt and θ are converged. 13: θ ← θ −∆θ + ηηtiny ∆θ . Rescale the total parameter update by η ηtiny

Review summaries:
The paper is concerned with biologically plausible models of learning. It takes equilibrium propagation -- where updates depend on local spatial information, in the sense that the information is available at each neuron -- and modifies the algorithm so updates are also local in time, thus obtaining C-EP. The key insight is that the updates in EP can be written as a telescoping sum over time points, eq (5).