# Scaling Equilibrium Propagation to Deep ConvNets
by Drastically Reducing its Gradient Estimator Bias

Axel Laborieux 1,∗, Maxence Ernoult1,2,∗, Benjamin Scellier3,†,
Yoshua Bengio3,4, Julie Grollier2, Damien Querlioz1
1Université Paris-Saclay, CNRS, C2N, 91120, Palaiseau, France
2Unité Mixte de Physique, CNRS, Thales, Université Paris-Saclay
3Mila, Université de Montréal
4Canadian Institute for Advanced Research
† Currently at Google
* Corresponding authors: {axel.laborieux<EMAIL_ADDRESS>
###### Abstract

Equilibrium Propagation (EP) is a biologically-inspired counterpart of
Backpropagation Through Time (BPTT) which, owing to its strong theoretical
guarantees and the locality in space of its learning rule, fosters the design
of energy-efficient hardware dedicated to learning. In practice, however, EP
does not scale to visual tasks harder than MNIST. In this work, we show that a
bias in the gradient estimate of EP, inherent in the use of finite nudging, is
responsible for this phenomenon and that cancelling it allows training deep
ConvNets by EP, including architectures with distinct forward and backward
connections. These results highlight EP as a scalable approach to compute
error gradients in deep neural networks, thereby motivating its hardware
implementation.

#### Introduction

How synapses in hierarchical neural circuits are adjusted throughout learning
a task remains a challenging question called the credit assignment problem
(Richards et al., 2019). Equilibrium Propagation (EP) (Scellier and Bengio,
2017) provides a biologically plausible solution to this problem in artificial
neural networks. EP is an algorithm for convergent RNNs which, by definition,
are given a static input and whose recurrent dynamics converge to a steady
state corresponding to the prediction of the network. EP proceeds in two
phases which both involve the same dynamics. First, the network relaxes to a
steady state, then the output layer is nudged towards a ground-truth target
until a second steady state is reached. During the second phase, the
perturbation at the output propagates to upstream layers, creating local error
signals that match exactly those computed by Backpropagation Through Time
(BPTT) (Ernoult et al., 2019). The spatial locality of the learning rule
prescribed by EP is highly attractive for designing energy-efficient
“neuromorphic” hardware implementations of gradient-based learning algorithms.
Deep learning models consume orders of magnitude more energy than the brain to
learn cognitive owing to the physical separation of memory and computation on
conventional hardware (Strubell et al., 2019; Markovic et al., 2020). On (non
Von Neumann) architectures where the memory is brought at the location of
computation, the implementation of EP could potentially be more energy
efficient than the one of Backpropagation on GPUs by at least two order of
magnitudes (Ambrogio et al., 2018).

However, previous works on EP (Scellier and Bengio, 2017; O’Connor et al.,
2018, 2019; Ernoult et al., 2019, 2020) limited their experiments to the MNIST
classification task and to shallow network architectures. Despite the
theoretical guarantees of EP, the literature suggests that no implementation
of EP has thus far succeeded to match the performance of standard deep
learning approaches to train deep networks on challenging visual tasks. In
this work, we show that performing the gradient computation phase of EP with
nudging strength of constant sign induces a systematic first order bias in the
EP gradient estimate which, once cancelled, unlocks the training of deep
ConvNets. By introducing a new method to estimate the gradient of the loss
based on three steady states instead of two, we achieve 11.68% test error on
CIFAR-10, with only 0.6 % accuracy degradation with respect to BPTT. Standard
EP with two steady states yields 86.64% test error. We also propose to
implement the neural network predictor as an external softmax readout,
subsequently allowing us to use the cross-entropy loss, contrary to previous
approaches using the squared error loss. Finally, based on ideas of Scellier
et al. (2018) and Kolen and Pollack (1994), we adapt the learning rule of EP
for architectures with distinct forward and backward connections, yielding
only 1.5% accuracy degradation on CIFAR-10 compared to bidirectional
connections.

Figure 1: (a), (b) Comparison between BP and EP in terms of computation.
$a_{L}$, $s_{L}$, $\delta_{L}$ denote respectively the pre-activation, the
activation and the error. (c) Schematic of the recurrent convolutional
architecture. (d) Comparison of the former ($\hat{\nabla}^{\rm EP}$) and of
the new ($\hat{\nabla}^{\rm EP}_{\rm sym}$) EP estimates with the BPTT
gradient ($\hat{\nabla}^{\rm BPTT}$) as a function of time during the second
phase (more instances in Appendix H).

#### Background

Learning in hardware. The challenge of implementing BP on a neuromorphic
system for training is due to its distinct modes of operation during forward
and backward passes (Fig. 1(a)). During the forward pass, neural activations
as well as their point derivatives have to be computed and stored for the
backward pass (Fumarola et al., 2016). Error signals computed with BP during
training also require high precision computation, which does not allow taking
advantage of neuromorphic hardware optimally (Narayanan et al., 2017).
Conversely, EP (Fig. 1(b)) uses the same mode of operation for the two phases:
the neural activations evolve along recurrent dynamics through a unique
bidirectional pathway. The error signal in the second phase is thus implicitly
propagated to upstream layers through the relaxation of neural activations.
Owing to these compelling assets, hardware designers are proposing novel
circuits for implementing energy-efficient on-chip learning with EP (Zoppo et
al., 2020; Kendall et al., 2020; Ji and Gross, 2020; Foroushani et al., 2020).

Equilibrium Propagation (EP) (Scellier and Bengio, 2017). We consider the
setting of supervised learning where $x$ is the static input and $y$ is the
target. The goal of learning is to optimize the network parameters $\theta$ to
minimize the loss $\mathcal{L}^{*}=\ell(s_{*},y)$ (Appendix A). Denoting
$s_{t}$ the state of the neurons at time step $t$, during the first (“free”)
phase, the RNN evolves according to the following dynamics:

$s_{t+1}=\frac{\partial\Phi}{\partial s}(x,s_{t},\theta),$ (1)

where $\Phi$ denotes the scalar primitive of the dynamics, and which converges
to a steady-state $s_{*}$. These discrete-time equations, introduced in
(Ernoult et al., 2019), are a simplification of the original version of EP for
real-time dynamical systems (Scellier and Bengio, 2017), where the dynamics
derive from an “energy function” $E$, which plays a similar role to the
primitive function $\Phi$. During the second (“nudged”) phase of EP, the
system evolves with an additional nudging term added to the dynamics, $\beta$
being a small scalar factor:

$s_{t+1}^{\beta}=\frac{\partial\Phi}{\partial
s}(x,s_{t}^{\beta},\theta)-\beta\frac{\partial\ell}{\partial
s}(s_{t}^{\beta},y),$ (2)

and converges to a steady-state $s_{*}^{\beta}$. The gradient estimate of EP
then reads

$\quad\widehat{\nabla}^{\rm{EP}}(\beta)\stackrel{{\scriptstyle\Delta}}{{=}}\frac{1}{\beta}\left(\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{\beta},\theta)-\frac{\partial\Phi}{\partial\theta}(x,s_{*},\theta)\right).$
(3)

A theorem derived by Ernoult et al. (2019), inspired from Scellier and Bengio
(2019), shows that, provided convergence in the first phase has been reached,
these gradient estimates given by EP truncated at any time step $t$ of the
second phase match, in the limit $\beta\to 0$, those computed by BPTT
truncated at time step $T-t$ (Appendix B).

Convolutional Architectures for Convergent RNNs. A convolutional architecture
for convergent RNNs with static input (Fig. 1(c)) was introduced by Ernoult et
al. (2019) and successfully trained with EP on the MNIST dataset, with
accuracy approaching the one of BPTT. Denoting $w_{n}$ a convolutional or
fully connected layer, $\star$ and $\mathcal{P}$ convolution and pooling
respectively (with $\tilde{w}$, $\mathcal{P}^{-1}$ denoting their respective
inverses) and $\sigma$ an activation function, the dynamics read (details in
Appendix E):

$\displaystyle\left\\{\begin{array}[]{ll}s^{n}_{t+1}&=\sigma\left(\mathcal{P}\left(w_{n-1}\star
s^{n-1}_{t}\right)+\tilde{w}_{n+1}\star\mathcal{P}^{-1}\left(s^{n+1}_{t}\right)\right),\qquad\mbox{for
convolutional layers,}\\\ s^{n}_{t+1}&=\sigma\left(w_{n-1}\cdot
s^{n-1}_{t}+w^{\top}_{n+1}\cdot s^{n+1}_{t}\right),\qquad\mbox{for fully
connected layers.}\end{array}\right.$ (6)

Equilibrium Propagation with unidirectional synaptic connections. Previous
works have proposed a more general formulation of EP, called _Vector Field_
(VF), where the dynamics need not derive from a primitive scalar function.
This allows training networks with distinct forward and backward connections
(Scellier et al., 2018; Ernoult et al., 2020) (see Appendix E.4 for an
explicit expression of the VF gradient estimate), further enhancing the
biological plausibility of EP.

#### Contributions of this work : Scaling EP Training

Reducing bias and variance in the gradient estimate of the loss function.
Scellier and Bengio (2017) motivate their two-phases procedure with $\beta=0$
and $\beta>0$ by using $\hat{\nabla}^{\rm EP}$ (Eq. 3) to compute an estimate
of the derivative
$\left.\frac{d}{d\beta}\right|_{\beta=0}\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{\beta},\theta)$
which they showed to be equal to
$-\frac{\partial\mathcal{L}_{*}}{\partial\theta}$. However, the use of
$\beta>0$ induces a systematic first order bias in this estimation, which we
propose to eliminate by performing a third phase with $-\beta$ as the nudging
factor (proof in Appendix C), keeping the first and second phases unchanged,
along with the following symmetric difference estimate:

$\widehat{\nabla}^{\rm{EP}}_{\rm
sym}(\beta)\stackrel{{\scriptstyle\Delta}}{{=}}\frac{1}{2\beta}\left(\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{\beta},\theta)-\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{-\beta},\theta)\right).$
(7)

We call $\widehat{\nabla}^{\rm{EP}}(\beta)$ and
$\widehat{\nabla}^{\rm{EP}}_{\rm sym}(\beta)$ the one-sided and symmetric EP
gradient estimates respectively. Note that the first-order term of
$\widehat{\nabla}^{\rm{EP}}(\beta)$ could also be cancelled out on average, by
choosing the sign of $\beta$ at random with even probability, which we call
“Random Sign”. We see qualitatively in Fig. 1(d) that
$\widehat{\nabla}^{\rm{EP}}_{\rm sym}(\beta)$ gives a better estimate of
$\hat{\nabla}^{\rm BPTT}$ (see Appendix B for a proof) than
$\widehat{\nabla}^{\rm{EP}}(\beta)$.

Changing the loss function. Previous implementations of EP used the squared
error loss as, in this setting, the output neurons can be viewed as a part of
the system which bidirectionally interacts with neighboring neurons and the
nudging term as an elastic force. We propose here an alternative approach,
where the output $\hat{y}_{t}$ is no longer a subset of neurons of the system
but is instead implemented as an external read-out, consisting of the
composition of a read-out weight matrix $w_{\rm out}$ with the softmax
function: $\widehat{y}_{t}=\mbox{softmax}(w_{\rm out}\cdot s_{t})$. This
approach enables the use of a cross-entropy cost function (see Appendix D for
implementation details).

Changing the learning rule of vector field EP. Until now, training experiments
of unidirectional weights EP have performed worse than bidirectional weights
EP (Ernoult et al., 2020). In this work, therefore, we tailor a new learning
rule for unidirectional weights (described in Appendix E.4), where the forward
and backward weights undergo the same weight updates, incorporating an equal
leakage term. This way, forward and backward weights, although they are
independently initialized, naturally converge to identical values throughout
the learning process, a technique which is adapted from Kolen and Pollack
(1994) and helped improve the performance of Feedback Alignment in Deep
ConvNets (Akrout et al., 2019). We adapt the three-phases procedure detailed
above to this setting to compute the common update of forward and backward
weights, thereby defining the gradient estimate $\hat{\nabla}^{\rm KP-VF}_{\rm
sym}$. We also define $\hat{\nabla}^{\rm VF}_{\rm sym}$ the gradient estimate
obtained by applying the three-phases procedure to the former VF approach of
Scellier et al. (2018); Ernoult et al. (2020).

Table 1: Accuracy comparison on CIFAR-10 between BPTT and EP with several gradient estimation schemes. Note that the different gradient estimates only apply to EP. We indicate over five trials the mean and standard deviation in parenthesis for the test error, and the mean train error. Loss Function | EP Gradient | EP Error (%) | BPTT Error (%)
---|---|---|---
Estimate | Test | Train | Test | Train
Squared Error | 2-Phase / $\hat{\nabla}^{\rm EP}$ |  $86.64$ $(5.82)$ | $84.90$ | $11.10$ $(0.21)$ | $3.69$
Random Sign |  $21.55$ $(20.00)$ | $20.01$
3-Phase / $\hat{\nabla}^{\rm EP}_{\rm sym}$ |  $12.45$ $(0.18)$ | $7.83$
Cross-Ent. | 3-Phase / $\hat{\nabla}^{\rm EP}_{\rm sym}$ |  $\mathbf{11.68}$ $\mathbf{(0.17)}$ | $\mathbf{4.98}$ |  $11.12$ $(0.21)$ | $2.19$
Cross-Ent. (Dropout) | 3-Phase / $\hat{\nabla}^{\rm EP}_{\rm sym}$ |  $11.87$ $(0.29)$ | $6.46$ |  $10.72$ $(0.06)$ | $2.99$
Cross-Ent. | 3-Phase / $\widehat{\nabla}^{\rm{VF}}_{\rm sym}$ |  $75.47$ $(4.72)$ | $78.04$ | $9.46$ $(0.17)$ | $0.80$
3-Phase / $\widehat{\nabla}^{\rm{KP-VF}}_{\rm sym}$ |  $\mathbf{13.15}$ $\mathbf{(0.49)}$ | $8.87$

Experimental results. In Table 1, we compare the accuracy achieved by the
ConvNet for each EP gradient estimate with the accuracy achieved by BPTT. On
bidirectional weight architectures, the one-sided gradient estimate
$\hat{\nabla}^{\rm EP}$ leads to unstable training behavior where the network
is unable to fit the data. When the bias in the gradient estimate is averaged
out by choosing at random the sign of $\beta$, the average test error over
five runs goes down to $21.55\%$. However, one run among the five yielded
instability similar to the one-sided estimate, whereas the four remaining runs
lead to $12.61\%$ test error and $8.64\%$ train error. This method for
estimating the loss gradient thus presents high variance, further experiments
in Appendix F confirm this trend. Conversely, the symmetric estimate
$\hat{\nabla}^{\rm EP}_{\rm sym}$ enables EP to consistently reach $12.45\%$
test error, with only $1.35\%$ degradation with respect to BPTT. Therefore,
removing the first-order error term in the gradient estimate is critical for
scaling to deeper architectures. However, proceeding to this end
deterministically (with three phases) rather than stochastically (random sign)
seems to be more reliable.

The new readout scheme introduced to optimize the cross-entropy loss function
enables EP to narrow the performance gap with BPTT down to $0.56\%$ while
outperforming the Squared Error setting by $0.77\%$. Finally, we adapted
dropout (Srivastava et al., 2014) to convergent RNNs (see Appendix F.3 for
implementation details) to see if the performance could be improved further.
However, we can observe in Table 1 that contrary to BPTT, the EP test error is
not improved by adding dropout, which we hypothesize is due to the residual
estimation bias of the BPTT gradients by EP.

In the situation when the architecture uses distinct forward and backward
weights, we find that the traditional estimate
$\widehat{\nabla}^{\rm{VF}}_{\rm sym}(\beta)$ leads to a poor accuracy with
$75.47\%$ test-error and simultaneously observe that forward and backward
weights do not align well. Conversely, when using our new estimate
$\widehat{\nabla}^{\rm{KP-VF}}_{\rm sym}(\beta)$, a good accuracy is recovered
with $1.5\%$ degradation with respect to the architecture with bidirectional
connections, and a $3\%$ degradation with respect to BPTT. In this case,
forward and backward weights are aligned by epoch 50, as observed in Appendix
G. These results suggest that enhancing forward and backward weights alignment
also helps EP training in deep ConvNets.

Discussion. In comparison with conventional implementations of EP, our results
unveil the necessity to compute better gradient estimates in order to scale EP
to deep ConvNets on hard visual tasks. Keeping the first order bias in the
gradient estimate of EP, as done traditionally, severely impedes the training
of these architectures and, conversely, removing it brings EP accuracy on
CIFAR-10 close to the one achieved by BPTT. Employing a new training technique
that preserves the spatial locality of EP computations, our results extend to
architectures with distinct forward and backward synaptic connections. We only
observe a $1.5\%$ performance degradation with respect to the bidirectional
architecture. Our three steady states-based gradient estimate comes at a
computational cost since one more phase is needed with regards to the
conventional implementation. In the longer run, the full potential of EP will
be best envisioned on neuromorphic hardware, which can sustain fast analog
device physics and use them to implement the dynamics of EP intrinsically
(Romera et al., 2018; Ambrogio et al., 2018). Our prescription to run two
nudging phases with opposite nudging strengths could be naturally implemented
on such systems, which often function differentially to cancel device inherent
biases (Bocquet et al., 2018).

### Broader Impact

This work may have a long-term impact on the design of energy-efficient
hardware leveraging the physics of the device to perform learning. The
demonstration that EP can scale to deep networks may also provide insights to
neuroscientists to understand the mechanisms of credit assignment in the
brain. Due to the long term nature of this impact, the positive and negative
outcomes of this work cannot yet be stated.

### References

  * Akrout et al. [2019] M. Akrout, C. Wilson, P. Humphreys, T. Lillicrap, and D. B. Tweed.  Deep learning without weight transport.  In _Advances in Neural Information Processing Systems_ , pages 974–982, 2019.
  * Ambrogio et al. [2018] S. Ambrogio, P. Narayanan, H. Tsai, R. M. Shelby, I. Boybat, C. di Nolfo, S. Sidler, M. Giordano, M. Bodini, N. C. Farinha, et al.  Equivalent-accuracy accelerated neural-network training using analogue memory.  _Nature_ , 558(7708):60–67, 2018.
  * Bocquet et al. [2018] M. Bocquet, T. Hirztlin, J.-O. Klein, E. Nowak, E. Vianello, J.-M. Portal, and D. Querlioz.  In-memory and error-immune differential rram implementation of binarized deep neural networks.  In _2018 IEEE International Electron Devices Meeting (IEDM)_ , pages 20–6. IEEE, 2018.
  * Ernoult et al. [2019] M. Ernoult, J. Grollier, D. Querlioz, Y. Bengio, and B. Scellier.  Updates of equilibrium prop match gradients of backprop through time in an rnn with static input.  _arXiv preprint arXiv:1905.13633_ , 2019.
  * Ernoult et al. [2020] M. Ernoult, J. Grollier, D. Querlioz, Y. Bengio, and B. Scellier.  Equilibrium propagation with continual weight updates.  _arXiv preprint arXiv:2005.04168_ , 2020.
  * Foroushani et al. [2020] A. N. Foroushani, H. Assaf, F. H. Noshahr, Y. Savaria, and M. Sawan.  Analog circuits to accelerate the relaxation process in the equilibrium propagation algorithm.  In _2020 IEEE International Symposium on Circuits and Systems (ISCAS)_ , pages 1–5. IEEE, 2020.
  * Fumarola et al. [2016] A. Fumarola, P. Narayanan, L. L. Sanches, S. Sidler, J. Jang, K. Moon, R. M. Shelby, H. Hwang, and G. W. Burr.  Accelerating machine learning with non-volatile memory: Exploring device and circuit tradeoffs.  In _2016 IEEE International Conference on Rebooting Computing (ICRC)_ , pages 1–8. Ieee, 2016.
  * He et al. [2015] K. He, X. Zhang, S. Ren, and J. Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  In _Proceedings of the IEEE international conference on computer vision_ , pages 1026–1034, 2015.
  * Ji and Gross [2020] Z. Ji and W. Gross.  Towards efficient on-chip learning using equilibrium propagation.  In _2020 IEEE International Symposium on Circuits and Systems (ISCAS)_ , pages 1–5. IEEE, 2020.
  * Kendall et al. [2020] J. Kendall, R. Pantone, K. Manickavasagam, Y. Bengio, and B. Scellier.  Training end-to-end analog neural networks with equilibrium propagation.  _arXiv preprint arXiv:2006.01981_ , 2020.
  * Kolen and Pollack [1994] J. F. Kolen and J. B. Pollack.  Backpropagation without weight transport.  In _Proceedings of 1994 IEEE International Conference on Neural Networks (ICNN’94)_ , volume 3, pages 1375–1380. IEEE, 1994.
  * Loshchilov and Hutter [2016] I. Loshchilov and F. Hutter.  Sgdr: Stochastic gradient descent with warm restarts.  _arXiv preprint arXiv:1608.03983_ , 2016.
  * Markovic et al. [2020] D. Markovic, A. Mizrahi, D. Querlioz, and J. Grollier.  Physics for neuromorphic computing.  _arXiv preprint arXiv:2003.04711_ , 2020.
  * Narayanan et al. [2017] P. Narayanan, A. Fumarola, L. L. Sanches, K. Hosokawa, S. C. Lewis, R. M. Shelby, and G. W. Burr.  Toward on-chip acceleration of the backpropagation algorithm using nonvolatile memory.  _IBM Journal of Research and Development_ , 61(4/5):11–1, 2017.
  * O’Connor et al. [2018] P. O’Connor, E. Gavves, and M. Welling.  Initialized equilibrium propagation for backprop-free training.  2018\.
  * O’Connor et al. [2019] P. O’Connor, E. Gavves, and M. Welling.  Training a spiking neural network with equilibrium propagation.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 1516–1523, 2019.
  * Paszke et al. [2017] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer.  Automatic differentiation in pytorch.  2017\.
  * Richards et al. [2019] B. A. Richards, T. P. Lillicrap, P. Beaudoin, Y. Bengio, R. Bogacz, A. Christensen, C. Clopath, R. P. Costa, A. de Berker, S. Ganguli, et al.  A deep learning framework for neuroscience.  _Nature neuroscience_ , 22(11):1761–1770, 2019\.
  * Romera et al. [2018] M. Romera, P. Talatchian, S. Tsunegi, F. A. Araujo, V. Cros, P. Bortolotti, J. Trastoy, K. Yakushiji, A. Fukushima, H. Kubota, et al.  Vowel recognition with four coupled spin-torque nano-oscillators.  _Nature_ , 563(7730):230–234, 2018.
  * Scellier and Bengio [2017] B. Scellier and Y. Bengio.  Equilibrium propagation: Bridging the gap between energy-based models and backpropagation.  _Frontiers in computational neuroscience_ , 11:24, 2017\.
  * Scellier and Bengio [2019] B. Scellier and Y. Bengio.  Equivalence of equilibrium propagation and recurrent backpropagation.  _Neural computation_ , 31(2):312–329, 2019.
  * Scellier et al. [2018] B. Scellier, A. Goyal, J. Binas, T. Mesnard, and Y. Bengio.  Generalization of equilibrium propagation to vector field dynamics.  _arXiv preprint arXiv:1808.04873_ , 2018.
  * Srivastava et al. [2014] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  _The journal of machine learning research_ , 15(1):1929–1958, 2014.
  * Strubell et al. [2019] E. Strubell, A. Ganesh, and A. McCallum.  Energy and policy considerations for deep learning in nlp.  _arXiv preprint arXiv:1906.02243_ , 2019.
  * Zoppo et al. [2020] G. Zoppo, F. Marrone, and F. Corinto.  Equilibrium propagation for memristor-based recurrent neural networks.  _Frontiers in neuroscience_ , 14:240, 2020.

## Appendix

###  Appendix A Pseudo code

####  A.1 Random one-sided estimation of the loss gradient

In this appendix, we define the random one-sided estimation used in this work
and by Scellier and Bengio [2017], Ernoult et al. [2020].

_Input_ : $x$, $y$, $\theta$, $\eta$.
_Output_ : $\theta$.

Algorithm 1 EP with random one-sided estimation of the loss gradient. We omit
the activation function $\sigma$ for clarity.

1:$s_{0}\leftarrow 0$

2:for $t=0$ to $T$ do $\triangleright$ First phase.

3: $s_{t+1}\leftarrow\frac{\partial\Phi}{\partial s}(x,s_{t},\theta)$

4:end for

5:$s_{*}\leftarrow s_{T}$

6:$\beta\leftarrow\beta\times{\rm{Bernoulli}(1,-1)}$ $\triangleright$ Random
sign.

7:$s_{0}^{\beta}\leftarrow s_{*}$

8:for $t=0$ to $K$ do $\triangleright$ Second phase.

9: $s_{t+1}^{\beta}\leftarrow\frac{\partial\Phi}{\partial
s}(x,s_{t}^{\beta},\theta)-\beta\frac{\partial\ell}{\partial
s}(s_{t}^{\beta},y)$

10:end for

11:$s_{*}^{\beta}\leftarrow s_{K}^{\beta}$

12:$\nabla_{\theta}^{\rm{EP}}\leftarrow\frac{1}{\beta}\left(\frac{\partial\Phi}{\partial\theta}(s_{*}^{\beta},\theta)-\frac{\partial\Phi}{\partial\theta}(s_{*},\theta)\right)$

13:$\theta\leftarrow\theta+\eta\nabla_{\theta}^{\rm{EP}}$

14:return $\theta$

####  A.2 Symmetric difference estimation of the loss gradient

In this appendix, we define the estimation procedure using a symmetric
difference estimate introduced in this work.

_Input_ : $x$, $y$, $\theta$, $\eta$.
_Output_ : $\theta$.

Algorithm 2 EP with symmetric difference estimation of the loss gradient. We
omit the activation function $\sigma$ for clarity.

1:$s_{0}\leftarrow 0$

2:for $t=0$ to $T$ do

3: $s_{t+1}\leftarrow\frac{\partial\Phi}{\partial s}(x,s_{t},\theta)$
$\triangleright$ First phase.

4:end for

5:$s_{*}\leftarrow s_{T}$ $\triangleright$ Store the free steady state.

6:$s_{0}^{\beta}\leftarrow s_{*}$

7:for $t=0$ to $K$ do

8: $s_{t+1}^{\beta}\leftarrow\frac{\partial\Phi}{\partial
s}(x,s_{t}^{\beta},\theta)-\beta\frac{\partial\ell}{\partial
s}(s_{t}^{\beta},y)$ $\triangleright$ Second phase.

9:end for

10:$s_{*}^{\beta}\leftarrow s_{K}^{\beta}$

11:$s_{0}^{-\beta}\leftarrow s_{*}$ $\triangleright$ Back to the free steady
state.

12:for $t=0$ to $K$ do

13: $s_{t+1}^{-\beta}\leftarrow\frac{\partial\Phi}{\partial
s}(x,s_{t}^{-\beta},\theta)+\beta\frac{\partial\ell}{\partial
s}(s_{t}^{-\beta},y)$ $\triangleright$ Third phase.

14:end for

15:$s_{*}^{-\beta}\leftarrow s_{K}^{-\beta}$

16:$\widehat{\nabla}_{\theta}^{\rm{EP}}\leftarrow\frac{1}{2\beta}\left(\frac{\partial\Phi}{\partial\theta}(s_{*}^{\beta},\theta)-\frac{\partial\Phi}{\partial\theta}(s_{*}^{-\beta},\theta)\right)$

17:$\theta\leftarrow\theta+\eta\widehat{\nabla}_{\theta}^{\rm{EP}}$

18:return $\theta$

###  Appendix B Training RNNs with BPTT

The convergent RNNs considered by EP can also be trained by Backpropagation
Through Time (BPTT). In this context, BPTT consists in performing the first
phase for $T$ time steps until the network reaches the steady state
$s_{T}=s_{*}$, computing the loss at the final time step and subsequently
backpropagating the gradients through the computational graph of the first
phase.

We write $\nabla^{\rm{BPTT}}(t)$ the gradient computed by BPTT truncated to
the last $t$ time steps ($T-t,\ldots,T$). To derive it in function of loss
gradients, let us rewrite Eq. (1) as $s_{t+1}=\frac{\partial\Phi}{\partial
s}(x,s_{t},\theta_{t}=\theta)$, where $\theta_{t}$ denotes the parameter at
time step $t$, the value $\theta$ being shared across all time steps. We
consider the loss after $T$ time steps $\mathcal{L}=\ell(s_{T},y)$. Rewriting
the dynamics in such a way enables us to define
$\frac{\partial\mathcal{L}}{\partial\theta_{t}}$ as the sensitivity of the
loss with respect to $\theta_{t}$, when
$\theta_{0},\ldots,\theta_{t-1},\theta_{t+1},\ldots,\theta_{T-1}$ remain fixed
(set to the value $\theta$). With these notations, the gradient computed by
BPTT truncated to the last $t$ time steps is

$\nabla^{\rm{BPTT}}(t)=\frac{\partial\mathcal{L}}{\partial\theta_{T-t}}+\ldots+\frac{\partial\mathcal{L}}{\partial\theta_{T-1}}.$
(8)

###  Appendix C Error terms in the estimates of the loss gradient

In this appendix, we prove Lemma 2 which shows that
$\widehat{\nabla}^{\rm{EP}}_{\rm sym}(\beta)$ is a better estimate of
$-\frac{\partial\mathcal{L}^{*}}{\partial\theta}$ than
$\widehat{\nabla}^{\rm{EP}}(\beta)$. First, we recall the theorem proved in
Scellier and Bengio [2017].

######  Theorem 1 (Scellier and Bengio [2017]).

$\left.\frac{d}{d\beta}\right|_{\beta=0}\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{\beta},\theta)=-\frac{\partial\mathcal{L}^{*}}{\partial\theta}.$
(9)

We also recall that the two estimates (one-sided and symmetric) are, by
definition:

$\displaystyle\widehat{\nabla}^{\rm{EP}}(\beta)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}\frac{1}{\beta}\left(\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{\beta},\theta)-\frac{\partial\Phi}{\partial\theta}(x,s_{*},\theta)\right),$
$\displaystyle\widehat{\nabla}^{\rm{EP}}_{\rm sym}(\beta)$
$\displaystyle\stackrel{{\scriptstyle\Delta}}{{=}}\frac{1}{2\beta}\left(\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{\beta},\theta)-\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{-\beta},\theta)\right).$

Finally we recall Lemma 2, for readability.

######  Lemma 2.

Provided the function
$\beta\mapsto\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{\beta},\theta)$ is
three times differentiable, we have, as $\beta\to 0$:

$\displaystyle\widehat{\nabla}^{\rm{EP}}(\beta)$
$\displaystyle=-\frac{\partial\mathcal{L}^{*}}{\partial\theta}+\frac{\beta}{2}\left.\frac{d^{2}}{d\beta^{2}}\right|_{\beta=0}\frac{\partial\Phi}{\partial\theta}(s_{*}^{\beta},\theta)+O(\beta^{2}),$
$\displaystyle\widehat{\nabla}^{\rm{EP}}_{\rm sym}(\beta)$
$\displaystyle=-\frac{\partial\mathcal{L}^{*}}{\partial\theta}+O(\beta^{2}).$

###### Proof of Lemma 2.

Let us define

$f(\beta)\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{\beta},\theta).$

The formula of Theorem 9 rewrites

$f^{\prime}(0)=-\frac{\partial\mathcal{L}^{*}}{\partial\theta}.$

As $\beta\to 0$, we have the Taylor expansion

$f(\beta)=f(0)+\beta
f^{\prime}(0)+\frac{\beta^{2}}{2}f^{\prime\prime}(0)+O(\beta^{3}).$ (10)

With these notations, the one-sided estimate reads

$\displaystyle\widehat{\nabla}^{\rm{EP}}(\beta)$
$\displaystyle=\frac{1}{\beta}\left(f(\beta)-f(0)\right)$
$\displaystyle=f^{\prime}(0)+\frac{\beta}{2}f^{\prime\prime}(0)+O(\beta^{2})$
$\displaystyle=-\frac{\partial\mathcal{L}^{*}}{\partial\theta}+\frac{\beta}{2}\left.\frac{d^{2}}{d\beta^{2}}\right|_{\beta=0}\frac{\partial\Phi}{\partial\theta}(x,s_{*}^{\beta},\theta)+O(\beta^{2}).$

We can also write a Taylor expansion around $0$ at the point $-\beta$. We have

$f(-\beta)=f(0)-\beta
f^{\prime}(0)+\frac{\beta^{2}}{2}f^{\prime\prime}(0)+O(\beta^{3}).$ (11)

Subtracting Eq. 11 from Eq. 10, we can rewrite the symmetric difference
estimate as

$\displaystyle\widehat{\nabla}^{\rm{EP}}_{\rm sym}(\beta)$
$\displaystyle=\frac{1}{2\beta}\left(f(\beta)-f(-\beta)\right)$
$\displaystyle=f^{\prime}(0)+O(\beta^{2})$
$\displaystyle=-\frac{\partial\mathcal{L}^{*}}{\partial\theta}+O(\beta^{2}).$

The derivative to the third order of $f$ is only used to get the
$O(\beta^{3})$ term in the expansion Eq. (10), it can be changed into
$o(\beta^{2})$ if we only assume $f$ twice differentiable. ∎

###  Appendix D Changing the loss function

We introduce a novel architecture to optimize the cross-entropy loss with EP,
narrowing the gap with conventional deep learning architectures for
classification tasks. In the next paragraph, we denote $\widehat{y}$ the set
of neurons that carries out the prediction of the neural network.

Figure 2: Free dynamics of the architectures used for the two loss functions
where the blue frame delimits the system. (a) Squared Error loss function. The
usual setting where the predictor $\hat{y}$ (in red) takes part in the free
dynamics of the neural network through bidirectional synaptic connections. (b)
Cross-Entropy loss function. The new approach proposed in this work where the
predictor $\hat{y}$ (also in red) is no longer involved in the system free
dynamics and is implemented as a softmax readout.

###### Squared Error loss function.

Previous implementations of EP used the squared error loss. Using this loss
function for EP is natural, as in this setting, the output $\widehat{y}$ is
viewed as a part of $s$ (the state variable of the network), which can
influence the state of the network through bidirectional synaptic connections
(see Fig. 2). The state of the network is of the form
$s=(s^{1},\dots,s^{N},\widehat{y})$ where $h=(s^{1},\dots,s^{N})$ represent
the “hidden layers”, and the corresponding cost function is

$\ell(\widehat{y},y)=\frac{1}{2}\left\|\widehat{y}-y\right\|^{2}.$ (12)

The second phase dynamics of the hidden state and output layer given by Eq.
(2) read, in this context:

$h_{t+1}^{\beta}=\frac{\partial\Phi}{\partial
h}(x,h_{t}^{\beta},\widehat{y}_{t}^{\beta},\theta),\qquad\widehat{y}_{t+1}^{\beta}=\frac{\partial\Phi}{\partial\widehat{y}}(x,h_{t}^{\beta},\widehat{y}_{t}^{\beta},\theta)+\beta\;(y-\widehat{y}_{t}^{\beta}).$
(13)

###### Softmax readout, Cross-Entropy loss function.

In this paper, we propose an alternative approach, where the output
$\widehat{y}$ is not a part of the state variable $s$ but is instead
implemented as a read-out (see Fig. 2), which is a function of $s$ and of a
weight matrix $w_{\rm out}$ of size $\dim(y)\times\dim(s)$. In practice,
$w_{\rm out}$ reads out the last convolutional layer. At each time step $t$ we
define:

$\widehat{y}_{t}=\mbox{softmax}(w_{\rm out}\cdot s_{t}).$ (14)

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

$\ell(s,y,w_{\rm out})=-\sum_{c=1}^{C}y_{c}\log(\textrm{softmax}_{c}(w_{\rm
out}\cdot s)).$ (15)

Using $\frac{\partial\ell}{\partial s}(s,y,w_{\rm out})=w_{\rm
out}^{\top}\cdot\left(\textrm{softmax}(w_{\rm out}\cdot s)-y\right)$, the
second phase dynamics given by Eq. (2) read in this context:

$s_{t+1}^{\beta}=\frac{\partial\Phi}{\partial
s}(x,s_{t}^{\beta},\theta)+\beta\;w_{\rm
out}^{\top}\cdot\left(y-\widehat{y}_{t}^{\beta}\right).$ (16)

Note here that the loss $\mathcal{L}^{*}=\ell(s_{*},y,w_{\rm out})$ also
depends on the parameter $w_{\rm out}$. Appendix E.3.2 provides the learning
rule applied to $w_{\rm out}$.

###  Appendix E Convolutional RNNs

Throughout this section, $N^{\rm conv}$ and $N^{\rm fc}$ denote respectively
the number of convolutional layers and fully connected layers in the
convolutional RNN, and $N^{\rm tot}\stackrel{{\scriptstyle\Delta}}{{=}}N^{\rm
conv}+N^{\rm fc}$. The neuron layers are denoted by $s$ and range from
$s^{0}=x$ the input to the output $s^{N^{\rm tot}}$ in the case of squared
error, or $s^{N^{\rm tot}-1}$ in the case of softmax read-out.

####  E.1 Definition of the operations

In this subsection we detail the operations involved in the dynamics of a
convolutional RNN.

  * •

The 2-D convolution between $w$ with dimension $(C_{\rm out},C_{\rm in},F,F)$
and an input $x$ of dimensions $(C_{\rm in},H_{\rm in},W_{\rm in})$ and stride
one is a tensor $y$ of size $(C_{\rm out},H_{\rm out},W_{\rm out})$ defined
by:

$y_{c,h,w}=(w\star x)_{c,h,w}=B_{c}+\sum_{i=0}^{C_{\rm
in}-1}\sum_{j=0}^{F-1}\sum_{k=0}^{F-1}w_{c,i,j,k}x_{i,j+h,k+w},$ (17)

where $B_{c}$ is a channel-wise bias.

  * •

The 2-D transpose convolution of $y$ by $\tilde{w}$ is then defined in this
work as the gradient of the 2-D convolution with respect to its input:

$(\tilde{w}\star y)\stackrel{{\scriptstyle\Delta}}{{=}}\frac{\partial(w\star
x)}{\partial x}\cdot y$ (18)

  * •

The dot product “$\bullet$” generalized to pairs of tensors of same shape
$(C,H,W)$:

$a\bullet
b=\sum_{c=0}^{C-1}\sum_{h=0}^{H-1}\sum_{w=0}^{W-1}a_{c,h,w}b_{c,h,w}.$ (19)

  * •

The pooling operation $\mathcal{P}$ with stride $F$ and filter size $F$ of
$x$:

$\mathcal{P}_{F}(x)_{c,h,w}=\underset{i,j\in[0,F-1]}{\rm
max}\left\\{x_{c,F(h-1)+1+i,F(w-1)+1+j}\right\\},$ (20)

with relative indices of maximums within each pooling zone given by:

${\rm ind}_{\mathcal{P}}(x)_{c,h,w}=\underset{i,j\in[0,F-1]}{\rm
argmax}\left\\{x_{c,F(h-1)+1+i,F(w-1)+1+j}\right\\}=(i^{*}(x,h),j^{*}(x,w)).$
(21)

  * •

The unpooling operation $\mathcal{P}^{-1}$ of $y$ with indices ${\rm
ind}_{\mathcal{P}}(x)$ is then defined as:

$\mathcal{P}^{-1}(y,{\rm
ind}_{\mathcal{P}}(x))_{c,h,w}=\sum_{i,j}y_{c,i,j}\cdot\delta_{h,F(i-1)+1+i^{*}(x,h)}\cdot\delta_{w,F(j-1)+1+j^{*}(x,w)},$
(22)

which consists in filling a tensor with the same dimensions as $x$ with the
values of $y$ at the indices ${\rm ind}_{\mathcal{P}}(x)$, and zeroes
elsewhere. For notational convenience, we omit to write explicitly the
dependence on the indices except when appropriate.

  * •

The flattening operation $\mathcal{F}$ is defined as reshaping a tensor of
dimensions $(C,H,W)$ to $(1,CHW)$. We denote by $\mathcal{F}^{-1}$ its
inverse.

####  E.2 Definition of the primitive function

For notational simplicity here, whether $w_{n}$ is a convolutional layer or a
fully connected layer is implied by the operator, respectively $\star$ for
convolutions and $\cdot$ for linear layers. The primitive function can
therefore be defined as:

$\Phi(x,\\{s^{n}\\})=\sum_{n=0}^{N^{\rm
conv}-1}s^{n+1}\bullet\mathcal{P}\left(w_{n+1}\star
s^{n}\right)+\sum_{n=N^{\rm conv}}^{N^{\rm tot}-1}s^{n\top}\cdot w_{n+1}\cdot
s^{n+1},$ (23)

where $\bullet$ is the Euclidean scalar product generalized to pairs of
tensors with same arbitrary dimension, and $\mathcal{P}$ is a pooling
operation.

####  E.3 Convolutional RNNs with bidirectional connections

In this section, we write explicitly the dynamics and the learning rules
applied for the convolutional architecture with bidirectional connections, for
the Squared loss function and the Cross-Entropy loss function, for the one-
sided and symmetric estimates.

#####  E.3.1 Squared Error loss

###### Equations of the dynamics.

In this case, the dynamics read:

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle
s^{n+1}_{t+1}=\sigma\left(\mathcal{P}(w_{n+1}\star
s^{n}_{t})+\tilde{w}_{n+2}\star\mathcal{P}^{-1}(s^{n+2}_{t})\right),\qquad\forall
n\in[0,N^{\rm conv}-2]\\\ \displaystyle s^{N^{\rm
conv}}_{t+1}=\sigma\left(\mathcal{P}(w_{N^{\rm conv}}\star s^{N^{\rm
conv}-1}_{t})+\mathcal{F}^{-1}({w_{N^{\rm conv}+1}}^{\top}\cdot s^{N^{\rm
conv}+1}_{t})\right),\\\ \displaystyle s^{N^{\rm
conv}+1}_{t+1}=\sigma\left(w_{N^{\rm conv}+1}\cdot\mathcal{F}(s^{N^{\rm
conv}}_{t})+{w_{N^{\rm conv}+2}}^{\top}\cdot s^{N^{\rm conv}+2}_{t}\right),\\\
\displaystyle s_{t+1}^{n+1}=\sigma\left(w_{n+1}\cdot
s^{n}_{t}+{w_{n+2}}^{\top}\cdot s^{n+2}_{t}\right),\qquad\forall n\in[N^{\rm
conv}+1,N^{\rm tot}-2]\\\ \displaystyle s_{t+1}^{N^{\rm
tot}}=\sigma\left(w_{N^{\rm tot}}\cdot s^{N^{\rm
tot}-1}_{t}\right)+\beta(y-s^{N^{\rm tot}}),\quad\text{with $\beta=0$ during
the first phase,}\end{array}\right.$ (29)

where we take the convention $s^{0}=x$. In this case, we have
$\hat{y}=s_{t+1}^{N^{\rm tot}}$. Considering the function:

$\displaystyle\Phi(x,s^{1},\cdots,s^{N^{\rm tot}})$
$\displaystyle=\sum_{n=N_{\rm conv}+2}^{N_{\rm tot}-1}s^{{n+1}^{\top}}\cdot
w_{n}\cdot s^{n}+s^{N_{\rm conv}+1}\cdot w_{N_{\rm
conv}+1}\cdot\mathcal{F}(s_{t}^{N_{\rm conv}})$
$\displaystyle+\sum_{n=1}^{N_{\rm
conv}-1}s^{n+1}\bullet\mathcal{P}\left(w_{n+1}\star
s^{n}\right)+s^{1}\bullet\mathcal{P}\left(w_{1}\star x\right),$

when ignoring the activation function, we have:

$\forall n\in[1,N^{\rm tot}]:\quad
s_{t}^{n}\approx\frac{\partial\Phi}{\partial s^{n}}.$ (30)

Note that in the case of the Squared Error loss function, the dynamics of the
output layer derive from $\Phi$ as it can be seen by Eq. (30).

###### Learning rules for the one-sided EP estimator.

In this case, the learning rules read:

$\displaystyle\left\\{\begin{array}[]{l}\forall n\in[N_{\rm conv}+2,N_{\rm
tot}-1]:\quad\Delta w_{n}=\frac{1}{\beta}\left(s_{*}^{n+1,\beta}\cdot
s_{*}^{{n,\beta}^{\top}}-s_{*}^{n+1}\cdot s_{*}^{{n}^{\top}}\right)\\\ \Delta
w_{N_{\rm conv}+1}=\frac{1}{\beta}\left(s_{*}^{N_{\rm
conv}+1,\beta}\cdot\mathcal{F}\left(s^{N_{\rm
conv},\beta}_{*}\right)^{\top}-s_{*}^{N_{\rm
conv}+1}\cdot\mathcal{F}\left(s^{N_{\rm conv}}_{*}\right)^{\top}\right)\\\
\forall n\in[1,N_{\rm conv}-1]:\quad\Delta
w_{n+1}=\frac{1}{\beta}\left(\mathcal{P}^{-1}(s^{n+1,\beta}_{*})\star
s^{n,\beta}_{*}-\mathcal{P}^{-1}(s^{n+1}_{*})\star s^{n}_{*}\right)\\\ \Delta
w_{1}=\frac{1}{\beta}\left(\mathcal{P}^{-1}(s^{1,\beta}_{*})\star
x-\mathcal{P}^{-1}(s^{1}_{*})\star x\right)\end{array},\right.$ (35)

###### Learning rules for the symmmetric EP estimator.

In this case, the learning rules read:

$\displaystyle\left\\{\begin{array}[]{l}\forall n\in[N_{\rm conv}+2,N_{\rm
tot}-1]:\quad\Delta w_{n}=\frac{1}{2\beta}\left(s_{*}^{n+1,\beta}\cdot
s_{*}^{{n,\beta}^{\top}}-s_{*}^{n+1,-\beta}\cdot
s_{*}^{{n,-\beta}^{\top}}\right)\\\ \Delta w_{N_{\rm
conv}+1}=\frac{1}{2\beta}\left(s_{*}^{N_{\rm
conv}+1,\beta}\cdot\mathcal{F}\left(s^{N_{\rm
conv},\beta}_{*}\right)^{\top}-s_{*}^{N_{\rm
conv}+1,-\beta}\cdot\mathcal{F}\left(s^{N_{\rm
conv},-\beta}_{*}\right)^{\top}\right)\\\ \forall n\in[1,N_{\rm
conv}-1]:\quad\Delta
w_{n+1}=\frac{1}{2\beta}\left(\mathcal{P}^{-1}(s^{n+1,\beta}_{*})\star
s^{n,\beta}_{*}-\mathcal{P}^{-1}(s^{n+1,-\beta}_{*})\star
s^{n,-\beta}_{*}\right)\\\ \Delta
w_{1}=\frac{1}{2\beta}\left(\mathcal{P}^{-1}(s^{1,\beta}_{*})\star
x-\mathcal{P}^{-1}(s^{1,-\beta}_{*})\star x\right)\end{array},\right.$ (40)

#####  E.3.2 Cross-Entropy loss

###### Equations of the dynamics.

In this case, the dynamics read:

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle
s^{n+1}_{t+1}=\sigma\left(\mathcal{P}(w_{n+1}\star
s^{n}_{t})+\tilde{w}_{n+2}\star\mathcal{P}^{-1}(s^{n+2}_{t})\right),\qquad\forall
n\in[0,N^{\rm conv}-2]\\\ \displaystyle s^{N^{\rm
conv}}_{t+1}=\sigma\left(\mathcal{P}(w_{N^{\rm conv}}\star s^{N^{\rm
conv}-1}_{t})+\mathcal{F}^{-1}({w_{N^{\rm conv}+1}}^{\top}\cdot s^{N^{\rm
conv}+1}_{t})\right),\\\ \displaystyle s^{N^{\rm
conv}+1}_{t+1}=\sigma\left(w_{N^{\rm conv}+1}\cdot\mathcal{F}(s^{N^{\rm
conv}}_{t})+{w_{N^{\rm conv}+2}}^{\top}\cdot s^{N^{\rm conv}+2}_{t}\right),\\\
\displaystyle s_{t+1}^{n+1}=\sigma\left(w_{n+1}\cdot
s^{n}_{t}+{w_{n+2}}^{\top}\cdot s^{n+2}_{t}\right),\qquad\forall n\in[N^{\rm
conv}+1,N^{\rm tot}-3]\\\ \displaystyle s_{t+1}^{N^{\rm
tot}-1}=\sigma\left(w_{N^{\rm tot}-1}\cdot s^{N^{\rm
tot}-2}_{t}\right)+\beta{w_{\rm out}}^{\top}\cdot(y-\hat{y})\quad\text{with
$\beta=0$ during the first phase,}\\\ \displaystyle\hat{y}={\rm
softmax}(w_{\rm out}\cdot s^{N^{\rm tot}-1}_{t}),\end{array}\right.$ (47)

where we keep again the convention $s^{0}=x$. Considering the function:

$\displaystyle\Phi(x,s^{1},\cdots,s^{N^{\rm tot}-1})$
$\displaystyle=\sum_{n=N_{\rm conv}+1}^{N_{\rm tot}-2}s^{{n+1}^{\top}}\cdot
w_{n}\cdot s^{n}+s^{N_{\rm conv}+1}\cdot w_{N_{\rm
conv}+1}\cdot\mathcal{F}(s_{t}^{N_{\rm conv}})$
$\displaystyle+\sum_{n=1}^{N_{\rm
conv}-1}s^{n+1}\bullet\mathcal{P}\left(w_{n+1}\star
s^{n}\right)+s^{1}\bullet\mathcal{P}\left(w_{1}\star x\right),$

when ignoring the activation function, we have:

$\forall n\in[1,N^{\rm tot}-1]:\quad
s_{t}^{n}\approx\frac{\partial\Phi}{\partial s^{n}},\qquad\hat{y}={\rm
softmax}(w_{\rm out}\cdot s^{N^{\rm tot}-1}_{t}).\\\ $ (48)

Note that in this case and contrary to the Squared Error loss function, the
dynamics of the output layer do not derive from the primitive function $\Phi$,
as it can be seen from Eq. (48)

###### Learning rules for the one-sided EP estimator.

In this case, the learning rules read:

$\displaystyle\left\\{\begin{array}[]{l}\Delta w_{\rm
out}=-\left(\widehat{y}_{*}^{\beta}-y\right)\cdot s_{*}^{\beta,N^{\top}}.\\\
\forall n\in[N_{\rm conv}+2,N_{\rm tot}-2]:\quad\Delta
w_{n}=\frac{1}{\beta}\left(s_{*}^{n+1,\beta}\cdot
s_{*}^{{n,\beta}^{\top}}-s_{*}^{n+1}\cdot s_{*}^{{n}^{\top}}\right)\\\ \Delta
w_{N_{\rm conv}+1}=\frac{1}{\beta}\left(s_{*}^{N_{\rm
conv}+1,\beta}\cdot\mathcal{F}\left(s^{N_{\rm
conv},\beta}_{*}\right)^{\top}-s_{*}^{N_{\rm
conv}+1}\cdot\mathcal{F}\left(s^{N_{\rm conv}}_{*}\right)^{\top}\right)\\\
\forall n\in[1,N_{\rm conv}-1]:\quad\Delta
w_{n+1}=\frac{1}{\beta}\left(\mathcal{P}^{-1}(s^{n+1,\beta}_{*})\star
s^{n,\beta}_{*}-\mathcal{P}^{-1}(s^{n+1}_{*})\star s^{n}_{*}\right)\\\ \Delta
w_{1}=\frac{1}{\beta}\left(\mathcal{P}^{-1}(s^{1,\beta}_{*})\star
x-\mathcal{P}^{-1}(s^{1}_{*})\star x\right)\end{array}.\right.$ (54)

###### Learning rules for the symmetric EP estimator.

In this case, the learning rules read:

$\displaystyle\left\\{\begin{array}[]{l}\Delta w_{\rm
out}=-\frac{1}{2}\left(\left(\widehat{y}_{*}^{\beta}-y\right)\cdot
s_{*}^{\beta,N^{\top}}+\left(\widehat{y}_{*}^{-\beta}-y\right)\cdot
s_{*}^{-\beta,N^{\top}}\right).\\\ \forall n\in[N_{\rm conv}+2,N_{\rm
tot}-2]:\quad\Delta w_{n}=\frac{1}{2\beta}\left(s_{*}^{n+1,\beta}\cdot
s_{*}^{{n,\beta}^{\top}}-s_{*}^{n+1,-\beta}\cdot
s_{*}^{{n,-\beta}^{\top}}\right)\\\ \Delta w_{N_{\rm
conv}+1}=\frac{1}{2\beta}\left(s_{*}^{N_{\rm
conv}+1,\beta}\cdot\mathcal{F}\left(s^{N_{\rm
conv},\beta}_{*}\right)^{\top}-s_{*}^{N_{\rm
conv}+1,-\beta}\cdot\mathcal{F}\left(s^{N_{\rm
conv,-\beta}}_{*}\right)^{\top}\right)\\\ \forall n\in[1,N_{\rm
conv}-1]:\quad\Delta
w_{n+1}=\frac{1}{2\beta}\left(\mathcal{P}^{-1}(s^{n+1,\beta}_{*})\star
s^{n,\beta}_{*}-\mathcal{P}^{-1}(s^{n+1,-\beta}_{*})\star
s^{n,-\beta}_{*}\right)\\\ \Delta
w_{1}=\frac{1}{2\beta}\left(\mathcal{P}^{-1}(s^{1,\beta}_{*})\star
x-\mathcal{P}^{-1}(s^{1,-\beta}_{*})\star x\right)\end{array}.\right.$ (60)

#####  E.3.3 Implementation details in PyTorch.

The equation of the dynamics as well as the EP estimates computation can be
expressed as derivatives of the primitive function $\Phi$. Therefore, the
automatic differentiation framework provided by PyTorch can be leveraged to
implement implicitly the equations of the dynamics and the EP estimates
computation by differentiating $\Phi$. Although this implementation is slower
than explicitly implementing the equations of the dynamics, it is more
flexible in terms of network architecture as $\Phi$ is relatively easy to
compute.

####  E.4 Convolutional RNNs with unidirectional connections

In this section, we write the explicit definition of the dynamics and the
learning rule of a convolutional architecture with unidirectional connections
where forward and backward connections are no longer constrained to be equal-
valued. In this setting, we use the Cross-Entropy loss function along with a
softmax readout to implement the output layer of the network. In this setting,
the dynamics of Eq. (1) is changed into the more general form:

$s_{t+1}=F(x,s_{t},\theta),$ (61)

and the original Vector Field learning rule reads:

$\Delta\theta=\eta\widehat{\nabla}^{\rm{VF}}(\beta),\qquad\text{where}\qquad\widehat{\nabla}^{\rm{VF}}(\beta)\stackrel{{\scriptstyle\Delta}}{{=}}\frac{1}{\beta}\frac{\partial
F}{\partial\theta}(x,s_{*},\theta)^{\top}\cdot\left(s^{\beta}_{*}-s_{*}\right),$
(62)

where VF stands for Vector Field [Scellier et al., 2018]. If the transition
function $F$ derives from a primitive function $\Phi$ (i.e., if
$F=\frac{\partial\Phi}{\partial s}$), then $\widehat{\nabla}^{\rm VF}(\beta)$
is equal to $\widehat{\nabla}^{\rm EP}(\beta)$ in the limit $\beta\to 0$ (
i.e. $\lim_{\beta\to 0}\widehat{\nabla}^{\rm VF}(\beta)=\lim_{\beta\to
0}\widehat{\nabla}^{\rm EP}(\beta)$).

###### Equations of the dynamics.

In this setting, the dynamics Eq. (47) have simply to be changed into:

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle
s^{n+1}_{t+1}=\sigma\left(\mathcal{P}(w^{\rm f}_{n+1}\star
s^{n}_{t})+\tilde{w}_{n+2}^{\rm
b}\star\mathcal{P}^{-1}(s^{n+2}_{t})\right),\qquad\forall n\in[0,N^{\rm
conv}-2]\\\ \displaystyle s^{N^{\rm
conv}}_{t+1}=\sigma\left(\mathcal{P}(w^{\rm f}_{N^{\rm conv}}\star s^{N^{\rm
conv}-1}_{t})+\mathcal{F}^{-1}({w^{\rm b}_{N^{\rm conv}+1}}^{\top}\cdot
s^{1}_{t})\right),\\\ \displaystyle s^{N^{\rm
conv}+1}_{t+1}=\sigma\left(w_{N^{\rm conv}+1}^{\rm
f}\cdot\mathcal{F}(s^{N^{\rm conv}}_{t})+{w_{N^{\rm conv}+2}}^{\rm
b^{\top}}\cdot s^{N^{\rm conv}+2}_{t}\right),\\\ \displaystyle
s_{t+1}^{n+1}=\sigma\left(w^{\rm f}_{n+1}\cdot s^{n}_{t}+{w^{\rm
b}_{n+2}}^{\top}\cdot s^{n+2}_{t}\right),\qquad\forall n\in[N^{\rm
conv}+1,N^{\rm tot}-3]\\\ \displaystyle s_{t+1}^{N^{\rm
tot}-1}=\sigma\left(w^{\rm f}_{N^{\rm tot}-1}\cdot s^{N^{\rm
tot}-2}_{t}\right)+\beta{w_{\rm out}}^{\top}\cdot(y-\hat{y}),\\\
\displaystyle\hat{y}={\rm softmax}(w_{\rm out}\cdot s^{N^{\rm
tot}-1}_{t}),\end{array}\right.$ (69)

where we distinguish now between forward and backward connections: $w^{\rm
f}_{n}\neq w^{\rm b}_{n}\quad\forall n\in[1,N_{\rm tot}-2]$.

###### Original Vector Field learning rule (VF).

The symmetric version of the original Vector Field learning rule is defined
as:

$\widehat{\nabla}^{\rm{VF}}_{\rm
sym}(\beta)\stackrel{{\scriptstyle\Delta}}{{=}}\frac{1}{2\beta}\frac{\partial
F}{\partial\theta}(x,s_{*},\theta)^{\top}\cdot\left(s^{\beta}_{*}-s^{-\beta}_{*}\right),$
(70)

which yields in the case of softmax read-out:

$\displaystyle\left\\{\begin{array}[]{l}\Delta w_{\rm
out}=-\frac{1}{2}\left(\left(\widehat{y}_{*}^{\beta}-y\right)\cdot
s_{*}^{\beta,N^{\top}}+\left(\widehat{y}_{*}^{-\beta}-y\right)\cdot
s_{*}^{-\beta,N^{\top}}\right).\\\ \forall n\in[N_{\rm conv}+2,N_{\rm
tot}-2]:\quad\Delta w_{n}^{\rm
f}=\frac{1}{2\beta}\left(s_{*}^{n+1,\beta}-s_{*}^{n+1,-\beta}\right)\cdot
s_{*}^{{n}^{\top}}\\\ \forall n\in[N_{\rm conv}+2,N_{\rm tot}-2]:\quad\Delta
w_{n}^{\rm
b}=\frac{1}{2\beta}s_{*}^{{n+1}}\cdot\left(s_{*}^{n,\beta}-s_{*}^{n,-\beta}\right)^{\top}\\\
\Delta w_{N_{\rm conv}+1}^{\rm f}=\frac{1}{2\beta}\left(s_{*}^{N_{\rm
conv}+1,\beta}-s_{*}^{N_{\rm
conv}+1,-\beta}\right)\cdot\mathcal{F}\left(s^{N_{\rm
conv}}_{*}\right)^{\top}\\\ \Delta w_{N_{\rm conv}+1}^{\rm
b}=\frac{1}{2\beta}s_{*}^{N_{\rm conv}+1}\cdot\left(\mathcal{F}\left(s^{N_{\rm
conv},\beta}_{*}\right)-\mathcal{F}\left(s^{N_{\rm
conv},-\beta}_{*}\right)\right)^{\top}\\\ \forall n\in[1,N_{\rm
conv}-1]:\quad\Delta w_{n+1}^{\rm
f}=\frac{1}{2\beta}\left(\mathcal{P}^{-1}(s^{n+1,\beta}_{*})-\mathcal{P}^{-1}(s^{n+1,-\beta}_{*})\right)\star
s^{n}_{*}\\\ \forall n\in[1,N_{\rm conv}-1]:\quad\Delta w_{n+1}^{\rm
b}=\frac{1}{2\beta}\mathcal{P}^{-1}(s^{n+1}_{*})\star\left(s^{n,\beta}_{*}-s^{n,-\beta}_{*}\right)\\\
\Delta w_{1}=\frac{1}{2\beta}\left(\mathcal{P}^{-1}(s^{1,\beta}_{*})\star
x-\mathcal{P}^{-1}(s^{1,-\beta}_{*})\star x\right)\end{array}.\right.$ (79)

Importantly, note that $\Delta w^{\rm f}_{n}\neq\Delta w^{\rm
b}_{n}\quad\forall n\in[1,N_{\rm tot}-2]$.

###### Kolen-Pollack algorithm.

When forward and backward weights have a common gradient estimate, and a
weight decay term $\lambda$, they converge to the same values. We recall the
proof, noting $t$ the iteration step, $\theta_{\rm f}$ and $\theta_{\rm b}$
respectively forward and backward weights, the update rule follows:

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle\theta_{\rm
f}(t+1)=\theta_{\rm f}(t)+\Delta\theta_{\rm f}\\\ \displaystyle\theta_{\rm
b}(t+1)=\theta_{\rm b}(t)+\Delta\theta_{\rm b}\end{array}\right..$

We can then write

$\displaystyle\theta_{\rm f}(t+1)-\theta_{\rm b}(t+1)$
$\displaystyle=\theta_{\rm f}(t)-\theta_{\rm b}(t)+\Delta\theta_{\rm
f}-\Delta\theta_{\rm b}$ $\displaystyle=\theta_{\rm f}(t)-\theta_{\rm
b}(t)-\eta\lambda\left(\theta_{\rm f}(t)-\theta_{\rm b}(t)\right)$
$\displaystyle=(1-\eta\lambda)\left(\theta_{\rm f}(t)-\theta_{\rm
b}(t)\right),$

where we use the fact that the estimates are the same for both parameters,
such that they cancel out. Then by recursion :

$\displaystyle\theta_{\rm f}(t)-\theta_{\rm b}(t)$
$\displaystyle=(1-\eta\lambda)^{t}\left(\theta_{\rm f}(0)-\theta_{\rm
b}(0)\right)\underset{t\to\infty}{\rightarrow}0,\quad\text{since}\quad|1-\eta\lambda|<1.$

###### Kolen-Pollack Vector Field learning rule (KP-VF).

Assuming general dynamics of the form of Eq. (61), we distinguish forward
connections $\theta_{\rm f}$ from backward connections $\theta_{\rm b}$ so
that $\theta=\\{\theta_{\rm f},\theta_{\rm b}\\}$, with $\theta_{\rm f}$ and
$\theta_{\rm b}$ having same dimension. Assuming a first phase, a second phase
with $\beta>0$ and a third phase with $-\beta$, we define:

$\forall{\rm i\in\\{f,b\\}},\qquad\overline{\nabla_{\theta_{\rm i}}^{\rm
VF}}(\beta)=\frac{1}{2\beta}\left(\frac{\partial F}{\partial\theta_{\rm
i}}^{\top}(x,s_{*}^{\beta},\theta)\cdot s_{*}^{\beta}-\frac{\partial
F}{\partial\theta_{\rm i}}^{\top}(x,s_{*}^{-\beta},\theta)\cdot
s_{*}^{-\beta}\right)$ (80)

and we propose the following update rules:

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle\Delta\theta_{\rm
f}=\eta\left(\widehat{\nabla}^{\rm{KP-VF}}_{\rm sym}(\beta)-\lambda\theta_{\rm
f}\right)\\\ \displaystyle\Delta\theta_{\rm
b}=\eta\left(\widehat{\nabla}^{\rm{KP-VF}}_{\rm sym}(\beta)-\lambda\theta_{\rm
b}\right)\end{array}\right.,\quad\mbox{with}\quad\widehat{\nabla}^{\rm{KP-
VF}}_{\rm sym}(\beta)=\frac{1}{2}(\overline{\nabla_{\theta_{\rm f}}^{\rm
VF}}(\beta)+\overline{\nabla_{\theta_{\rm b}}^{\rm VF}}(\beta))$ (83)

where $\eta$ is the learning rate and $\lambda$ a leakage parameter. The
estimate $\widehat{\nabla}^{\rm{KP-VF}}_{\rm sym}(\beta)$ can be thought of a
generalization of Eq. (7), as highlighted in Appendix E.4 with an explicit
application of Eq. (83) to a ConvNet. More specifically, applying Eq. (80) to
Eq. (69) yields:

$\displaystyle\left\\{\begin{array}[]{l}\forall n\in[N_{\rm conv}+2,N_{\rm
tot}-2]:\\\ \qquad\qquad\qquad\overline{\nabla_{w_{n}^{\rm f}}^{\rm
VF}}(\beta)=\overline{\nabla_{w_{n}^{\rm b}}^{\rm
VF}}(\beta)=\frac{1}{2\beta}\left(s_{*}^{n+1,\beta}\cdot
s_{*}^{{n,\beta}^{\top}}-s_{*}^{n+1,-\beta}\cdot
s_{*}^{{n,-\beta}^{\top}}\right)\\\ \overline{\nabla_{w_{N_{\rm conv}+1}^{\rm
f}}^{\rm VF}}(\beta)=\overline{\nabla_{w_{N_{\rm conv}+1}^{\rm b}}^{\rm
VF}}(\beta)=\\\ \qquad\qquad\qquad\frac{1}{2\beta}\left(s_{*}^{N_{\rm
conv}+1,\beta}\cdot\mathcal{F}\left(s^{N_{\rm
conv},\beta}_{*}\right)^{\top}-s_{*}^{N_{\rm
conv}+1,-\beta}\cdot\mathcal{F}\left(s^{N_{\rm
conv},-\beta}_{*}\right)^{\top}\right)\\\ \forall n\in[1,N_{\rm conv}-1]:\\\
\overline{\nabla_{w_{n+1}^{\rm f}}^{\rm
VF}}(\beta)=\frac{1}{2\beta}\left(\mathcal{P}^{-1}\left(s^{n+1,\beta}_{*},{\rm
ind}_{\mathcal{P}}\left(w_{n+1}^{\rm f}\star
s^{n,\beta}_{*}\right)\right)\star s^{n,\beta}_{*}\right.\\\
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\left.-\mathcal{P}^{-1}\left(s^{n+1,-\beta}_{*},{\rm
ind}_{\mathcal{P}}\left(w_{n+1}^{\rm f}\star
s^{n,-\beta}_{*}\right)\right)\star s^{n,-\beta}_{*}\right)\\\ \forall
n\in[1,N_{\rm conv}-1]:\\\ \overline{\nabla_{w_{n+1}^{\rm b}}^{\rm
VF}}(\beta)=\frac{1}{2\beta}\left(\mathcal{P}^{-1}\left(s^{n+1,\beta}_{*},{\rm
ind}_{\mathcal{P}}\left(w_{n+1}^{\rm b}\star
s^{n,\beta}_{*}\right)\right)\star s^{n,\beta}_{*}\right.\\\
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\left.-\mathcal{P}^{-1}\left(s^{n+1,-\beta}_{*},{\rm
ind}_{\mathcal{P}}\left(w_{n+1}^{\rm b}\star
s^{n,-\beta}_{*}\right)\right)\star s^{n,-\beta}_{*}\right)\\\
\end{array}.\right.$ (94)

Combining Eqs. (94) with Eq. (83) gives the associated parameter updates. The
updates for $w_{1}$ and $w_{\rm out}$ are the same than those of Eq. (79).
Importantly, note that while $\forall n\in[N_{\rm conv}+1,N_{\rm
tot}-2]:\quad\overline{\nabla_{w_{n}^{\rm f}}^{\rm
VF}}(\beta)=\overline{\nabla_{w_{n}^{\rm b}}^{\rm VF}}(\beta)$, we have
$\forall n\in[1,N_{\rm conv}-1]:\quad\overline{\nabla_{w_{n}^{\rm f}}^{\rm
VF}}(\beta)\neq\overline{\nabla_{w_{n}^{\rm b}}^{\rm VF}}(\beta)$ because of
inverse pooling. In other words, the updates of the convolutional filters do
not solely depend on the pre and post synaptic activations but also on the
location of the maximal elements within each pooling window, itself depending
on the filter considered. Hence the motivation to average
$\overline{\nabla_{w_{n}^{\rm f}}^{\rm VF}}(\beta)$ and
$\overline{\nabla_{w_{n}^{\rm b}}^{\rm VF}}(\beta)$ and use this quantity to
update to $w_{n}^{\rm b}$ and $w_{n}^{\rm f}$ and apply the Kolen-Pollack
technique.

###### Implementation details in PyTorch.

The dynamics in the case of unidirectional connections does not derive from a
primitive function $\Phi$. Therefore, it is not possible to implicitly get the
dynamics by differentiating one primitive function. A way around is to get the
VF dynamics by differentiating one quantity $\tilde{\Phi}^{n}$ by layer. This
quantity is not a primitive function and is especially designed to get the
right equations once differentiated. We define $\tilde{\Phi}^{n}(w_{n}^{\rm
f},w_{n+1}^{\rm b},s^{n-1},s^{n})$ by:

$\displaystyle\left\\{\begin{array}[]{l}\forall n\in[1,N_{\rm
conv}-1]:\tilde{\Phi}^{n}=s^{n}\bullet\mathcal{P}(w_{n}^{\rm f}\star
s^{n-1})+s^{n+1}\bullet\mathcal{P}(w_{n+1}^{\rm b}\star s^{n})\\\
\tilde{\Phi}^{N_{\rm conv}}=s^{N_{\rm conv}}\bullet\mathcal{P}(w_{N_{\rm
conv}}^{\rm f}\star s^{N_{\rm conv}-1})+s^{N_{\rm conv}+1}\cdot w_{N_{\rm
conv}+1}^{\rm b}\cdot\mathcal{F}(s^{N_{\rm conv}})\\\ \forall n\in[N_{\rm
conv}+1,N_{\rm tot}-1]:\tilde{\Phi}^{n}=s^{n}\cdot w_{n}^{\rm f}\cdot
s^{n-1}+s^{n+1}\cdot w_{n+1}^{\rm b}\cdot s^{n}\\\ \tilde{\Phi}^{N_{\rm
tot}-1}=s^{N_{\rm tot}-1}\cdot w_{N_{\rm tot}-1}^{\rm f}\cdot s^{N_{\rm
tot}-2}+\beta\ell(s^{N_{\rm tot}-1},y,w_{\rm out})\end{array},\right.$ (99)

where $\ell$ is defined by Eq. (15), and $\beta=0$ in the first phase. Then,
$\forall n\in[1,N^{\rm tot}-1]$, the dynamics of Eq. (69) read:

$s^{n}_{t+1}=\sigma\left(\frac{\partial\tilde{\Phi}^{n}}{\partial
s^{n}}(w_{n}^{\rm f},w_{n+1}^{\rm b},s^{n-1}_{t},s^{n}_{t})\right).$ (100)

The original VF update of Eq. (79) can be written as $\forall n\in[1,N^{\rm
tot}-1],\forall{\rm i}\in\\{{\rm f},{\rm b}\\}$:

$\Delta w_{n}^{\rm
i}=\frac{1}{2\beta}\left(\frac{\partial\tilde{\Phi}^{n}}{\partial w_{n}^{\rm
i}}(s^{n,\beta}_{*},s^{n-1}_{*})-\frac{\partial\tilde{\Phi}^{n}}{\partial
w_{n}^{\rm i}}(s^{n,-\beta}_{*},s^{n-1}_{*})\right),$ (101)

and Eq. (94) as:

$\overline{\nabla_{w_{n}^{\rm i}}^{\rm
VF}}(\beta)=\frac{1}{2\beta}\left(\frac{\partial\tilde{\Phi}^{n}}{\partial
w_{n}^{\rm
i}}(s^{n,\beta}_{*},s^{n-1,\beta}_{*})-\frac{\partial\tilde{\Phi}^{n}}{\partial
w_{n}^{\rm i}}(s^{n,-\beta}_{*},s^{n-1,-\beta}_{*})\right).$ (102)

###  Appendix F Experimental details

Figure 3: Train (dashed) and test (solid) errors on CIFAR-10 with the Squared
Error loss function. The curves are averaged over 5 runs and shadows stand for
$\pm 1$ standard deviation. Figure 4: Train (dashed) and test (solid) errors
on CIFAR-10 with the Cross-Entropy loss function. The curves are averaged over
5 runs and shadows stand for $\pm 1$ standard deviation.

####  F.1 Environment and hyper-parameters

The experiments are run using PyTorch 1.4.0 and torchvision 0.5.0. [Paszke et
al., 2017]. The convolutional architecture used in the CIFAR-10 experiment
consists of four $3\times 3$ convolutional layers of respective feature maps
128 - 256 - 512 - 512. We use a stride of one for each convolutional layer,
and zero-padding of one for each layer except for the last layer. Each layer
is followed by a $2\times 2$ Max Pooling operation with a stride of two. The
resulting flattened feature vector is of size 512. The weights are initialized
using the default initialization of PyTorch, which is the uniform Kaiming
initialization introduced by He et al. [2015]. The data is normalized and
augmented with random horizontal flips and random crops. The training is
performed with stochastic gradient descent with momentum and weight decay. We
use the learning rate scheduler introduced by Loshchilov and Hutter [2016] to
speed up convergence. The simulations were carried across several servers
consisting of 14 GPUs in total. Each run was performed on a single GPU for an
average run time of 2 days.

For the unidirectional architecture, the backward weights are defined for all
convolutional layers except the first convolutional layer connected to the
static input. The forward and backward weights are initialized independently
at the beginning of training. The backward weights have no bias contrary to
their forward counterparts. The hyper-parameters such as learning rate, weight
decay and momentum are shared between forward and backward weights.

Table 2: Hyper-parameters used for the CIFAR-10 experiments. Hyper-parameter | Squared Error | Cross-Entropy
---|---|---
$T$ | 250 | 250
$K$ | 30 | 25
$\beta$ | 0.5 | 1.0
Batch Size | 128 | 128
| Initial learning rates
---
(Layer-wise)
0.25 - 0.15 - 0.1 - 0.08 - 0.05 | 0.25 - 0.15 - 0.1 - 0.08 - 0.05
Final learning rates | $10^{-5}$ | $10^{-5}$
| Weight decay
---
(All layers)
$3\cdot 10^{-4}$ | $3\cdot 10^{-4}$
Momentum | 0.9 | 0.9
Epoch | 120 | 120
| Cosine Annealing
---
Decay time (epochs)
100 | 100

####  F.2 Random-sign estimate variance

Although not explicitly stated in this purpose, the use of such randomization
has been reported in some earlier publications on the MNIST task [Scellier and
Bengio, 2017, Ernoult et al., 2020]. However, in this work, we show that this
method exhibits high variance in the training procedure. The results presented
in Table 1 consists of five runs. In the case of the EP random-sign estimate,
one run among the five collapses to random guess similar to the one-sided
estimate. In order to test the frequency of such a phenomenon, we performed
another five runs with both symmetric and random-sign estimates. The results
for each run presented in Table 3 show that two trials among ten are unstable,
confirming further the high variance nature of the random-sign estimate.

Figure 5: Test error curve of each run with the Squared Error loss function and random-sign estimate. The two collapsed runs among the ten trials are steady to 90% because in such cases the network typically outputs the same class for each data point. Table 3: Best test error comparison between random-sign and symmetric estimates, for ten runs. Run index | EP random-sign | EP symmetric
---|---|---
$1$ | $12.97$ | $12.24$
$2$ | $12.72$ | $12.31$
$3$ | $12.30$ | $12.68$
$4$ | $12.45$ | $12.43$
$5$ | $12.78$ | $12.57$
$6$ | $12.66$ | $12.55$
$7$ | $12.84$ | $12.44$
$8$ | $12.59$ | $12.52$
$9$ | $57.32$ | $12.85$
$10$ | $89.98$ | $12.60$
Mean | $24.86$ | $\mathbf{12.52}$
w/o collapse | $\mathbf{12.66}$ | N.A

####  F.3 Adding dropout

We adapt dropout [Srivastava et al., 2014] for convergent RNNs by shutting
some units to zero with probability $p<1$ when computing
$\Phi(x,s_{t},\theta)$. We multiply the remaining active units by the factor
$\frac{1}{1-p}$ to keep the same neural activity on average, so that the
learning rule is rescaled by $\left(\frac{1}{1-p}\right)^{2}$. The dropped out
units are the same within one training iteration but they differ across the
examples of one mini batch. In our experiments we use $p=0.1$ on the last
convolutional layer before the linear classifier. The results are reported in
Table 1.

####  F.4 Changing the activation function

Previous implementations of EP used a shifted hard sigmoid activation
function:

$\sigma(x)=\max(0,\min(x,1)).$ (103)

In their experiments with ConvNets on MNIST, Ernoult et al. [2019] observed
saturating units that cannot pass error signals during the second phase. In
this work, to mitigate this effect, we have rescaled by a factor $1/2$ the
slope of the activation function to ease signal propagation and prevent
saturation, therefore changing Eq. (103) into:

$\sigma(x)=\max\left(0,\min\left(\frac{x}{2},1\right)\right).$ (104)

###  Appendix G Weight alignment for unidirectional connections

The angle $\alpha$ between forward and backward weights is defined as :

$\alpha=\frac{180}{\pi}{\rm Acos}\left(\frac{w^{\rm b}\bullet w^{\rm
f}}{\|w^{\rm b}\|\|w^{\rm
f}\|}\right),\quad\text{where}\quad\|w\|=\sqrt{w\bullet w}.$ (105)

Fig. 6 shows the angle between forward and backward weights during training on
CIFAR-10 for both the original VF learning rule (dashed) and the new learning
rule inspired by Kolen and Pollack [1994].

Figure 6: Angle $\alpha$ between forward and backward weights for the new
estimate $\widehat{\nabla}^{\rm{KP-VF}}_{\rm sym}$ introduced (solid) and
$\widehat{\nabla}^{\rm{VF}}_{\rm sym}$ (dashed).

###  Appendix H Layer-wise comparison of EP estimates

In this section we show on Fig. 7 more instances of Fig. 1(d) for each layer
of the convolutional architecture.

Figure 7: Layer-wise comparison between EP gradient estimates and BPTT
gradients for 5 layers deep CNN on CIFAR-10 Data. Layer index increases from
top to bottom, left to right, top-left being the first layer.