Title: EQUILIBRIUM PROPAGATION WITH CONTINUAL WEIGHT UPDATES

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.
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.
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 θ

1. INTRODUCTION:

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).

and evolving according to the dynamics:

EQUATION

F is the transition function of the system. Assuming convergence of the dynamics before time step T , we have s T = 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 s T = s * for the current value of θ. The goal is to optimize the parameter θ in order to minimize a loss:

EQUATION

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 s 1 , s 2 , s 3 . . . 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,

EQUATION

3. EQUILIBRIUM PROPAGATION WITH CONTINUAL WEIGHT UPDATES (C-EP):

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.

3.1. FROM EP TO C-EP: AN INTUITION BEHIND CONTINUAL WEIGHT UPDATES:

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:

EQUATION

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. -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.

3.2. DESCRIPTION OF THE C-EP ALGORITHM:

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 s t . The dynamics of s t 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:

EQUATION

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.

3.3. GRADIENT-DESCENDING DYNAMICS (GDD):

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 = (s T , y) after T time steps, computed with BPTT: 

EQUATION

More details about L and the gradients ∇ BPTT s (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:

EQUATION

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 s 0 , s 1 , . . . , s T be the convergent sequence of states and denote s * = s T the steady state. Further assume that there exists some step K where 0 < K ≤ T such that s * = s T = s T -1 = . . . s T -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) = -∇ BPTT s (t), lim β→0 lim η→0 ∆ C-EP θ (β, η, t) = -∇ BPTT θ (t). (9) Theorem 1 rewrites s β,η t+1 ≈ s β,η t -β ∂L ∂s T -t and θ β,η t+1 ≈ θ β,η t -η ∂L ∂θ T -t ,

showing that in the second p ĥ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. 2 (s t + θ), the steady state is denoted s * and the loss is L * = 1 2 s 2 * . For completeness, we also include the corresponding gradients of Recurrent Backpropagation (RBP) and the normalized updates of EP, denoted ∇ RBP and ∆ EP respectively. The equivalence between C-EP, EP, RBP and BPTT holds in the general setting: see Appendix A for a thorough study of their relationship.

4. NUMERICAL EXPERIMENTS:

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) .

4.1. MODELS: VANILLA RNNS WITH SYMMETRIC AND ASYMMETRIC WEIGHTS:

Vanilla RNN with symmetric weights trained by C-EP. The first phase dynamics is defined as:

EQUATION

)

where σ is an activation function, W is a symmetric weight matrix connecting the layers s and W x 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 s t+1 ≈ ∂Φ ∂s (s t , W ) with Φ(s, W ) = 1 2 s • W • s if we ignore the activation function σ. Following Eq. ( 6) and Eq. ( 7), we define the normalized updates of this model as:

EQUATION

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 s 0 , s 1 , ..., s N , with W n,n+1 connecting the layers s n+1 and s n in both directions, the corresponding primitive function is

Φ = n (s n ) • W n,n+1 • s n+1 + s N • W x • 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:

EQUATION

)

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 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.

4.3. CONTINUAL VECTOR FIELD (C-VF) TRAINING EXPERIMENTS:

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.

5. DISCUSSION:

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) . 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.

Appendix:



A PROOF OF THEOREM 1:

In this appendix, we prove Theorem 1, which we recall here. Theorem 1 (GDD Property). Let s 0 , s 1 , . . . , s T be the convergent sequence of states and denote s * = s T the steady state. Further assume that there exists some step K where 0 < K ≤ T such that s * = s T = s T -1 = . . . s T -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 : 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).

EQUATION

A.2 SKETCH OF THE PROOF:

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:

EQUATION

)

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:

EQUATION

)

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.

s T -K = s T -K+1 = • • • = s T -1 = s T = s * .

Then the first K gradients of BPTT are equal to the first K gradient of RBP, i.e. ∀t = 0, 1, . . . , K :

EQUATION

)

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) . 

EQUATION

)

We have also defined the normalized updates of C-EP:

EQUATION

)

We also recall the dynamics of EP in the second phase:

EQUATION

as well as the normalized updates of EP, as defined in Ernoult et al. (2019) :

EQUATION

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:

EQUATION

)

Proof of Lemma 2. We want to compute the limits of ∆ C-EP s (β, η, 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,

EQUATION

)

Now, taking η = 0 in the bottom equation of Eq. ( 16) yields the recurrence relation θ β,0 t+1 = θ β,0 t , so that θ β,0 t = θ β,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

EQUATION

It follows from Eq. ( 20) and Eq. ( 21) that

EQUATION

)

Now let us compute lim η→0 (η>0) ∆ C-EP s (β, η, t). Using Eq. ( 16), we have

EQUATION

)

Similarly as before, for fixed t, ∂Φ ∂θ x, s β,η t , θ β,η t is a continuous function of η. Therefore

EQUATION

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, s t+1 , θ 0 ) - ∂Φ ∂θ (x, s t+1 , θ 0 ) (definition of ∆ EP θ , Eq. (19)) = η β ∂Φ ∂θ (x, s β K , θ 0 ) - ∂Φ ∂θ (x, s * , θ 0 )

B EQUIVALENCE OF BPTT AND RBP:

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.

B.1 SKETCH OF THE PROOF OF LEMMA 4:

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.

s T -K = s T -K+1 = • • • = s T -1 = s T = s * .

Then the first K gradients of BPTT are equal to the first K gradient of RBP, i.e. ∀t = 0, 1, . . . , K :

EQUATION

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 = (s T , y),

• RBP computes the gradients of the loss at the steady state, i.e. L * = (s * , y).

B.2 BACKPROPAGATION THROUGH TIME (BPTT):

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 s T after T time steps, denoted L = (s T , 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 s t+1 = F (x, s t , θ 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:

EQUATION

)

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:

EQUATION

)

Proposition 5 (Gradients of BPTT). The gradients ∇ BPTT s (t) and ∇ BPTT θ (t) satisfy the recurrence relationship

EQUATION

EQUATION

EQUATION

B.3 FROM BACKPROP THROUGH TIME (BPTT) TO RECURRENT BACKPROP (RBP):

In general, to apply BPTT, it is necessary to store in memory the history of past hidden states s 1 , s 2 , . . . , s T in order to compute the gradients ∇ BPTT s (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 s T -K = s T -K+1 = • • • = s T -1 = s T = 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 s T , s T -1 , . . . , s T -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 ∇ RBP s (t) and ∇ RBP θ (t), which we call 'gradients of RBP', iteratively for t = 0, 1, 2, . . . Definition 6 (Gradients of RBP). The gradients ∇ RBP s (t) and ∇ RBP θ (t) are defined and computed iteratively as follows:

EQUATION

EQUATION

EQUATION

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: s 0 ← 0 2: for t = 0 to T -1 do 3: s t+1 ← F (x, s t , θ) 4: end for 5: ∇ BPTT s (0) ← ∂ ∂s (s T , y) 6: for t = 1 to T do 7: ∇ s (t) ← ∂F ∂s (x, s T -t , θ) • ∇ s (t -1) 8: ∇ θ (t) ← ∂F ∂θ (x, s T -t , θ) • ∇ s (t -1) 9: end for 10: ∇ BPTT θ (tot) ← T -1 t=0 ∇ BPTT θ (t)

Algorithm 4 RBP Input: x, y, θ. Output: θ.

1: s 0 ← 0 2: repeat 3:

s t+1 ← F (x, s t , θ) 4: until s t = s * 5: ∇ RBP s (0) ← ∂ ∂s (s * , y) 6: repeat 7: ∇ RBP s (t) ← ∂F ∂s (x, s * , θ) • ∇ RBP s (t -1) 8: ∇ RBP θ (t) ← ∂F ∂θ (x, s * , θ) • ∇ RBP s (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 s T -t is required to compute ∂F ∂s (x, s T -t , θ) and ∂F ∂θ (x, s T -t , θ) ; thus it is necessary to store in memory the sequence of states s 1 , s 2 , . . . , s T . 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.

B.4 WHAT 'GRADIENTS' ARE THE GRADIENTS OF RBP?:

In this subsection we motivate the name of 'gradients' for the quantities ∇ RBP s (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.

EQUATION

∇ RBP s (t) and ∇ RBP θ (t) can also be expressed as gradients of L t = (s t , y), the cost after t time steps. In the terminology of Scellier and Bengio (2019) , L t was named the projected cost. For t = 0, L 0 is simply the cost of the initial state s 0 . For t > 0, L t is the cost of the state projected a duration t in the future. 

EQUATION

Rearranging the terms, and using the Taylor expansion

EQUATION

EQUATION

)

Proof of Proposition 8. By the chain rule of differentiation we have

EQUATION

)

Evaluation this expression for s 0 = s * we get 

EQUATION

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 θ:

EQUATION

where s * is the steady state of the dynamics (it is easy to see that the solution is s * = θ). The dynamics rewrites s t+1 = F (s t , θ) with the transition function F (s, θ) = 1 2 (s + θ), and the loss rewrites L * = (s * ) with the cost function (s) = 1 2 s 2 . Furthermore, a primitive function of the system 1 is Φ(s, θ) = 1 4 (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 s T , and we approximate L * (the loss at the steady state) by the loss after T time steps L = (s T ).

In order to compute (i.e. 'backpropagate') the gradients of BPTT, Proposition 5 tells us that we need to compute ∂ ∂s (s T ) = s T , ∂F ∂s (s t , θ) = 1 2 and ∂F ∂θ (s t , θ) = 1 2 . We get

EQUATION

)

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

EQUATION

)

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:

EQUATION

)

This linear dynamical system can be solved analytically:

EQUATION

)

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

EQUATION

)

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:

EQUATION

)

1 The primitive function Φ is determined up to a constant.

First, rearranging the terms in the second equation, we get

EQUATION

)

It follows that

EQUATION

)

Therefore, all we need to do is to compute ∆ C-EP s (β, η, t). Second, by iterating the second equation over all indices from t = 0 to t -1 we get

EQUATION

)

Using s * = θ and plugging this into the first equation we get

EQUATION

)

Solving this linear dynamical system, and using the initial condition s β,η 0 = θ we get

EQUATION

)

Finally:

EQUATION

D COMPLEMENT ON GRADIENT-DESCENDING DYNAMICS (GDD)

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.

D.1 GRADIENT-DESCENDING UPDATES (GDU) OF ERNOULT ET AL. ( 2019)

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.

D.2 A GENERALISATION OF THE GDD PROPERTY:

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:

EQUATION

)

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:

EQUATION

Now let us define the transition function F (s, W ) = σ(W • s), so that the dynamics of the first phase rewrites

EQUATION

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 :

EQUATION

Now, recall the definition of the normalized updates of C-VF, as well as the gradients of the loss L = (s T , y) after T time steps, computed with BPTT:

EQUATION

The loss L and the gradients ∇ BPTT s (t) and ∇ BPTT θ (t) are defined formally in Appendix B.2.

Theorem 9 (Generalisation of the GDD Property). Let s 0 , s 1 , . . . , s T be the convergent sequence of states and denote s * = s T the steady state. Further assume that there exists some step K where 0 < K ≤ T such that s * = s T = s T -1 = . . . s T -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 :

EQUATION

A few remarks need to be made:

EQUATION

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

EQUATION

In this case the extra assumption in Theorem 9 is automatically satisfied.

D.3 PROOF OF THEOREM 9:

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 ∇ BPTT s (t) and ∇ BPTT θ (t) satisfy the recurrence relationship

EQUATION

EQUATION

EQUATION

Lemma 10 (Updates of C-VF). Define the (normalized) neural and weight updates of C-VF in the limit η → 0 and β → 0:

EQUATION

They satisfy the recurrence relationship

EQUATION

EQUATION

EQUATION

The proof of Lemma 10 is similar to the one provided in Ernoult et al. (2019) .

are defined as:

∀t ∈ [0, T ] :

EQUATION

In the second phase the dynamics reads:

EQUATION

As usual, y denotes the target. Consider the function:

EQUATION

We can compute, for example:

EQUATION

)

Comparing Eq. ( 76) and Eq. ( 79), and ignoring the activation function σ, we can see that

EQUATION

And similarly for the layers s 0 and s 2 .

According to the definition of ∆ C-EP θ in Eq. ( 19), for every layer and every t ∈ [0, K]:

EQUATION

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 = (s 2 , s 1 , s 0 ) and we define the weight matrices W and W x as:

EQUATION

Note that Eq. ( 76) and Eq. ( 78) can be vectorized into:

EQUATION

EQUATION

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 ] :

EQUATION

and those of the second phase as:

EQUATION

where y denotes the target. Defining:

EQUATION

ignoring the activation function σ, Eq. ( 85) rewrites:

EQUATION

According to the definition of ∆ C-EP θ in Eq. ( 19), for every layer W nn+1 and every t ∈ [0, K]:

EQUATION

Defining s = (s N , s N -1 , . . . , s 0 ) and:

EQUATION

Eq. ( 85) and Eq. ( 87) can also be vectorized into:

EQUATION

)

Thereafter we introduce the other models in this general case. Equations. Recall that we consider the layered architecture of Fig. 7 , where s 0 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 ] :

   s 0 t+1 = σ W 01 • s 1 t s n t+1 = σ W nn+1 • s n+1 t + W n-1n • s n-1 t ∀n ∈ [1, N -1] s N t+1 = σ W N,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. W nn-1 = W n-1n . In the second phase the dynamics reads:

EQUATION

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 W nn+1 and every t ∈ [0, K]:

EQUATION

E.3 REAL-TIME RNN WITH SYMMETRIC WEIGHTS TRAINED WITH C-EP:

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:

EQUATION

so that the equations of motion read:

∀t ∈ [0, T ] :

   s 0 t+1 = (1 -)s 0 t + W 01 • σ(s 1 t ) s n t+1 = (1 -)s n t + (W nn+1 • σ s n+1 t + W n-1n • σ(s n-1 t )) ∀n ∈ [1, N -1] s N t+1 = (1 -)s N t + (W N,x • σ (x) + W N -1N • σ(s N -1 t ))

In the second phase:

∀t ∈ [0, K] :

EQUATION

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: Equations. Recall that we consider the layered architecture of Fig. 7 , where s 0 denotes the output layer. The dynamics of the first phase in C-VF are defined as:

     ∆ C-EP W N,x (β, η, t) = 1 β σ s N,β,η t+1 • σ (x) -σ s N,β,η t • σ (x) ∆ C-EP Wnn+1 (β, η, t) = 1 β σ s n,β,η t+1 • σ s n+1,β,η t+1 -σ s n,β,η t • σ s n+1,β,η t ∀n ∈ [0, N -1] E.4

∀t ∈ [0, T ] :

   s 0 t+1 = σ(W 01 • s 1 t ) s n t+1 = σ(W nn+1 • s n+1 t + W nn-1 • s n-1 t ) ∀n ∈ [1, N -1] s N t+1 = σ W N,x • x + W N N -1 • s N -1 t (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] :

EQUATION

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 = (s 0 , s 1 , . . . , s N -1 ) and provided that we define W and W x as:

EQUATION

For all layers W nn+1 and W n+1n , and every t ∈ [0, K], we define: E .6 for precise hyperparameters.

       ∆ C-VF W N,x (β, η, t) = 1 β (s N,β,η t+1 -s N,β,η t ) • x ∆ C-VF Wnn+1 (β, η, t) = 1 β (s n,β,η t+1 -s n,β,η t ) • s n+1,β,η t ∆ C-VF Wn+1n (β, η, t) = 1 β (s n+1,β,η t+1 -s n+1,β,η t ) • s n,β

Equations. For this model, the dynamics of the first phase are defined as:

∀t ∈ [0, T ] :

   s 0 t+1 = (1 -)s 0 t + W 01 • σ s 1 t s n t+1 = (1 -)s n t + W nn+1 • σ s n+1 t + W nn-1 • σ s n-1 t ∀n ∈ [1, N -1] s N t+1 = (1 -)s N t + W N,x • σ (x) + W N N -1 • σ s N -1 t

where is the time-discretization parameter. Again, as in the discre-time version of C-VF, the feedforward and feedback connections W nn-1 and W n-1n are unconstrained. In the second phase, the dynamics reads:

EQUATION

101) where y denotes the target, as usual. For every feedforward connection matrix W nn+1 and every feedback connection matrix W n+1n , and for every time step t ∈ [0, K] in the second phase, we define

           ∆ C-VF W N,x (β, η, t) = 1 β (s N,β,η t+1 -s N,β,η t ) • σ(x) ∆ C-VF Wnn+1 (β, η, t) = 1 β s n,β,η t+1 -s n,β,η t • σ s n+1,β,η t ∆ C-VF Wn+1n (β, η, t) = 1 β s n+1,β,η t+1 -s n+1,β,η t • σ s n,β,η t E.6 FIGURES FOR THE GDD EXPERIMENTS

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. F EXPERIMENTAL DETAILS F.1 TRAINING EXPERIMENTS (TABLE 4.1) 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. 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

EQUATION

)

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 1 2 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) = i A ii 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:

EQUATION

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: 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:

cos(Ψ(p)) = p × - Tr(θ f • θ f ) Tr(θ f • θ f ) + (1 -p) × Tr(θ f • θ f ) Tr(θ f • θ f ) = 1 -2p

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: θ. 

:



Figures and Tables:

Figure 2: Gradient-Descending Dynamics (GDD, Theorem 1). In the second phase of Continual Equilibrium Prop (C-EP), the dynamics of neurons and synapses descend the gradients of BPTT, i.e. ∆ C-EP (t) = -∇ BPTT (t). The colors illustrate when corresponding computations are realized in C-EP and BPTT. Top left. 1 st phase of C-EP with static input x and target y. The final state s T is the steady state s * . Bottom left. Backprop through time (BPTT). Bottom right. 2 nd phase of C-EP. The starting state s β,η 0 is the final state of the forward-time pass, i.e. the steady state s * .

Figure 3: Illustration of Theorem 1 on a simple model. The variables s and θ are scalars, the first phase equation is s t+1 = 12 (s t + θ), the steady state is denoted s * and the loss isL * = 1 2 s 2 * .For completeness, we also include the corresponding gradients of Recurrent Backpropagation (RBP) and the normalized updates of EP, denoted ∇ RBP and ∆ EP respectively. The equivalence between C-EP, EP, RBP and BPTT holds in the general setting: see Appendix A for a thorough study of their relationship.

Figure 4: Left: Training results on MNIST with EP, C-EP and C-VF. "#h" stands for the number of hidden layers. We indicate over 5 trials the mean and standard deviation for the test error (mean train error in parenthesis). T (resp. K) is the number of iterations in the 1 st (resp. 2 nd ) phase. For C-VF results, the initial angle between forward (θ f ) and backward (θ b ) weights is Ψ(θ f , θ b ) = 0 • . Right: Test error rate on MNIST achieved by C-VF as a function of the initial Ψ(θ f , θ b ).

Figure 5: Three versions of EP: standard Equilibrium Propagation (EP), Continual Equilibrium Propagation (C-EP) and Continual Vector Field EP (C-VF). #-h denotes the number of hidden layers. (a): test error rate on MNIST as a function of the initial angle Ψ between the total normalized update of EP and the total gradient of BPTT. (b): Dashed and continuous lines respectively represent the normalized updates ∆ θ (t) (i.e. ∆ EP θ (t), ∆ C-EP θ

EQUIVALENCE OF C-EP AND EP First, recall the dynamics of C-EP in the second phase: starting from s β,η 0 = s * and θ β,η 0 = θ we have ∀t ≥ 0 :

DISCRETE-TIME RNN WITH SYMMETRIC WEIGHTS TRAINED WITH C-EP Context of use. This model is used for training experiments in Section 4.2 and Table 4.1.

-TIME RNN WITH ASYMMETRIC WEIGHTS TRAINED WITH C-VF Context of use. This model is used to generate Fig. 5 -see Table

Figure 8: Discrete-Time RNN with symmetric weights. Left: ∆ C-EP s (t) normalized updates (η = 0) and -∇ BPTTs

Figure 9: Discrete-Time RNN with symmetric weights. Left: ∆ C-EP θ (t) normalized updates (η = 0) and -∇ BPTTθ

Figure 10: Real-Time RNN with symmetric weights. Left: ∆ C-EP s (t) normalized updates (η = 0) and -∇ BPTTs

Figure 11: Real-Time RNN with symmetric weights. Left∆ C-EP θ (t) normalized updates (η = 0) and -∇ BPTTθ

Figure 12: Discrete-Time RNN with asymmetric weights. Left: ∆ C-VF s (t) normalized updates (η = 0) and -∇ BPTTs

Figure 13: Discrete-Time RNN with asymmetric weights. Left: ∆ C-VF θ (t) normalized updates (η = 0) and -∇ BPTTθ

Figure 14: Real-Time RNN with asymmetric weights. Left: ∆ C-VF s (t) normalized updates (η = 0) and -∇ BPTTs

Figure 15: Real-Time RNN with asymmetric weights. Left: ∆ C-VF θ (t) normalized updates (η = 0) and -∇ BPTTθ

Data set. Training experiments were carried out on the MNIST data set. Training set and test set include 60000 and 10000 samples respectively.

Figure 16: Train and test error achieved on MNIST with Continual Equilibrium Propagation (C-EP) on the Discrete-Time RNN model with symmetric weights. Plain lines indicate mean, shaded zones delimiting mean plus/minus standard deviation over 5 trials. Left: C-EP on the fully connected layered architecture with one hidden layer (784-512-10) without beta randomization. Middle: C-EP on the fully connected layered architecture with one hidden layer (784-512-10) with beta randomization. Right: C-EP on the fully connected layered architecture with two hidden layers (784-512-512-10) without beta randomization.

Figure 17: Train and test error achieved on MNIST by Continual Vector Field Equilibrium Propagation (C-VF) on the Discrete-Time RNN model with asymmetric weights with one hidden layer (784-512-10) for different initialization for the angle between forward and backward weights (Ψ). Plain lines indicate mean, shaded zones delimiting mean plus/minus standard deviation over 5 trials.

Figure 18: Train and test error achieved on MNIST by Continual Vector Field Equilibrium Propagation (C-VF) on the vanilla RNN model with asymmetric weights with two hidden layers (784-512-512-10) for different initialization for the angle between forward and backward weights (Ψ). Plain lines indicate mean, shaded zones delimiting mean plus/minus standard deviation over 5 trials.

Proposition 8 (Gradients of RBP are Gradients of the Projected Cost). The 'RBP gradients' ∇ RBPwhere the initial state s 0 is the steady state s * .

Table content: <table><tr><td colspan="13">Proof of Proposition 7. First of all, by Definition 6 (Eq. 32-34) it is straightforward to see that</td></tr><tr><td colspan="2">∀t ≥ 0,</td><td colspan="2">∇ RBP s</td><td colspan="4">(t) =</td><td/><td colspan="2">∂F ∂s</td><td colspan="2">(x, s  *  , θ)</td><td>t</td><td>•</td><td>∂ ∂s</td><td>(s  *  , y) ,</td><td>(37)</td></tr><tr><td colspan="2">∀t ≥ 1,</td><td colspan="2">∇ RBP θ</td><td colspan="4">(t) =</td><td colspan="2">∂F ∂θ</td><td colspan="3">(x, s  *  , θ) •</td><td>∂F ∂s</td><td>(x, s  *  , θ)</td><td>t-1</td><td>•</td><td>∂ ∂s</td><td>(s  *  , y) .</td><td>(38)</td></tr><tr><td colspan="6">Second, recall that the loss L  *  is</td><td/><td/><td/><td/><td/><td colspan="2">L  *  = (s  *  , y) ,</td><td>(39)</td></tr><tr><td>where</td><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td/><td/><td colspan="3">s  *  = F (x, s  *  , θ) .</td><td>(40)</td></tr><tr><td colspan="13">By the chain rule of differentiation, the gradient of L  *  (Eq. 39) is</td></tr><tr><td/><td/><td/><td/><td/><td/><td/><td/><td colspan="3">∂L  *  ∂θ</td><td>=</td><td>∂ ∂s</td><td>(s  *  , y) •</td><td>∂s  *  ∂θ</td><td>.</td><td>(41)</td></tr><tr><td colspan="13">In order to compute ∂s *  ∂θ , we differentiate the steady state condition (Eq. 40) with respect to θ, which yields</td></tr><tr><td/><td/><td/><td colspan="2">∂s  *  ∂θ</td><td>=</td><td colspan="3">∂F ∂s</td><td colspan="4">(x, s  *  , θ) •</td><td>∂s  *  ∂θ</td><td>+</td><td>∂F ∂θ</td><td>(x, s</td></tr><tr><td>and ∇ RBP θ</td><td colspan="12">(t) can be expressed as gradients of the projected cost:</td><td>s</td><td>(t)</td></tr><tr><td/><td colspan="2">∀t ≥ 0,</td><td colspan="4">∇ RBP s</td><td colspan="4">(t) =</td><td colspan="2">∂L t ∂s 0 s0=s *</td><td>,</td><td>∇ RBP θ</td><td>(t) =</td><td>∂L t ∂θ 0 s0=s *</td><td>(36)</td></tr></table>

DISCRETE-TIME RNN WITH ASYMMETRIC WEIGHTS TRAINED WITH C-VF Context of use. This model is used for training experiments in Section 4.3 and Table 4.1.

Table content: <table/>

Table of hyperparameters used to demonstrate Theorem 1.

Table content: <table><tr><td/><td colspan="3">Figure Angle Ψ ( • ) Activation</td><td>T</td><td>K</td><td>β</td><td>Learning rates</td></tr><tr><td>C-EP</td><td>5</td><td>0</td><td>tanh</td><td colspan="4">800 80 0.01 0.08</td><td>0 -0</td></tr><tr><td>C-VF</td><td>5</td><td>45</td><td>tanh</td><td colspan="4">800 80 0.01 0.08</td><td>0 -0</td></tr><tr><td>C-EP</td><td>5</td><td>0</td><td>tanh</td><td colspan="4">800 80 0.01 0.08 1.510 -5 -1.510 -5</td></tr><tr><td>C-VF</td><td>5</td><td>45</td><td>tanh</td><td colspan="4">800 80 0.01 0.08 1.510 -5 -1.510 -5</td></tr><tr><td colspan="2">C-VF 14-15</td><td>0</td><td>tanh</td><td colspan="4">800 80 0.005 0.08</td><td>0 -0</td></tr><tr><td colspan="2">C-VF 14-15</td><td>0</td><td>tanh</td><td colspan="4">800 80 0.005 0.08</td><td>2.10 -5 -2.10 -5</td></tr><tr><td colspan="2">C-VF 12-13</td><td>0</td><td>tanh</td><td colspan="3">150 10 0.01</td><td>-</td><td>0 -0</td></tr><tr><td colspan="2">C-VF 12-13</td><td>0</td><td>tanh</td><td colspan="3">150 10 0.01</td><td>-</td><td>2.10 -5 -2.10 -5</td></tr><tr><td colspan="2">C-EP 10-11</td><td>0</td><td>tanh</td><td colspan="4">800 80 0.05 0.08</td><td>0 -0</td></tr><tr><td colspan="2">C-EP 10-11</td><td>0</td><td>tanh</td><td colspan="4">800 80 0.05 0.08</td><td>2.10 -5 -2.10 -5</td></tr><tr><td>C-EP</td><td>8-9</td><td>0</td><td>tanh</td><td colspan="3">150 10 0.01</td><td>-</td><td>0 -0</td></tr><tr><td>C-EP</td><td>8-9</td><td>0</td><td>tanh</td><td colspan="3">150 10 0.01</td><td>-</td><td>2.10 -5 -2.10 -5</td></tr></table>

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 Table4.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 ) ( • ) WHY C-EP DOES NOT PERFORM AS WELL AS STANDARD EP ?

Table content: <table><tr><td>Error (%)</td><td>T</td><td>K Random β Epochs</td></tr><tr><td>Test</td><td>Train</td><td/></tr></table>

Peer review (iclr2020):
I think it is an intriguing paper, but unfortunately left me a bit confused. I have to admit is not a topic I*m really versed in, so it might be that this affected my evaluation of the work. But also, as a paper submitted to ICLR I would expect the paper to be self-contained and be able to provide all the details needed. I do appreciate the authors providing in the appendix the proofs of the other theorems even if they come form other works. The paper introduces C-EP, an extension of a previously introduced algorithm EP, such that it becomes biologically plausible. In particular EP is local in space but not in time (you need the steady state of the recurrent state after the first stage at the end of the second stage to get your gradients). I think this is fair, and the need for biological plausibility is well motivated in the beginning of the work. My first issue is with the proof for the equivalence between EP and C-EP. This is done by taking the limit of eta goes to 0. I think I must be missing something. But the proof relies on eta being small enough such that 	heta_i = 	heta (i.e. theta does not change). Given this state evolves the same way as for EP, because we basically not changing theta. Yet the crux of my issue is exactly here. The proof relies on the fact that we don*t change theta. So then when you converged on the second phase, isn*t theta the same as theta_0? So you haven*t actually learned anything!? Basically looking at the delta isn*t this just misleading? Ok lets assume that on the last step you allow yourself to change eta to be non-zero. (I.e. we are just after the delta in theta, and what to show we can get the same delta in theta as EP which is how the proof is phrased). Then in that difference aren*t you looking at s_{t+1} and s_t rather than s_{t+1} and s_0, which is what EP would do? In EP you have s^_x0008_eta_* - s_*. This is not what you get if you don*t update theta and apply C-EP? I think there might be something I*m missing about the mathematical argument here. At a higher-level question, we talk about the transition function F as being a gradient vector field, i.e. there exist a phi such that F is d phi/d theta. Why is this assumption biologically plausable ? Parametrizing gradient vector fields in general is far from trivial, and require very specific structure of the neural implementation of F to be true. Several works have looked at parametrizing gradient vector fields (https://arxiv.org/abs/1906.01563, https://arxiv.org/pdf/1608.05343.pdf) and the answer is that without parametrizing it by actually taking the gradient of a function there is not much of a choice. Incidentally, here we exploit that F = sigma (Wx), with W symmetric. This is a paramtrization of a gradient vector field, i.e. of xU, where UU^T =W I think. But if you want to make F deep than it becomes non-trivial to restrict it to gradient vector field. Is the assumption that we never want to move away from vanilla RNNs? And W symmetric is also not biologically plausible. In C-EP you say is not needed to be symmetric, but that implicitly means there exist no phi and everything that follows breaks, no? I*m also confused by how one has access to d phi / ds and d phi / d theta. Again I feel like I*m missing information and the formalism is not introduced in a way that it is easy to parse. My understand is that you have an RNN that updates the state s. And the transfer function of this RNN is meant to be d phi / ds, which is trues if the recurrent weight is symmetric. Fine. But then why do we have access to d phi/ dtheta? Who is this function? Is the assumption that d s / dtheta is something we can compute in a biologically plausible way? Is this something that is obvious?