# Improving equilibrium propagation without weight symmetry through Jacobian
homeostasis

Axel Laborieux1 &Friedemann Zenke1,2<EMAIL_ADDRESS>
1 Friedrich Miescher Institute for Biomedical Research, Basel, Switzerland
2 Faculty of Natural Sciences, University of Basel, Basel, Switzerland

###### Abstract

EP is a compelling alternative to the BP for computing gradients of neural
networks on biological or analog neuromorphic substrates. Still, the algorithm
requires weight symmetry and infinitesimal equilibrium perturbations, i.e.,
nudges, to estimate unbiased gradients efficiently. Both requirements are
challenging to implement in physical systems. Yet, whether and how weight
asymmetry affects its applicability is unknown because, in practice, it may be
masked by biases introduced through the finite nudge. To address this
question, we study generalized equilibrium propagation (EP), which can be
formulated without weight symmetry, and analytically isolate the two sources
of bias. For complex-differentiable non-symmetric networks, we show that the
finite nudge does not pose a problem, as exact derivatives can still be
estimated via a Cauchy integral. In contrast, weight asymmetry introduces bias
resulting in low task performance due to poor alignment of EP’s neuronal error
vectors compared to BP. To mitigate this issue, we present a new homeostatic
objective that directly penalizes functional asymmetries of the Jacobian at
the network’s fixed point. This homeostatic objective dramatically improves
the network’s ability to solve complex tasks such as ImageNet 32$\times$32\.
Our results lay the theoretical groundwork for studying and mitigating the
adverse effects of imperfections of physical networks on learning algorithms
that rely on the substrate’s relaxation dynamics.

##  1 Introduction

Virtually all state-of-the-art artificial intelligence models are trained
using the backpropagation of error algorithm (BP) [1]. In BP, neural
activations are evaluated during the forward pass, and gradients are obtained
efficiently through a corresponding backward pass. While BP can be implemented
efficiently on digital hardware which _simulates_ neural networks, it is less
suitable for computing gradients in physical neural networks such as brains or
analog neuromorphic substrates. This limitation is mainly due to two technical
requirements of the backward pass in BP. First, while the forward pass
consists of a nonlinear propagation of neural activation, the backward pass is
a linear operation as it amounts to solving a linear system to the erros, and
these errors do not affect the neuronal activity [2, 3]. Second, back-
propagation requires the transpose of connectivity matrices, which is also
often referred to as weight symmetry or “weight transport” [4]. In physical
networks both linear signal propagation and weight symmetry are difficult to
achieve. For instance, the biochemical signaling in the brain is intrinsically
stochastic and nonlinear.

EP [5] is an alternative algorithm for rigorously computing loss gradients
using _only_ the network’s dynamics, unlike BP. To accomplish this feat,
classic EP requires an energy-based model (EBM), such as a differentiable
Hopfield network [6], which relaxes to an equilibrium for any given input. EP
computes gradients by comparing neuronal activity at the free equilibrium to
the activity at a second equilibrium that is “nudged” toward a desired target
output. Unlike BP in feed-forward networks, EP uses the _same_ network
dynamics to propagate errors and feed-forward activity. EP _implicitly_
computes gradients without the need to linearly propagate error vectors, and
could thus operate on analog and biological neural networks.

Nevertheless, classic EP still requires an infinitesimally small nudge and
weight symmetry to avoid unwanted biases in the gradient estimates. These
requirements impede an interpretation in the biological context and imposes
additional constraints for analog neuromorphic hardware. While bias from a
large nudge can be, in principle, offset by dint of holomorphic EP [7], a
complex valued generalization, this approach has only been demonstrated using
symmetric weights. Although a generalization of EP to non-symmetric dynamical
systems exists [8], the approach has only been successfully applied to
relatively simple tasks like MNIST [8, 9, 10]. In particular, generalized EP
fails to train on CIFAR-10, even with an infinitesimal nudge, suggesting that
weight asymmetry induces bias [11]. More comprehensive studies on the effect
of non-symmetric feedback connections on training have only been done in the
context of BP [12, 13, 14, 15, 16]. While alternative approaches studied
learning beyond approximate BP [17, 18, 19, 20, 21], here we focus on
generalized EP due to its strong theoretical foundation and because it enables
physical systems to learn with their dynamics.

In summary, it remains an open question whether and how generalized EP can
compute suitable gradient estimates for finite nudges and asymmetric weights
and how well the resulting networks perform on complex real-world tasks for
varying degrees of weight asymmetry. Overall, the main contributions of this
work are:

  * •

A comprehensive analysis of the individual sources of bias in the gradient
estimate from weight asymmetry and finite-size nudge in generalized EP.

  * •

An extension of the integral estimate of holomorphic EP [7] to non symmetric
complex-differentiable dynamical systems, and show that it can be estimated in
continuous time.

  * •

A new homeostatic loss that reduces the asymmetry of the Jacobian at the free
equilibrium point without enforcing perfect weight symmetry.

  * •

An empirical demonstration that generalized EP with weight asymmetry and
homeostatic loss scale to ImageNet 32$\times$32, with a small performance gap
compared to the symmetric case.

##  2 Background

Before we start our analysis, we briefly review how generalized EP computes
gradients in dynamical systems with fixed points. Further we will introduce
neuronal error vectors as a useful notion on which we rest our analysis.

###  2.1 Gradients in converging dynamical systems

Let us first recall how gradients are computed in converging dynamical systems
[22]. To that end, let us consider a dynamical system with state variable
$\boldsymbol{u}$ and parameters $\boldsymbol{\theta}$. A subset of units
receive static input currents $\boldsymbol{x}$. The dynamics of the system are
governed by a differentiable vector field $F$:

$\frac{\mathrm{d}\boldsymbol{u}}{\mathrm{d}t}=F(\boldsymbol{\theta},\boldsymbol{u},\boldsymbol{x}).$
(1)

We further assume that the dynamics of Eq. (1) converge to a fixed point of
the state $\boldsymbol{u}$. This phase of the dynamics is defined as the “free
phase”, because there are no contributions from any target $\boldsymbol{y}$
associated to $\boldsymbol{x}$. While stable fixed points provably exists for
EBMs [6], there is no guarantee to converge for arbitrary dynamical systems.
However, a fixed point exists if $F$ is a contraction map [23, 24]. In the
following, we just assume the existence of a fixed point and denote it by
$\boldsymbol{u}^{\ast}_{0}$:

$0=F(\boldsymbol{\theta},\boldsymbol{u}^{\ast}_{0},\boldsymbol{x}).$ (2)

In $\boldsymbol{u}^{\ast}_{0}$, the superscript $\ast$ stands for the
convergence to a fixed point, and the subscript $0$ stands for the “free”
phase, but its meaning will be rigorously defined in the next subsection. We
further assume that $\boldsymbol{u}^{\ast}_{0}$ is an implicit function of
$\boldsymbol{\theta}$, such that its derivative exist by means of the implicit
function theorem (IFT). Suppose now that an objective function $\mathcal{L}$
measures the proximity of a subset of $\boldsymbol{u}^{\ast}$ (output units)
to a target $\boldsymbol{y}$. The gradient of $\mathcal{L}$ with respect to
$\boldsymbol{\theta}$ is obtained by the chain rule through the fixed point
$\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{\theta}}=\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}\frac{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}{\mathrm{d}\boldsymbol{\theta}}$,
where $\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}$ is
the error at the output units, and
$\frac{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}{\mathrm{d}\boldsymbol{\theta}}$ is
obtained by differentiating Eq. (2) with respect to $\boldsymbol{\theta}$:

$\displaystyle 0$ $\displaystyle=\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0})+J_{F}(\boldsymbol{u}^{\ast}_{0})\cdot\frac{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}{\mathrm{d}\boldsymbol{\theta}},$
$\displaystyle\frac{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}{\mathrm{d}\boldsymbol{\theta}}$
$\displaystyle=-J_{F}(\boldsymbol{u}^{\ast}_{0})^{-1}\cdot\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0}),$ (3)

Where $J_{F}(\boldsymbol{u}^{\ast}_{0}):=\frac{\partial
F}{\partial\boldsymbol{u}}(\boldsymbol{\theta},\boldsymbol{u}^{\ast}_{0},\boldsymbol{x})$
is the Jacobian of the network at the fixed point (Fig. 1). Replacing
$\frac{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}{\mathrm{d}\boldsymbol{\theta}}$ in
the expression of the gradient
$\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{\theta}}$ and transposing
the whole expression yields [22]:

$\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{\theta}}^{\top}=\underbrace{\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0})^{\top}}_{\text{``pre-
synaptic"}}\cdot\underbrace{\left(-J_{F}(\boldsymbol{u}^{\ast}_{0})^{-\top}\cdot\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}^{\top}\right)}_{\text{``post-
synaptic"}~{}\boldsymbol{\delta}}.$ (4)

We refer to the part of the expression in parenthesis as the neuronal error
vector $\boldsymbol{\delta}$. In the widely used reverse-mode automatic
differentiation (AD), a linear system involving the transpose of the Jacobian
is solved to compute $\boldsymbol{\delta}$ from the output error
$\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}$, and
finally multiply with $\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0})$, which is a sparse
matrix involving pre-synaptic variables. In feed forward networks, the lower-
triangular structure of the Jacobian allows computing $\boldsymbol{\delta}$
recursively layer by layer from the output (Fig. 1c) as exploited in BP [1].
In converging dynamical systems, recurrent backpropagation (RBP) and variants
[25, 26, 24] obtain $\boldsymbol{\delta}$ as the fixed point of an auxiliary
linear dynamical system using $J_{F}(\boldsymbol{u}^{\ast}_{0})^{\top}$. In
deep equilibrium models [27] $\boldsymbol{\delta}$ is found with a root-
finding algorithm. In the next section, we review how $\boldsymbol{\delta}$ is
obtained in EP.

Figure 1: Different neural network architectures and their Jacobian after the
free phase. a) A continuous layered Hopfield network. The existence of an
energy function enforces Jacobian symmetry with reciprocal and tied
connectivity $w_{ij}=w_{ji}$. b) A layered network with reciprocal
connectivity but independent forward and backward connections, which make the
Jacobian not symmetric. c) A discrete feed forward network without feedback
connections. The Jacobian is lower triangular, allowing the explicit back
propagation of error layer by layer.

###  2.2 Neuronal errors in Equilibrium Propagation

The EP [5] gradient formula for EBMs can be derived without explicit use of
the energy function formalism [8], by noticing a similar pattern between the
derivative of $\boldsymbol{u}^{\ast}_{0}$ with respect to parameters
$\boldsymbol{\theta}$ (Eq. (3)) and $\boldsymbol{\delta}$ in Eq. (4). While
$\frac{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}{\mathrm{d}\boldsymbol{\theta}}$ is
obtained by inversion of $J_{F}(\boldsymbol{u}^{\ast}_{0})$,
$\boldsymbol{\delta}$ is obtained by inversion of
$J_{F}(\boldsymbol{u}^{\ast}_{0})^{\top}$. However, in EBMs for which there
exists a scalar energy function $E$ such that $F=-\frac{\partial
E}{\partial\boldsymbol{u}}$, we have precisely
$J_{F}(\boldsymbol{u}^{\ast}_{0})=J_{F}(\boldsymbol{u}^{\ast}_{0})^{\top}$
(Fig. 1a) due to the Schwarz theorem on symmetric second derivatives. The
remaining difference between Eqs. (3) and (4) is the quantity to which the
inverted Jacobian is applied: $\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0})$ in Eq. (3) and
$\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}^{\top}$ in
Eq. (4). The solution is to define a new scalar parameter $\beta$ such that
these quantities are the same:

######  Definition 1 (Nudge parameter $\beta$).

Let $\beta$ be a scalar parameter that is equal to zero during the free phase,
justifying the notation $\boldsymbol{u}^{\ast}_{0}$, and such that
$\frac{\partial
F}{\partial\beta}(\boldsymbol{u}^{\ast}_{0})=\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}^{\top}$.
When $\beta\neq 0$, the target $\boldsymbol{y}$ contributes to the dynamics
through $\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}$.
The term “teaching amplitude” refers to $|\beta|$.

By applying the same derivation as Eq. (3) for the new parameter $\beta$, we
have:

$\left.\frac{\mathrm{d}\boldsymbol{u}^{\ast}}{\mathrm{d}\beta}\right|_{\beta=0}=-J_{F}(\boldsymbol{u}^{\ast}_{0})^{-1}\cdot\frac{\partial
F}{\partial\beta}(\boldsymbol{u}^{\ast}_{0})=-\underbrace{J_{F}(\boldsymbol{u}^{\ast}_{0})^{-\top}}_{\text{for
\acsp{EBM}}}\cdot\underbrace{\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}^{\top}}_{\text{by
Def.~{}\ref{def:beta}}}=\boldsymbol{\delta}.$ (5)

Therefore, the EP neuronal error vector is the derivative of the fixed point
$\boldsymbol{u}^{\ast}_{0}$ with respect to $\beta$ in $\beta=0$, which we
denote by $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ for short. Note that the
derivation assumes that $\boldsymbol{u}^{\ast}_{0}$ is an implicit function of
$\beta$. In practice, $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ is
approximated by finite differences after letting the network settle to a
second equilibrium
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}\approx(\boldsymbol{u}^{\ast}_{\beta}-\boldsymbol{u}^{\ast}_{0})/\beta$,
which contains a bias due to the finite $\beta\neq 0$ used in the second
phase. The EP gradient estimate can then be written without using the energy
function by replacing $\boldsymbol{\delta}$ by
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ in Eq. (4):

$\tilde{\nabla}_{\boldsymbol{\theta}}:=\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0})^{\top}\cdot\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}.$
(6)

Importantly, the quantities $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ and
$\tilde{\nabla}_{\boldsymbol{\theta}}$ are defined for _any_ differentiable
dynamical system [8], not just for EBMs, but they coincide with
$\boldsymbol{\delta}$ and the actual gradient
$\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{\theta}}$ (Eq. (4)) only
when the Jacobian $J_{F}(\boldsymbol{u}^{\ast}_{0})$ is symmetric [28], which
applies only to EBMs. We provide more details on the relation between Eq. (6)
and the usual gradient formula involving the energy function [5] in Appendix
B. Intuitively, $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ can be seen as the
replacement for $\boldsymbol{\delta}$ that the system can obtain using its own
dynamics. When $J_{F}(\boldsymbol{u}^{\ast}_{0})\neq
J_{F}(\boldsymbol{u}^{\ast}_{0})^{\top}$ (Fig. 1b), the asymmetric part of the
Jacobian contributes an additional bias to the estimate of
$\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{\theta}}$. However, its
contribution has never been studied independently because it is masked by the
finite nudge $\beta$ in practice. It is worth nothing that
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}=0$ for feed forward networks due to
the absence of feedback connections (Fig. 1c).

##  3 Theoretical results

###  3.1 Generalized holomorphic EP

To study individual bias contributions we build on Holomorphic EP [7], which,
in the case of weight symmetry, allows computing unbiased gradients with
finite nudges. As we will see, the property generalizes to asymmetric weights
and $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ can be exactly computed with
finite teaching amplitude $|\beta|>0$ for any complex differentiable
(holomorphic) vector field $F$. We use the fact that for holomorphic
functions, derivatives at a given point are equal to the Cauchy integral over
a path around the point of interest $\beta=0$ [29] (see Appendix A.1 for more
details):

$\left.\frac{\mathrm{d}\boldsymbol{u}^{\ast}}{\mathrm{d}\beta}\right|_{\beta=0}=\frac{1}{2\mathrm{i}\pi}\oint_{\gamma}\frac{\boldsymbol{u}^{\ast}_{\beta}}{\beta^{2}}\mathrm{d}\beta=\frac{1}{T|\beta|}\int_{0}^{T}\boldsymbol{u}^{\ast}_{\beta(t)}e^{-2\mathrm{i}\pi
t/T}\mathrm{d}t.$ (7)

Here we used the teaching signal $\beta(t)=|\beta|e^{2\mathrm{i}\pi
t/T},t\in[0,T]$ for the path $\gamma$, and performed the change of variable
$\mathrm{d}\beta=\beta(t)\frac{2\mathrm{i}\pi}{T}\mathrm{d}t$, where
$\mathrm{i}\in\mathbb{C}$ is the imaginary unit. Moreover, Eq. (7) assumes
that the fixed point $\boldsymbol{u}^{\ast}$ is an implicit function of
$\beta$ and well-defined on the entire path $\gamma$. In other words, we
assume that the system is at equilibrium for all $t\in[0,T]$, which can be
achieved when the timescale of the dynamical system is shorter than the one of
the time-dependent nudge $\beta(t)$. In numerical simulations, we estimate the
integral with a discrete number of $N$ points:
$\nicefrac{{t}}{{T}}=\nicefrac{{k}}{{N}},k\in[0,...,N-1]$ (see Appendix D.1).
Regardless of whether the dynamical system has a symmetric Jacobian at the
free fixed point or not (Fig. 1), the exact value of
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ can be computed for a finite range
of teaching amplitudes $|\beta|$ (Fig. 2a,b).

###  3.2 Estimating gradients in continuous time

Although Eq. (7) implies that $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ is
obtained at the end of one period of the teaching signal ($t=T$), it can also
estimated continuously by letting the integral run over multiple oscillation
cycles of the teaching signal $\beta(t)$ (Fig. 2 f):

$\frac{1}{t|\beta|}\int_{0}^{t}\boldsymbol{u}^{\ast}_{\beta(\tau)}e^{-2\mathrm{i}\pi\tau/T}\mathrm{d}\tau\underset{t\to\infty}{\longrightarrow}\left.\frac{\mathrm{d}\boldsymbol{u}^{\ast}}{\mathrm{d}\beta}\right|_{\beta=0}.$
(8)

This relation is made obvious by splitting the integral over all the completed
periods, such that the relative contribution of the remainder integral
asymptotically approaches $0$ as $t\to\infty$ (see Appendix C.1). The
advantage of this approach for physical systems is that
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ need not be accessed at the precise
times $t=T$, thereby dispensing with the need for separate phases [30].

Moreover, the free fixed point does not need to be evaluated to obtain
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$. However, the pre-synaptic term
$\frac{\partial F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0})$ in
Eq. (6) a priori still requires evaluating the free fixed point. We will see
next that we can also estimate it continuously in time thanks to the Mean
Value Theorem from complex analysis [29] (Fig. 2e), Appendix C.1):

$\frac{1}{t}\int_{0}^{t}\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{\beta(\tau)})\mathrm{d}\tau\underset{t\to\infty}{\longrightarrow}\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0}).$ (9)

Overall, the gradient estimate $\tilde{\nabla}_{\boldsymbol{\theta}}$ in the
time-continuous form corresponds to the average of pre-synaptic activation
multiplied by the amplitude of the post-synaptic potential, consistent with
models of STDP [31, 32]. In the next section, we isolate the contribution to
the bias brought by the asymmetry of the Jacobian
$J_{F}(\boldsymbol{u}^{\ast}_{0})$.

Figure 2: Separating the sources of bias due to finite nudge and Jacobian
asymmetry. a,b) Cosine similarity between the neuronal error vector
$\boldsymbol{\delta}$ computed by RBP and classic EP (grey) [5, 8],
holomorphic EP (purple and pink) for different number of points to estimate
Eq. (7), in function of the teaching amplitude $|\beta|$, for symmetric (a))
and non symmetric (b)) equilibrium Jacobian. c,d) Kernel density estimate of
the cosine distribution between both neuronal error vectors in a two-hidden
layer MLP. Alignment worsens with increasing asymmetry and depth. e)
Continuous-time estimate of the free fixed point through oscillations induced
by the teaching signal $\beta(t)$. f) Continuous-time estimate of the neuronal
error vector of generalized hEP.

###  3.3 Isolating bias from Jacobian asymmetry

We are now in the position to quantify the bias originating from asymmetric
Jacobians in the computation of the neuronal error vector of hEP. By isolating
$\frac{\partial F}{\partial\beta}$ and
$\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}\boldsymbol{u}^{\ast}_{0}}$ on both
ends of Eq. (5), and because they are equal by Definiton 1, we see that both
the neuronal error vectors of hEP and RBP are linked by the following
relation:

$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}=J_{F}(\boldsymbol{u}^{\ast}_{0})^{-1}J_{F}(\boldsymbol{u}^{\ast}_{0})^{\top}\boldsymbol{\delta},$
(10)

Thus, for asymmetric Jacobians, $\boldsymbol{\delta}$ and
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ are not aligned which introduces
bias in the resulting gradient estimates.

To further quantify this bias, we introduce $S$ and $A$ the symmetric and
skew-symmetric parts of $J_{F}$, such that
$J_{F}(\boldsymbol{u}^{\ast}_{0})=S+A$ and
$J_{F}(\boldsymbol{u}^{\ast}_{0})^{\top}=S-A$. Then we can show that when the
Frobenius norm
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|S^{-1}A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\to
0$, $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ and $\boldsymbol{\delta}$ are
linked by (see Appendix C.2 for a proof):

$\displaystyle\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}=\boldsymbol{\delta}-2S^{-1}A\boldsymbol{\delta}+o(S^{-1}A\boldsymbol{\delta}).$
(11)

This expression highlights the direct dependence of the bias on $A$ and
further suggests a way to improve the Jacobian’s symmetry, thus reducing the
bias by directly decreasing the norm of $A$.

###  3.4 Increasing functional symmetry through Jacobian homeostasis.

To reduce the norm of $A$ we note that
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}=\operatorname{Tr}(A^{\top}A)=\tfrac{1}{2}\operatorname{Tr}(J_{F}^{\top}J_{F})-\tfrac{1}{2}\operatorname{Tr}(J_{F}^{2})$,
where the dependence of $J_{F}$ on $\boldsymbol{u}^{\ast}$ is omitted. Similar
to work by Bai et al. [33] on deep equilibrium models, we can use the
classical Hutchinson trace estimator [34] (see Appendix A.2) to obtain the
following objective for minimizing Jacobian asymmetry and, therefore,
improving “functional symmetry”:

$\mathcal{L}_{\text{homeo}}:=\underset{\boldsymbol{\varepsilon}\sim\mathcal{N}(0,I)}{\mathbb{E}}\left[\|J_{F}\boldsymbol{\varepsilon}\|^{2}-\boldsymbol{\varepsilon}^{\top}J_{F}^{2}\boldsymbol{\varepsilon}\right].$
(12)

If we break down the two terms present in this objective, we see that
$\|J_{F}\boldsymbol{\varepsilon}\|^{2}$ should be minimized, which corresponds
to making the network more robust to random perturbations by minimizing the
Frobenius norm of the Jacobian. Interestingly, Stock et al. [35] showed how a
local heterosynaptic learning rule could optimize this objective locally
through synaptic balancing (see Eqs. (12)–(14) in [35]), whereby neurons try
to balance their net inputs with their output activity. The second term in Eq.
(12) comes with a minus sign and should therefore be maximized. It is, thus,
akin to a reconstruction loss, similar to the training objective of feedback
connections in target propagation (TP) [36, 37, 38, 39, 40], and reminiscent
of predictive slow feature learning [41, 42, 43]. Overall,
$\mathcal{L}_{\text{homeo}}$ is a plausible objective for increasing and
maintaining functional symmetry of the Jacobian in any dynamical system.
Finally, it is important to note that while perfect weight symmetry implies
functional symmetry, the converse is not true. For instance, Eq. (12) can be
optimized in arbitrary systems without reciprocal connectivity. This
distinction makes $\mathcal{L}_{\text{homeo}}$ more general than alternative
strategies acting on weight symmetry directly [44, 45].

##  4 Experiments

In the following experiments, we used the setting of convergent recurrent
neural networks [46], as well as the linear readout for optimizing the cross-
entropy loss [11] (see appendix D.2). All simulations were implemented in JAX
[47] and Flax [48] (Apache License 2.0). All datasets were obtained through
the Tensorflow Datasets API [49]. We release the code on
GitHub111https://github.com/Laborieux-Axel/generalized-holo-ep and report the
specific hyperparameters of each simulation in the Appendix D.3.

#### Holomorphic EP matches automatic differentiation at computing
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$.

We sought to understand the amount of bias due to the finite nudge $\beta$
when using different ways of estimating the neuronal error vector
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$. To this end, we trained a two-
hidden layer network with independent forward and backward connections (Fig.
1b) on the Fashion MNIST dataset [50] using stochastic gradient descent for 50
epochs. To investigate the evolution of weight symmetry during training, the
forward and backward connections were initialized symmetrically. We compared
different ways of estimating the neuronal error vector
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$. As the ideal “ceiling” case, we
computed the ground truth using forward-mode automatic differentiation in JAX.
Additionally, we computed the classic one-sided EP estimate [8] denoted by
$N=1$ in Table 1, as well as the Cauchy integral estimate (Eq. 7) computed
with various number $N$ of points (see Appendix D.1). In all cases we ran
simulations with two teaching amplitudes $|\beta|=0.05$ and $|\beta|=0.5$. We
report the final validation errors in Table 1. We observed that the bias of
the nudge is already noticeable at low teaching amplitude when comparing to
better estimates obtained using higher $N$ or the ground truth. The bias
became more noticeable when the teaching amplitude was increased. The lowest
average validation error was obtained with $N=6$ points and $|\beta|=0.5$,
matching the ground truth within statistical uncertainty, and consistent with
our theory. Finally, to test the potential of running hEP in continuous time,
we implemented the continuous-time estimates of Eqs (8) and (9) where the
average spanned five oscillation periods and $N=4$ points were used.
Crucially, there was no free phase in this setting. Despite longer numerical
simulation times in this case (see Appendix D.4), we observed that the
continuous-time estimate only suffered from moderate performance drop with
respect to the ground truth $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ despite
the complete absence of the free phase, thereby confirming our predictions.

Table 1: Validation error in % on Fashion MNIST for different biases level in the estimate of $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$. All values are averages ($n=5$) $\pm$ stddev. Teach. amp. |  $N=1$ [8] | $N=2$ | $N=4$ | $N=6$ | Continuous-time | True $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$
---|---|---|---|---|---|---
$|\beta|=0.05$ | 15.8 $\pm$ 0.6 | 14.6 $\pm$ 0.4 | 14.7 $\pm$ 0.7 | 14.5 $\pm$ 0.6 | 15.8 $\pm$ 0.8 | 14.7 $\pm$ 0.6
$|\beta|=0.5$ | 38.4 $\pm$ 6.2 | 16.3 $\pm$ 0.7 | 14.8 $\pm$ 0.8 | 14.3 $\pm$ 0.6 | 17.3 $\pm$ 1.0
Figure 3: Comparison of training dynamics of a two-hidden-layer MLP on Fashion
MNIST using RBP (top row) and generalized hEP (bottom row). a,b) Evolution of
the angle between forward and backward connections during training on for
varying initial angle. Training tends to reduce weight symmetry. c,d) Learning
curves. The training loss is dependent on initial angle only for generalized
hEP. e,f) Evolution of the cosine similarity between neuronal error vectors of
RBP and hEP over training. Curves are averaged over three seeds and shaded
areas denote $\pm$ 1 stddev.

#### Jacobian asymmetry leads to different learning dynamics.

To see whether and how Jacobian asymmetry influences training, we studied the
dynamics of the neuronal error vectors
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$, obtained through hEP, and
$\boldsymbol{\delta}$ from RBP respectively. We used the same experimental
setting as in the previous section and trained neural networks with varying
degrees of initial weight asymmetry. The degree of initial weight asymmetry
was controlled by setting the backward weights $w_{b}$ as
$w_{b}\leftarrow\sin(\alpha)w_{b}+\cos(\alpha)w_{f}^{\top}$, where $w_{f}$ are
the forward weights. Throughout training, we recorded the evolution of the
angle $\alpha$, and the cosine between both neuronal error vectors (Fig. 3).

We observed that networks trained with RBP all performed comparably in terms
of loss (Fig. 3c) whereas their weight symmetry reduced quicker than in
networks trained with hEP (Fig. 3a,b). Furthermore, RBP trained networks
settled in a parameter regime in which the cosine between
$\boldsymbol{\delta}$ and $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ was
negative for all except the output layer (Fig. 3e,f). In contrast, the
performance of networks trained with hEP strongly depended on the initial
degree of symmetry (Fig. 3d), and although the cosine angle with RBP decreased
during training, it remained positive throughout. These findings could suggest
that learning with hEP alone reaches a parameter regime in which the cosine
with $\boldsymbol{\delta}$ is too low to further decrease the loss while
maintaining the current performance level. In summary, even networks that
start with symmetric weights, lose this symmetry during training, which leads
to increasing bias with respect to the ground-truth neuronal error vector and
impaired task performance. In the next paragraph, we investigate whether and
to what extent the homeostatic loss introduced above (cf Eq. (12)) can
mitigate this problem.

#### Jacobian homeostasis improves functional symmetry and enables training on
larger datasets.

To study the effect of the homeostatic loss (Eq. (12)), we added it scaled by
the hyper parameter $\lambda_{\text{homeo}}$ to the cross entropy loss (see
Appendix D.3). For all experiments we estimated the average in Eq. (12) using
5 stochastic realizations of a Gaussian noise $\boldsymbol{\varepsilon}$ per
sample in the mini-batch. With these settings, we trained the same network as
in the previous paragraphs with the homeostatic loss on Fashion MNIST (Fig.
4a) and found a small reduction in validation error. We further observed that
the added homeostatic loss increased the symmetry of the Jacobian, as defined
by the symmetry measure
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|S\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}/({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|S\right|\kern-1.07639pt\right|\kern-1.07639pt\right|})+{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|})$
(Section 3), over the course of training. In line with this observation, the
alignment between $\boldsymbol{\delta}$ and
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ (Fig. 4c-d) also increased for each
layer.

Because the homeostatic loss is defined on the Jacobian instead of the
weights, we hypothesized that its effect may be more general than merely
increasing weight symmetry. To test this idea, we introduced a network
architecture without feedback connections between adjacent layers, thereby
precluding the possibility for symmetric feedback. Specifically, we trained a
multi-layer architecture in which the output units fed back to the first
hidden layer directly (Fig. 4e) [10, 51, 20]. After training this
antisymmetric architecture, we observed the same key effects on the validation
error, Jacobian symmetry, and error vector alignment (Fig. 4g,h) as in the
reciprocally connected network, thereby confirming our hypothesis that the
homeostatic loss is more general than weight alignment [44].

Figure 4: The homeostatic loss improves hEP training of arbitrary dynamical
systems by acting directly on the Jacobian. Each row corresponds to a specific
architecture-dataset pair. a) Multi-layer architecture where layers are
reciprocally connected. Evolution of the validation error b), Jacobian
symmetry measure c), and layer-wise cosine between neuronal error vectors of
hEP and RBP d) during training of Fashion MNIST. e-h) Same plots as a-d for an
architecture where the output layer directly feeds back to the first layer. i)
Recurrent convolutional architecture used on CIFAR-10. Evolution of validation
error j), angle between weights k) and layer-wise cosine between the neuronal
error vectors of hEP and RBP l). Curves are averaged over five random seeds
for Fashion MNIST and three for CIFAR-10 and shaded areas denote $\pm$ 1
standard error.

Since, the improvement in validation error was measurable, but small in the
network architectures we studied (Fig. 4b,f), we wondered whether the positive
effect of the homeostatic loss would be more pronounced on more challenging
tasks. To this end, we extended a recurrent convolutional architecture (Fig.
4i); 7) to asymmetric feedback weights (see Appendix D.2 for further details).
We trained this architecture using generalized hEP on CIFAR-10, CIFAR-100 [52]
vision benchmarks both with and without homeostatic loss.

Without the homeostatic loss, the network trained on CIFAR-10 reached 60.4%
validation accuracy (Table 2). With the homeostatic loss, the performance
increased to 84.3 (Table 2), with only approximate weight symmetry (Fig. 4k).
Strikingly, this is only a reduction by 4.3% in accuracy in comparison to the
symmetric architecture [7]. When using $N=2$ instead of the ground truth for
estimating $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$, we observed an
additional drop by 2.9% points due to the finite nudge bias. We also noticed
that the homeostatic loss has no measurable effect on performance when
training the asymmetric network with RBP (Table 2), suggesting that the
homeostatic loss is only beneficial while not restricting model performance.
Finally, we reproduced similar findings on CIFAR-100 and ImageNet $32\times
32$, with a remaining 5% gap on Top-5 validation accuracy from the perfectly
symmetric architecture. Together these findings suggest that the addition of
homeostatic objectives is increasingly important for training on larger
datasets and further highlights the need additional homeostatic processes in
the brain [53, 35].

Table 2: Validation accuracy of asymmetric networks trained with generalized hEP (ours) and our implementation of RBP. The symmetric architecture of [7] is reported for comparison. All values are averages ($n=3$) $\pm$ stddev. | CIFAR-10 | CIFAR-100 | ImageNet $32\times 32$
---|---|---|---
| Top-1 (%) | Top-1 (%) | Top-5 (%) | Top-1 (%) | Top-5 (%)
hEP w/o $\mathcal{L}_{\text{homeo}}$ | 60.4 $\pm$ 0.4 | 35.2 $\pm$ 0.3 | 64.4 $\pm$ 0.5 | – | –
hEPN=2,|β|=1 | 81.4 $\pm$ 0.1 | 51.1 $\pm$ 0.8 | 79.2 $\pm$ 0.4 | – | –
hEP (True $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$) | 84.3 $\pm$ 0.1 | 53.8 $\pm$ 0.8 | 81.0 $\pm$ 0.3 | 31.4 $\pm$ 0.1 | 55.2 $\pm$ 0.1
RBP | 87.8 $\pm$ 0.3 | 60.8 $\pm$ 0.2 | 84.4 $\pm$ 0.1 | – | –
RBP w/o $\mathcal{L}_{\text{homeo}}$ | 87.7 $\pm$ 0.2 | 61.0 $\pm$ 0.2 | 85.2 $\pm$ 0.3 | – | –
Sym. hEPN=2,|β|=1 [7] | 88.6 $\pm$ 0.2 | 61.6 $\pm$ 0.1 | 86.0 $\pm$ 0.1 | 36.5 $\pm$ 0.3 | 60.8 $\pm$ 0.4

##  5 Discussion

We have studied a generalized form of EP in the absence of perfect weight
symmetry and for finite nudge amplitudes, a more realistic scenario in
biological neural networks and for constrained neuromorphic systems. We found
that relaxing the strong algorithmic assumptions underlying EP rapidly
degraded its ability to estimate objective function gradients in larger
networks and on real-world benchmark problems. We proceeded to understand the
theoretical origins of this degradation. Since no energy function exists in
the asymmetric setting, we based our analysis on a generalization of EP to
general vector fields with stable fixed points [8]. This analysis identified
two main sources of bias affecting performance that had not been studied in
isolation: The finite nudge amplitude and the Jacobian asymmetry. We then
introduced two effective strategies to counteract these influences. First, we
showed that holomorphic EP [7] generalizes to the asymmetric setting and thus
provides one effective way to diminish bias due to finite nudge amplitude.
Second, we introduce a form of Jacobian homeostasis that, acts directly on the
Jacobian, to ensure functional symmetry. We show that this addition
substantially reduces bias due to weight asymmetry and, more generally,
improves training of networks in which weight symmetry is impossible due to
architectural constraints on the feedback connections (cf. Fig. 4), a setting
particularly interesting for neurobiology. When applied together, both
strategies ultimately allowed us to train deep dynamical networks without
perfect weight symmetry on tasks such as ImageNet $32\times 32$ with only a
small gap in performance to idealized symmetric networks, a crucial
requirement when assessing biological plausibility [54].

The study of asymmetric weights has received a lot of attention both in early
and recent work on artificial neural networks [44, 12, 14, 13, 55, 56]. Our
findings are in keeping with this body of literature as we also found that
learning with fully asymmetric weights is possible on small datasets through
sub-optimal but non-zero alignment with back propagation. However, reminiscent
of studies on Feedback Alignment [12], training becomes more challenging in
deeper networks and on larger datasets where additional objectives are
required [57, 58, 54]. Although our findings do not provide a solution how to
completely dispense with the symmetry requirement for gradient estimate, they
are more general than previous work since the _same_ dynamics and learning
rule is applied to both feed-forward and feedback weights in the network.
Getting rid of the symmetry requirement might require different paradigms to
circumvent the need for feedback altogether, such as [59, 18, 17, 21, 43, 60].

Nevertheless, several questions pertaining to the plausible implementation of
EP remain open. For instance, removing the bias due to finite nudge with hEP
requires the neurons to oscillate in the complex plane. The physical
interpretation of complex-valued neurons is not straightforward, but it could
be implemented through phase coding [61, 62] in addition to firing rates,
which would require the addition of a fast carrier wave, similar to the gamma
rhythm in neurobiology [63]. A detailed study of this potential link to
neuronal oscillations found in the brain is left for future work. Another
limitation of EP is the required convergence to an equilibrium. Since most
current work on EP is performed on digital computers, which simulate the
substrate physics, and not performed on analog substrates on which this
relaxation is achieved through device physics “for free,” [64, 65] the
relaxation phase currently limits the size and depth of models trainable with
EP. This simulation cost is presumably also the reason for EP’s lower task
performance compared to other BP alternatives [17, 18, 55, 56]. It is not
clear whether slow convergence would be an issue in neurobiology, where
theoretical work on spiking neural networks suggests that they may be ideally
suited for rapid relaxation dynamics toward a quasi-equilibrium through
efficient coding [66, 67, 68, 69, 70, 71, 72]. Moreover, the potential issue
of fast convergence is shared broadly among the class of EBMs [73], and has
sparked research towards fundamentally new ways of performing credit
assignment [20, 15, 19]. Finally, EBM architectures such as modern Hopfield
networks [74, 75], or different ways of reaching equilibrium [76, 67, 68, 77]
could help make the algorithm faster and more plausible.

In summary, our work strengthens EP as a theoretical framework for
understanding learning processes in physical neural networks, a key
requirement for designing power efficient systems that learn with their own
relaxation dynamics.

## Acknowledgments and Disclosure of Funding

We thank all members of the Zenke Group for comments and discussions. We thank
Maxence Ernoult, Nicolas Zucchet, Jack Kendall and Benjamin Scellier for
helpful feedback. This project was supported by the Swiss National Science
Foundation [grant numbers PCEFP3_202981 and TMPFP3_210282] and the Novartis
Research Foundation. The authors declare no competing interests.

## References

  * Rumelhart et al. [1986] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning representations by back-propagating errors.  _Nature_ , 323(6088):533–536, 1986.
  * Lillicrap et al. [2020] Timothy P Lillicrap, Adam Santoro, Luke Marris, Colin J Akerman, and Geoffrey Hinton.  Backpropagation and the brain.  _Nature Reviews Neuroscience_ , 21(6):335–346, 2020.
  * Whittington and Bogacz [2019] James CR Whittington and Rafal Bogacz.  Theories of error back-propagation in the brain.  _Trends in cognitive sciences_ , 23(3):235–250, 2019.
  * Grossberg [1987] Stephen Grossberg.  Competitive learning: From interactive activation to adaptive resonance.  _Cognitive science_ , 11(1):23–63, 1987.
  * Scellier and Bengio [2017] Benjamin Scellier and Yoshua Bengio.  Equilibrium propagation: Bridging the gap between energy-based models and backpropagation.  _Frontiers in computational neuroscience_ , 11:24, 2017\.
  * Hopfield [1984] John J Hopfield.  Neurons with graded response have collective computational properties like those of two-state neurons.  _Proceedings of the national academy of sciences_ , 81(10):3088–3092, 1984.
  * Laborieux and Zenke [2022] Axel Laborieux and Friedemann Zenke.  Holomorphic equilibrium propagation computes exact gradients through finite size oscillations.  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, _Advances in Neural Information Processing Systems_ , volume 35, pages 12950–12963. Curran Associates, Inc., 2022.  URL https://proceedings.neurips.cc/paper_files/paper/2022/file/545a114e655f9d25ba0d56ea9a01fc6e-Paper-Conference.pdf.
  * Scellier et al. [2018] Benjamin Scellier, Anirudh Goyal, Jonathan Binas, Thomas Mesnard, and Yoshua Bengio.  Generalization of equilibrium propagation to vector field dynamics.  _arXiv preprint arXiv:1808.04873_ , 2018.
  * Ernoult et al. [2020] Maxence Ernoult, Julie Grollier, Damien Querlioz, Yoshua Bengio, and Benjamin Scellier.  Equilibrium propagation with continual weight updates.  _arXiv preprint arXiv:2005.04168_ , 2020.
  * Kohan et al. [2018] Adam A. Kohan, Edward A. Rietman, and Hava T. Siegelmann.  Error Forward-Propagation: Reusing Feedforward Connections to Propagate Errors in Deep Learning, August 2018.
  * Laborieux et al. [2021] Axel Laborieux, Maxence Ernoult, Benjamin Scellier, Yoshua Bengio, Julie Grollier, and Damien Querlioz.  Scaling equilibrium propagation to deep convnets by drastically reducing its gradient estimator bias.  _Frontiers in neuroscience_ , 15:129, 2021.
  * Lillicrap et al. [2016] Timothy P Lillicrap, Daniel Cownden, Douglas B Tweed, and Colin J Akerman.  Random synaptic feedback weights support error backpropagation for deep learning.  _Nature communications_ , 7(1):1–10, 2016.
  * Nøkland [2016] Arild Nøkland.  Direct feedback alignment provides learning in deep neural networks.  _Advances in neural information processing systems_ , 29, 2016.
  * Launay et al. [2020] Julien Launay, Iacopo Poli, François Boniface, and Florent Krzakala.  Direct feedback alignment scales to modern deep learning tasks and architectures.  _Advances in neural information processing systems_ , 33:9346–9360, 2020.
  * Clark et al. [2021] David G. Clark, L. F. Abbott, and SueYeon Chung.  Credit Assignment Through Broadcasting a Global Error Vector, October 2021.
  * Refinetti et al. [2021] Maria Refinetti, Stéphane d’Ascoli, Ruben Ohana, and Sebastian Goldt.  Align, then memorise: The dynamics of learning with feedback alignment, June 2021.
  * Ren et al. [2023] Mengye Ren, Simon Kornblith, Renjie Liao, and Geoffrey Hinton.  Scaling forward gradient with local losses.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=JxpBP1JM15-.
  * Journé et al. [2023] Adrien Journé, Hector Garcia Rodriguez, Qinghai Guo, and Timoleon Moraitis.  Hebbian deep learning without feedback.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=8gd4M-_Rj1.
  * Hinton [2022] Geoffrey Hinton.  The forward-forward algorithm: Some preliminary investigations.  _arXiv preprint arXiv:2212.13345_ , 2022.
  * Dellaferrera and Kreiman [2022] Giorgia Dellaferrera and Gabriel Kreiman.  Error-driven Input Modulation: Solving the Credit Assignment Problem without a Backward Pass.  In _Proceedings of the 39th International Conference on Machine Learning_ , pages 4937–4955. PMLR, June 2022.
  * Löwe et al. [2019] Sindy Löwe, Peter O’Connor, and Bastiaan Veeling.  Putting an end to end-to-end: Gradient-isolated learning of representations.  _Advances in neural information processing systems_ , 32, 2019.
  * Baldi [1995] Pierre Baldi.  Gradient descent learning algorithm overview: A general dynamical systems perspective.  _IEEE Transactions on neural networks_ , 6(1):182–195, 1995.
  * Scarselli et al. [2008] Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini.  The graph neural network model.  _IEEE transactions on neural networks_ , 20(1):61–80, 2008.
  * Liao et al. [2018] Renjie Liao, Yuwen Xiong, Ethan Fetaya, Lisa Zhang, KiJung Yoon, Xaq Pitkow, Raquel Urtasun, and Richard Zemel.  Reviving and improving recurrent back-propagation.  In _International Conference on Machine Learning_ , pages 3082–3091. PMLR, 2018.
  * Almeida [1990] Luis B Almeida.  A learning rule for asynchronous perceptrons with feedback in a combinatorial environment.  In _Artificial neural networks: concept learning_ , pages 102–111. IEEE Computer Society Press, 1990.
  * Pineda [1987] Fernando Pineda.  Generalization of back propagation to recurrent and higher order neural networks.  In _Neural information processing systems_ , 1987.
  * Bai et al. [2019] Shaojie Bai, J Zico Kolter, and Vladlen Koltun.  Deep equilibrium models.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Scellier and Bengio [2019] Benjamin Scellier and Yoshua Bengio.  Equivalence of equilibrium propagation and recurrent backpropagation.  _Neural computation_ , 31(2):312–329, 2019.
  * Appel [2007] Walter Appel.  _Mathematics for physics and physicists_.  Princeton University Press, 2007.
  * Williams et al. [2023] Ezekiel Williams, Colin Bredenberg, and Guillaume Lajoie.  Flexible phase dynamics for bio-plausible contrastive learning.  _arXiv preprint arXiv:2302.12431_ , 2023.
  * Dan and Poo [2004] Yang Dan and Mu-ming Poo.  Spike timing-dependent plasticity of neural circuits.  _Neuron_ , 44(1):23–30, 2004.
  * Clopath and Gerstner [2010] Claudia Clopath and Wulfram Gerstner.  Voltage and spike timing interact in stdp–a unified model.  _Frontiers in synaptic neuroscience_ , 2:25, 2010.
  * Bai et al. [2021] Shaojie Bai, Vladlen Koltun, and J Zico Kolter.  Stabilizing equilibrium models by jacobian regularization.  _arXiv preprint arXiv:2106.14342_ , 2021.
  * Hutchinson [1989] Michael F Hutchinson.  A stochastic estimator of the trace of the influence matrix for laplacian smoothing splines.  _Communications in Statistics-Simulation and Computation_ , 18(3):1059–1076, 1989.
  * Stock et al. [2022] Christopher H Stock, Sarah E Harvey, Samuel A Ocko, and Surya Ganguli.  Synaptic balancing: A biologically plausible local learning rule that provably increases neural network noise robustness without sacrificing task performance.  _PLOS Computational Biology_ , 18(9):e1010418, 2022.
  * Meulemans et al. [2020] Alexander Meulemans, Francesco Carzaniga, Johan Suykens, João Sacramento, and Benjamin F Grewe.  A theoretical framework for target propagation.  _Advances in Neural Information Processing Systems_ , 33:20024–20036, 2020.
  * Ernoult et al. [2022] Maxence Ernoult, Fabrice Normandin, Abhinav Moudgil, Sean Spinney, Eugene Belilovsky, Irina Rish, Blake Richards, and Yoshua Bengio.  Towards scaling difference target propagation by learning backprop targets.  _arXiv preprint arXiv:2201.13415_ , 2022.
  * Lee et al. [2015] Dong-Hyun Lee, Saizheng Zhang, Asja Fischer, and Yoshua Bengio.  Difference target propagation.  In _Joint european conference on machine learning and knowledge discovery in databases_ , pages 498–515. Springer, 2015.
  * Bengio [2014] Yoshua Bengio.  How auto-encoders could provide credit assignment in deep networks via target propagation.  _arXiv preprint arXiv:1407.7906_ , 2014.
  * Le Cun et al. [1988] Yann Le Cun, Conrad Galland, and Geoffrey E Hinton.  Gemini: gradient estimation through matrix inversion after noise injection.  _Advances in neural information processing systems_ , 1, 1988.
  * Wiskott and Sejnowski [2002] Laurenz Wiskott and Terrence J Sejnowski.  Slow feature analysis: Unsupervised learning of invariances.  _Neural computation_ , 14(4):715–770, 2002.
  * Lipshutz et al. [2020] David Lipshutz, Charles Windolf, Siavash Golkar, and Dmitri Chklovskii.  A Biologically Plausible Neural Network for Slow Feature Analysis.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 14986–14996. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper/2020/hash/ab73f542b6d60c4de151800b8abc0a6c-Abstract.html.
  * Halvagal and Zenke [2022] Manu Srinath Halvagal and Friedemann Zenke.  The combination of hebbian and predictive plasticity learns invariant object representations in deep sensory networks.  _bioRxiv_ , 2022.
  * Kolen and Pollack [1994] John F Kolen and Jordan 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.
  * Akrout et al. [2019] Mohamed Akrout, Collin Wilson, Peter Humphreys, Timothy Lillicrap, and Douglas B Tweed.  Deep Learning without Weight Transport.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Ernoult et al. [2019] Maxence Ernoult, Julie Grollier, Damien Querlioz, Yoshua Bengio, and Benjamin Scellier.  Updates of equilibrium prop match gradients of backprop through time in an rnn with static input.  _Advances in neural information processing systems_ , 32, 2019.
  * Bradbury et al. [2018] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang.  JAX: composable transformations of Python+NumPy programs, 2018\.  URL http://github.com/google/jax.
  * Heek et al. [2020] Jonathan Heek, Anselm Levskaya, Avital Oliver, Marvin Ritter, Bertrand Rondepierre, Andreas Steiner, and Marc van Zee.  Flax: A neural network library and ecosystem for JAX, 2020.  URL http://github.com/google/flax.
  * Abadi et al. [2015] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dan Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  URL http://tensorflow.org/.  Software available from tensorflow.org.
  * Xiao et al. [2017] Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.  _arXiv preprint arXiv:1708.07747_ , 2017.
  * Kohan et al. [2023] Adam Kohan, Edward A. Rietman, and Hava T. Siegelmann.  Signal Propagation: The Framework for Learning and Inference in a Forward Pass.  _IEEE Transactions on Neural Networks and Learning Systems_ , pages 1–12, 2023.  ISSN 2162-237X, 2162-2388.  doi: 10.1109/TNNLS.2022.3230914.
  * Krizhevsky [2009] Alex Krizhevsky.  Learning multiple layers of features from tiny images.  Technical report, University of Toronto, 2009.
  * Zenke et al. [2017] Friedemann Zenke, Wulfram Gerstner, and Surya Ganguli.  The temporal paradox of Hebbian learning and homeostatic plasticity.  _Curr Opin Neurobiol_ , 43:166–176, April 2017.  ISSN 0959-4388.  doi: 10.1016/j.conb.2017.03.015.  URL http://www.sciencedirect.com/science/article/pii/S0959438817300910.
  * Bartunov et al. [2018] Sergey Bartunov, Adam Santoro, Blake Richards, Luke Marris, Geoffrey E Hinton, and Timothy Lillicrap.  Assessing the scalability of biologically-motivated deep learning algorithms and architectures.  _Advances in neural information processing systems_ , 31, 2018.
  * Payeur et al. [2021] Alexandre Payeur, Jordan Guerguiev, Friedemann Zenke, Blake A Richards, and Richard Naud.  Burst-dependent synaptic plasticity can coordinate learning in hierarchical circuits.  _Nature neuroscience_ , 24(7):1010–1019, 2021\.
  * Greedy et al. [2022] Will Greedy, Heng Wei Zhu, Joseph Pemberton, Jack Mellor, and Rui Ponte Costa.  Single-phase deep learning in cortico-cortical networks.  _Advances in Neural Information Processing Systems_ , 35:24213–24225, 2022.
  * Xiao et al. [2018] Will Xiao, Honglin Chen, Qianli Liao, and Tomaso Poggio.  Biologically-plausible learning algorithms can scale to large datasets, December 2018.
  * Liao et al. [2016] Qianli Liao, Joel Z. Leibo, and Tomaso Poggio.  How Important is Weight Symmetry in Backpropagation?, February 2016.
  * Frenkel et al. [2021] Charlotte Frenkel, Martin Lefebvre, and David Bol.  Learning Without Feedback: Fixed Random Learning Signals Allow for Feedforward Training of Deep Neural Networks.  _Frontiers in Neuroscience_ , 15:629892, February 2021.  ISSN 1662-453X.  doi: 10.3389/fnins.2021.629892.
  * Moraitis et al. [2022] Timoleon Moraitis, Dmitry Toichkin, Adrien Journé, Yansong Chua, and Qinghai Guo.  Softhebb: Bayesian inference in unsupervised hebbian soft winner-take-all networks.  _Neuromorphic Computing and Engineering_ , 2(4):044017, 2022.
  * Frady and Sommer [2019] E Paxon Frady and Friedrich T Sommer.  Robust computation with rhythmic spike patterns.  _Proceedings of the National Academy of Sciences_ , 116(36):18050–18059, 2019.
  * Bybee et al. [2022] Connor Bybee, E Paxon Frady, and Friedrich T Sommer.  Deep learning in spiking phasor neural networks.  _arXiv preprint arXiv:2204.00507_ , 2022.
  * Engel et al. [2001] Andreas K Engel, Pascal Fries, and Wolf Singer.  Dynamic predictions: oscillations and synchrony in top–down processing.  _Nature Reviews Neuroscience_ , 2(10):704–716, 2001.
  * Yi et al. [2023] Su-in Yi, Jack D Kendall, R Stanley Williams, and Suhas Kumar.  Activity-difference training of deep neural networks using memristor crossbars.  _Nature Electronics_ , 6(1):45–51, 2023.
  * Kendall et al. [2020] Jack Kendall, Ross Pantone, Kalpana Manickavasagam, Yoshua Bengio, and Benjamin Scellier.  Training end-to-end analog neural networks with equilibrium propagation.  _arXiv preprint arXiv:2006.01981_ , 2020.
  * Denève and Machens [2016] Sophie Denève and Christian K. Machens.  Efficient codes and balanced networks.  _Nat Neurosci_ , 19(3):375–382, March 2016.  ISSN 1097-6256.  doi: 10.1038/nn.4243.  URL http://www.nature.com/neuro/journal/v19/n3/full/nn.4243.html.
  * Haider et al. [2021] Paul Haider, Benjamin Ellenberger, Laura Kriener, Jakob Jordan, Walter Senn, and Mihai A Petrovici.  Latent equilibrium: Arbitrarily fast computation with arbitrarily slow neurons.  _Advances in Neural Information Processing Systems_ , 34:17839–17851, 2021.
  * Senn et al. [2023] Walter Senn, Dominik Dold, Akos F Kungl, Benjamin Ellenberger, Jakob Jordan, Yoshua Bengio, Joao Sacramento, and Mihai A Petrovici.  A neuronal least-action principle for real-time learning in cortical circuits.  _bioRxiv_ , pages 2023–03, 2023.
  * Sacramento et al. [2018] João Sacramento, Rui Ponte Costa, Yoshua Bengio, and Walter Senn.  Dendritic cortical microcircuits approximate the backpropagation algorithm.  _Advances in neural information processing systems_ , 31, 2018.
  * Meulemans et al. [2021] Alexander Meulemans, Matilde Tristany Farinha, Javier García Ordóñez, Pau Vilimelis Aceituno, João Sacramento, and Benjamin F Grewe.  Credit assignment in neural networks through deep feedback control.  _Advances in Neural Information Processing Systems_ , 34:4674–4687, 2021.
  * Meulemans et al. [2022a] Alexander Meulemans, Nicolas Zucchet, Seijin Kobayashi, Johannes Von Oswald, and João Sacramento.  The least-control principle for local learning at equilibrium.  _Advances in Neural Information Processing Systems_ , 35:33603–33617, 2022a.
  * Meulemans et al. [2022b] Alexander Meulemans, Matilde Tristany Farinha, Maria R Cervera, João Sacramento, and Benjamin F Grewe.  Minimizing control for credit assignment with strong feedback.  In _International Conference on Machine Learning_ , pages 15458–15483. PMLR, 2022b.
  * Zucchet and Sacramento [2022] Nicolas Zucchet and João Sacramento.  Beyond backpropagation: bilevel optimization through implicit differentiation and equilibrium propagation.  _Neural Computation_ , 34(12):2309–2346, 2022\.
  * Krotov and Hopfield [2016] Dmitry Krotov and John J Hopfield.  Dense associative memory for pattern recognition.  _Advances in neural information processing systems_ , 29, 2016.
  * Ramsauer et al. [2020] Hubert Ramsauer, Bernhard Schäfl, Johannes Lehner, Philipp Seidl, Michael Widrich, Thomas Adler, Lukas Gruber, Markus Holzleitner, Milena Pavlović, Geir Kjetil Sandve, et al.  Hopfield networks is all you need.  _arXiv preprint arXiv:2008.02217_ , 2020.
  * Luczak et al. [2022] Artur Luczak, Bruce L McNaughton, and Yoshimasa Kubo.  Neurons learn by predicting future activity.  _Nature Machine Intelligence_ , pages 1–11, 2022.
  * Zucchet et al. [2021] Nicolas Zucchet, Simon Schug, Johannes von Oswald, Dominic Zhao, and João Sacramento.  A contrastive rule for meta-learning.  _arXiv preprint arXiv:2104.01677_ , 2021.
  * Stergiou et al. [2021] A Stergiou, R Poppe, and G Kalliatakis.  Refining activation downsampling with softpool. arxiv 2021.  _arXiv preprint arXiv:2101.00440_ , 2021.
  * Elfwing et al. [2018] Stefan Elfwing, Eiji Uchibe, and Kenji Doya.  Sigmoid-weighted linear units for neural network function approximation in reinforcement learning.  _Neural Networks_ , 107:3–11, 2018.
  * Loshchilov and Hutter [2017] Ilya Loshchilov and Frank Hutter.  Decoupled weight decay regularization.  _arXiv preprint arXiv:1711.05101_ , 2017.

##  Appendix A Mathematical background

###  A.1 Complex analysis

Here we provide intuitions and the minimal theoretical background on complex
analysis to allow appreciating the results presented in Sections 3.1 and 3.2.
We refer the interested reader to [29] for a more complete introduction.

The notion of complex differentiability in a point $z_{0}\in\mathbb{C}$ is
defined for a function $f:z\in\mathbb{C}\mapsto f(z)\in\mathbb{C}$ by the
existence of the limit:

$\lim_{z\to z_{0}}\frac{f(z)-f(z_{0})}{z-z_{0}}.$ (13)

The limit is written $f^{\prime}(z_{0})$, and the function $f$ is said to be
“holomorphic” in $z_{0}$. Although the definition is similar to
$\mathbb{R}$-differentiability, the fact that the limit is taken with a
complex number that can go to $z_{0}$ from all directions in the complex plane
makes the variation of $f$ more “rigid”. However, most usual functions on real
numbers can be extended to complex inputs and outputs and are holomorphic,
such as the exponential, the sigmoid, all polynomials, trigonometric
functions, softmax, and logarithm. Non-holomorphic functions include functions
that, for instance, use the absolute value or min/max operators.

In exchange for the rigidness of holomorphic functions, one obtains the
possibility to express their derivatives at any order through integrals as
described by the Cauchy formulas.

$f^{(n)}(z_{0})=\frac{n!}{2\mathrm{i}\pi}\oint_{\gamma}\frac{f(z)}{(z-z_{0})^{n+1}}\mathrm{d}z,$
(14)

where $n!$ is the factorial product of integers up to $n$, $\mathrm{i}$ is the
imaginary unit, and the integral is taken over a path
$\gamma\subset\mathbb{C}$ going around $z_{0}$ once and counterclockwise, and
provided that $f$ is holomorphic on a set that includes the path. For
instance, for $n=0$ and the path $\gamma:\theta\in[-\pi,\pi]\mapsto
z_{0}+re^{\mathrm{i}\theta}$, doing the change of variable
$\mathrm{d}z=r\mathrm{i}e^{\mathrm{i}\theta}\mathrm{d}\theta$ yields the mean
value property linking the value of $f$ in $z_{0}$ to its variation on a
circle of radius $r$ going around $z_{0}$:

$\displaystyle f(z_{0})$
$\displaystyle=\frac{1}{2\mathrm{i}\pi}\int_{-\pi}^{\pi}\frac{f(z_{0}+re^{i\theta})}{z_{0}+re^{\mathrm{i}\theta}-z_{0}}r\mathrm{i}e^{\mathrm{i}\theta}\mathrm{d}\theta,$
$\displaystyle f(z_{0})$
$\displaystyle=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(z_{0}+re^{i\theta})\mathrm{d}\theta.$
(15)

And $f^{\prime}(z_{0})$ can be computed by plugging $n=1$ in Eq. (14), and the
same change of variable:

$\displaystyle f^{\prime}(z_{0})$
$\displaystyle=\frac{1}{2\mathrm{i}\pi}\oint_{\gamma}\frac{f(z)}{(z-z_{0})^{2}}\mathrm{d}z,$
$\displaystyle f^{\prime}(z_{0})$
$\displaystyle=\frac{1}{2\mathrm{i}\pi}\int_{-\pi}^{\pi}\frac{f(z_{0}+re^{i\theta})}{(z_{0}+re^{\mathrm{i}\theta}-z_{0})^{2}}r\mathrm{i}e^{\mathrm{i}\theta}\mathrm{d}\theta,$
$\displaystyle f^{\prime}(z_{0})$ $\displaystyle=\frac{1}{2\pi
r}\int_{-\pi}^{\pi}f(z_{0}+re^{i\theta})e^{-\mathrm{i}\theta}\mathrm{d}\theta.$
(16)

Importantly, these formulas are integrals over a non-vanishing radius $r>0$
where the integrands only involve the function $f$, which is a promising
conceptual step toward biological plausiblity and hardware design since
integrals are easier to implement in noisy substrates than finite differences.

###  A.2 Hutchinson trace estimator

The Hutchinson trace estimator [34] refers to the following “trick” to
estimate the trace of a square matrix $M=[m_{ij}]$ : suppose we have a random
vector $\boldsymbol{\varepsilon}\in\mathbb{R}^{n}$ such that
$\mathbb{E}[\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^{\top}]=I$, where
$I$ is the identity matrix. Then we have:

$\displaystyle\mathbb{E}[\boldsymbol{\varepsilon}^{\top}M\boldsymbol{\varepsilon}]$
$\displaystyle=\mathbb{E}\left[\sum_{ij}m_{ij}\varepsilon_{i}\varepsilon_{j}\right]$
$\displaystyle=\sum_{ij}m_{ij}\mathbb{E}\left[\varepsilon_{i}\varepsilon_{j}\right]$
$\displaystyle=\sum_{ij}m_{ij}I_{ij}$
$\displaystyle=\sum_{i}m_{ii}=\operatorname{Tr}(M).$ (17)

Since the Frobenius norm can be expressed with the trace operator, it can be
estimated in the same way. For instance, the Frobenius norm of the asymmetric
part $A$ of the Jacobian $J_{F}$, $A=\tfrac{1}{2}(J_{F}-J_{F}^{\top})$ is
given by:

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}$
$\displaystyle=\operatorname{Tr}(A^{\top}A)$
$\displaystyle=\frac{1}{4}\operatorname{Tr}\left((J_{F}^{\top}-J_{F})(J_{F}-J_{F}^{\top})\right)$
$\displaystyle=\frac{1}{4}\operatorname{Tr}\left(J_{F}^{\top}J_{F}-J_{F}^{\top}J_{F}^{\top}-J_{F}J_{F}+J_{F}J_{F}^{\top}\right)$
$\displaystyle=\frac{1}{2}\operatorname{Tr}(J_{F}^{\top}J_{F})-\frac{1}{2}\operatorname{Tr}(J_{F}^{2})$
$\displaystyle=\frac{1}{2}\mathbb{E}[\boldsymbol{\varepsilon}^{\top}J_{F}^{\top}J_{F}\boldsymbol{\varepsilon}]-\frac{1}{2}\mathbb{E}[\boldsymbol{\varepsilon}^{\top}J_{F}^{2}\boldsymbol{\varepsilon}]\quad\text{using
Eq.~{}\eqref{eq:hutchinson}}$
$\displaystyle=\frac{1}{2}\mathbb{E}\left[\|J_{F}\boldsymbol{\varepsilon}\|^{2}-\boldsymbol{\varepsilon}^{\top}J_{F}^{2}\boldsymbol{\varepsilon}\right].$
(18)

To obtain the fourth line, we used the linear property of the trace operator,
the fact that the trace is unchanged by the matrix transposition, and fact
that $\operatorname{Tr}(M_{1}M_{2})=\operatorname{Tr}(M_{2}M_{1})$ for any two
matrices $M_{1}$ and $M_{2}$, which is a consequence of the symmetry of the
Frobenius inner product. The final quantity in Eq. (A.2) can be estimated
efficiently with vector-Jacobian product routines available in automatic
differentiation frameworks such as JAX or PyTorch, without instantiating the
full Jacobian matrix.

##  Appendix B Relation between the energy and vector field gradient formulas

Here we briefly provide an intuition for how the EP gradient formula obtained
with the energy function formalism [5] relates to the EP gradient formula
obtained directly with the vector field (Eq. (6), [8]). Although the formulas
seem different, they are in fact equivalent. The main takeaway is that the
connections $w_{ij}$ and $w_{ji}$ are fused into a single parameter in the
energy function formalism, whereas the vector field formalism distinguishes
the two. For a continuous Hopfield model with two-body interaction terms of
the form $-\sigma(u_{i})w_{ij}\sigma(u_{j})$, the energy-based (EB) gradient
of the loss with respect to $w_{ij}$ is [5]:

$\displaystyle\left.\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}w_{ij}}\right|_{\text{EB}}=\left.\frac{\mathrm{d}\sigma(u_{i}^{\ast})\sigma(u_{j}^{\ast})}{\mathrm{d}\beta}\right|_{\beta=0}.$

If we use the product rule of derivatives, we obtain:

$\left.\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}w_{ij}}\right|_{\text{EB}}=\sigma(u_{i}^{\ast})\sigma^{\prime}(u_{j}^{\ast})\left.\frac{\mathrm{d}u_{j}^{\ast}}{\mathrm{d}\beta}\right|_{\beta=0}+\sigma(u_{j}^{\ast})\sigma^{\prime}(u_{i}^{\ast})\left.\frac{\mathrm{d}u_{i}^{\ast}}{\mathrm{d}\beta}\right|_{\beta=0}.$

Both terms in the sum corresponds respectively to the vector field (VF)
gradient formula for $w_{ij}$ and $w_{ji}$. Here $\frac{\partial F}{\partial
w_{ij}}=\sigma(u_{i}^{\ast})\sigma^{\prime}(u_{j}^{\ast})$. Therefore, we
have:

$\left.\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}w_{ij}}\right|_{\text{EB}}=\left.\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}w_{ij}}\right|_{\text{VF}}+\left.\frac{\mathrm{d}\mathcal{L}}{\mathrm{d}w_{ji}}\right|_{\text{VF}}.$

##  Appendix C Theoretical proofs

###  C.1 Convergence of the continuous-time estimate

We recall the limit:

$\frac{1}{t|\beta|}\int_{0}^{t}\boldsymbol{u}^{\ast}_{\beta(\tau)}e^{-2\mathrm{i}\pi\tau/T}\mathrm{d}\tau\underset{t\to\infty}{\longrightarrow}\left.\frac{\mathrm{d}\boldsymbol{u}^{\ast}}{\mathrm{d}\beta}\right|_{\beta=0}.$

We split the running time $t$ into an integer $k$ amount of periods
$t=kT+t^{\prime}$. Then by splitting the integral we have:

$\displaystyle\frac{1}{t|\beta|}\int_{0}^{t}\boldsymbol{u}^{\ast}_{\beta(\tau)}e^{-2\mathrm{i}\pi\tau/T}\mathrm{d}\tau$
$\displaystyle=\frac{1}{t|\beta|}\int_{0}^{kT+t^{\prime}}\boldsymbol{u}^{\ast}_{\beta(\tau)}e^{-2\mathrm{i}\pi\tau/T}\mathrm{d}\tau$
$\displaystyle=\frac{1}{t|\beta|}\left(k\int_{0}^{T}\boldsymbol{u}^{\ast}_{\beta(\tau)}e^{-2\mathrm{i}\pi\tau/T}\mathrm{d}\tau+\int_{0}^{t^{\prime}}\boldsymbol{u}^{\ast}_{\beta(\tau)}e^{-2\mathrm{i}\pi\tau/T}\mathrm{d}\tau\right)$
$\displaystyle=\frac{1}{t|\beta|}\left(kT|\beta|\left.\frac{\mathrm{d}\boldsymbol{u}^{\ast}}{\mathrm{d}\beta}\right|_{\beta=0}+\int_{0}^{t^{\prime}}\boldsymbol{u}^{\ast}_{\beta(\tau)}e^{-2\mathrm{i}\pi\tau/T}\mathrm{d}\tau\right)$
$\displaystyle=\frac{t-t^{\prime}}{t}\left.\frac{\mathrm{d}\boldsymbol{u}^{\ast}}{\mathrm{d}\beta}\right|_{\beta=0}+\frac{1}{t|\beta|}\int_{0}^{t^{\prime}}\boldsymbol{u}^{\ast}_{\beta(\tau)}e^{-2\mathrm{i}\pi\tau/T}\mathrm{d}\tau.$
(19)

Since $0<t^{\prime}<T$, the first term converges toward
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$, and the second term is a bounded
integral divided by $t$, and therefore goes to $0$ as $t\to\infty$.

For the pre-synaptic term we recall the statement:

$\frac{1}{t}\int_{0}^{t}\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{\beta(\tau)})\mathrm{d}\tau\underset{t\to\infty}{\longrightarrow}\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0}).$

The proof is the same as above, once we have shown that:

$\frac{1}{T}\int_{0}^{T}\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{\beta(\tau)})\mathrm{d}\tau=\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0}).$ (20)

By applying the Mean Value Theorem (Eq. (A.1)) [29] to the function
$\varphi:\beta\in\mathcal{D}_{\boldsymbol{u}^{\ast}}\mapsto\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{\beta})$ in $\beta=0$,
we have that:

$\frac{\partial
F}{\partial\boldsymbol{\theta}}(\boldsymbol{u}^{\ast}_{0})=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{\partial
F}{\partial\boldsymbol{\theta}}\left(\boldsymbol{u}^{\ast}_{re^{\mathrm{i}\alpha}}\right)\mathrm{d}\alpha,$

where $r$ is a radius such that the circle is contained in the domain of
definition $\mathcal{D}_{\boldsymbol{u}^{\ast}}$ of $\varphi$. Then we apply
the equation in $r=|\beta|$ and do the change of variable $\alpha=2\pi\tau/T$,
which gives Eq. (20).

###  C.2 Bias term from Jacobian asymmetry

We recall that $\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ and
$\boldsymbol{\delta}$ are linked by:

$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}=J_{F}(\boldsymbol{u}^{\ast}_{0})^{-1}J_{F}(\boldsymbol{u}^{\ast}_{0})^{\top}\boldsymbol{\delta}.$
(21)

Let $S$ and $A$ be the symmetric and skew-symmetric parts of $J_{F}$, such
that
$S=(J_{F}(\boldsymbol{u}^{\ast}_{0})+J_{F}(\boldsymbol{u}^{\ast}_{0})^{\top})/2$
and
$A=(J_{F}(\boldsymbol{u}^{\ast}_{0})-J_{F}(\boldsymbol{u}^{\ast}_{0})^{\top})/2$.
Since $J_{F}(\boldsymbol{u}^{\ast}_{0})$ is assumed invertible and
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
is assumed small, we have that $S$ is invertible and:

$\displaystyle
J_{F}(\boldsymbol{u}^{\ast}_{0})^{-1}J_{F}(\boldsymbol{u}^{\ast}_{0})^{\top}$
$\displaystyle=(S+A)^{-1}(S-A)$ $\displaystyle=(I+S^{-1}A)^{-1}S^{-1}(S-A)$
$\displaystyle=(I+S^{-1}A)^{-1}(I-S^{-1}A)$
$\displaystyle=\left(\sum_{k=0}^{\infty}(-1)^{k}(S^{-1}A)^{k}\right)(I-S^{-1}A)$
$\displaystyle=\sum_{k=0}^{\infty}(-1)^{k}(S^{-1}A)^{k}-\sum_{k=0}^{\infty}(-1)^{k}(S^{-1}A)^{k+1}$
$\displaystyle=\sum_{k=0}^{\infty}(-1)^{k}(S^{-1}A)^{k}+\sum_{k=1}^{\infty}(-1)^{k}(S^{-1}A)^{k}$
$\displaystyle=I-2S^{-1}A+2\sum_{k=2}^{\infty}(-1)^{k}(S^{-1}A)^{k}.$

We used the Neumann Series to compute the inverse $(I+S^{-1}A)^{-1}$, which
converges given the assumption
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|S^{-1}A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\to
0$. From Eq (21) we can write:

$\displaystyle\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}=\boldsymbol{\delta}-2S^{-1}A\boldsymbol{\delta}+o(S^{-1}A\boldsymbol{\delta})$

##  Appendix D Simulation details

###  D.1 Numerical estimate of the Cauchy integral

For evaluating the integral in Eq. (7), we set an integer number of $N$ values
of $\beta$:

$\displaystyle\beta_{k}=|\beta|e^{2\mathrm{i}\pi k/N},\quad
k\in[0,1,...,N-1].$

The integral is then approximated by:

$\frac{1}{T|\beta|}\int_{0}^{T}\boldsymbol{u}^{\ast}_{\beta(\tau)}e^{-2\mathrm{i}\pi\tau/T}\mathrm{d}\tau\approx\frac{1}{N|\beta|}\sum_{k=0}^{N}\boldsymbol{u}^{\ast}_{\beta_{k}}e^{-2\mathrm{i}\pi
k/N}.$ (22)

The leading term of the bias when using $N$ points is $|\beta|^{N}/(N+1)!$,
where $!$ is the factorial product. However, unless specified we used forward
mode automatic differentiation to compute the ground truth
$\mathrm{d}_{\beta}\boldsymbol{u}^{\ast}$ without bias.

###  D.2 Architecture and dynamics

#### Multi-layer networks.

The disrete dynamics of the multi-layer networks of experiments in Figs. 2, 3,
4 on Fashion MNIST read:

$\displaystyle\left\\{\begin{array}[]{lll}\boldsymbol{u}_{1}&\leftarrow&\mathbf{w}_{in}\boldsymbol{x}+\sum_{l^{\prime}\text{
to
}1}\mathbf{w}_{1,l^{\prime}}\sigma(\boldsymbol{u}_{l^{\prime}})+\mathbf{b}_{1},\\\\[8.5359pt]
\boldsymbol{u}_{l}&\leftarrow&\sum_{l^{\prime}\text{ to
}l}\mathbf{w}_{l,l^{\prime}}\sigma(\boldsymbol{u}_{l^{\prime}})+\mathbf{b}_{l},\\\\[8.5359pt]
\boldsymbol{u}_{L}&\leftarrow&\mathbf{w}_{L,L-1}^{f}\sigma(\boldsymbol{u}_{L-1})+\mathbf{b}_{L}+\beta\mathbf{w}_{ro}^{\top}\left(\mathbf{y}-{\rm
softmax}\left(\mathbf{w}_{ro}\sigma(\boldsymbol{u}_{L})+\mathbf{b}_{ro}\right)\right),\end{array}\right.$

where $\mathbf{w}_{l,l^{\prime}}$ is the weight matrix going from layer
$l^{\prime}$ to $l$, and $\mathbf{b}_{l}$ the bias of layer $l$.
($\mathbf{w}_{ro},\mathbf{b}_{ro}$) is a $10\times 10$ readout layer of the
parameterized cross-entropy loss of [11]. In particular, it does not belong to
the network connections since it does not influence the free dynamics (when
$\beta=0$). The activation $\sigma$ is the shifted sigmoid $x\mapsto
1/(1+e^{-4x+2})$.

#### Convolutional networks.

For the convolutional neural network experiments, the discrete dynamics read:

$\displaystyle\left\\{\begin{array}[]{lll}\displaystyle\boldsymbol{u}_{1}&\leftarrow&\mathcal{P}(\mathbf{w}_{1}^{f}\ast\boldsymbol{x})+\tilde{\mathcal{P}}(\boldsymbol{u}_{1},\mathbf{w}_{2}^{b},\boldsymbol{u}_{2})+\mathbf{b}_{1},\\\\[8.5359pt]
\boldsymbol{u}_{2}&\leftarrow&\mathcal{P}(\mathbf{w}_{2}^{f}\ast\sigma(\boldsymbol{u}_{1}))+\tilde{\mathcal{P}}(\boldsymbol{u}_{2},\mathbf{w}_{3}^{b},\boldsymbol{u}_{3})+\mathbf{b}_{2},\\\\[8.5359pt]
\boldsymbol{u}_{3}&\leftarrow&\mathcal{P}(\mathbf{w}_{3}^{f}\ast\sigma(\boldsymbol{u}_{2}))+\tilde{\mathcal{P}}(\boldsymbol{u}_{3},\mathbf{w}_{4}^{b},\boldsymbol{u}_{4})+\mathbf{b}_{3},\\\\[8.5359pt]
\boldsymbol{u}_{4}&\leftarrow&\mathcal{P}(\mathbf{w}_{4}^{f}\ast\sigma(\boldsymbol{u}_{3}))+\mathbf{b}_{4}+\beta\mathbf{w}_{ro}^{\top}\left(\mathbf{y}-{\rm
softmax}\left(\mathbf{w}_{ro}\sigma(\boldsymbol{u}_{L})+\mathbf{b}_{ro}\right)\right).\end{array}\right.$

Where the superscripts $f$, $b$ respectively mean forward and backward
weights. $\mathcal{P}$ denotes Softmax pooling [78]. The backward
convolutional module going from $\boldsymbol{u}_{l+1}$ to $\boldsymbol{u}_{l}$
is defined as:

$\tilde{\mathcal{P}}(\boldsymbol{u}_{l},\mathbf{w}^{b},\boldsymbol{u}_{l+1}):=\frac{\partial}{\partial\sigma(\boldsymbol{u}_{l})}\left[\boldsymbol{u}_{l+1}\cdot\mathcal{P}(\mathbf{w}^{b}\ast\sigma(\boldsymbol{u}_{l}))\right].$

The activation $\sigma$ is a sigmoid-weighted linear unit [79] defined by:

$\sigma(x):=\left(\frac{x}{2}\right)\frac{1}{1+e^{-x}}+\left(1-\frac{x}{2}\right)\frac{1}{1+e^{-x+2}}.$

The shapes of
$\boldsymbol{u}_{1},\boldsymbol{u}_{2},\boldsymbol{u}_{3},\boldsymbol{u}_{4}$
are respectively $(16,16,128),(8,8,256),(4,4,512),(1,1,512)$. The weights are
all $3\times 3$ kernel with no strides. The paddings are all ‘same’ except for
the last layer. The Softmax pooling has stride 2 and window-size $2\times 2$.
It is defined by [78]:

$y=\sum_{i\in\mathbf{R}}\left(\frac{e^{x_{i}}}{\sum_{j\in\mathbf{R}}e^{x_{j}}}\right)x_{i},$

where $\mathbf{R}$ is the $2\times 2$ window.

###  D.3 Hyperparameters

#### Choice of hyperparameters.

The learning rate was the main hyperparameter searched for by coarse grid
search. The number of time steps was tuned such that the networks have enough
time to equilibrate by measuring the norm of the vector field. The coefficient
for the homeostatic loss was searched by grid search over the range of 0.1,
1.0, and 10.0. Hyperparameters for the CIFAR-10/100 and ImageNet $32\times 32$
were chosen based on those reported in the literature for similar experiments
(for comparison) [11, 7].

Table 3: Hyperparameters used for the Fashion MNIST training experiments. Hyperparameter | Fig. 3, Table 1, Fig. 4e-h) | Fig. 4a-d)
---|---|---
Batch size | 50 | 50
Optimizer | SGD | Adam
Learning rate | 1e-2 | 1e-4
Momentum | 0.9 | –
Epochs | 50 | 50
$\lambda_{\text{homeo}}$ | 0.1 (Fig. 4e-h) | 1.0
$T_{\text{free}}$ | 150 | 150
$T_{\text{nudge}}$ | 20 | 20

In the case of continuous-time estimates (Table 1, Eqs. (8), (9)), we used 60
time steps for each of the $N=4$ $\beta$ values, and five periods, for a total
to 1200 time steps per batch.

Table 4: Hyperparameters used for the CIFAR-10, CIFAR-100 and ImageNet $32\times 32$ training experiments. Hyperparameter | Fig. 4i-l) | Table 2 CIFAR-(10/100) | Table 2 ImageNet $32\times 32$
---|---|---|---
Batch size | 256 | 256 | 256
Optimizer | SGD | SGD | SGD
Learning rate | 5e-3 | 7e-3 | 5e-3
Weight decay | 1e-4 | 1e-4 | 1e-4
Momentum | 0.9 | 0.9 | 0.9
Epochs | 90 | 150 | 90
$\lambda_{\text{homeo}}$ | 10 or 0 | 15 | 10
$T_{\text{free}}$ | 250 | 250 | 250
$T_{\text{nudge}}$ | 40 | 40 | 40
* Learning rates were scaled layer-wise by [25, 15, 10, 8, 5] and decayed with cosine schedule [80].

###  D.4 Hardware and simulation times

Simulations on Fashion MNIST were run on single RTX 5000 GPU, for 10mins per
run and 1 hour for the continuous-time estimate. The convolutional network
simulations were run on an in-house cluster consisting of 5 nodes with 4 v100
NVIDIA GPUs each, one node with 4 A100 NVIDIA GPUs, and one node with 8 A40
NVIDIA GPUs. Runs on CIFAR-10/100 took about 8 hours each, and the Imagenet
$32\times 32$ runs took 72 hours each. Each run was parallelized over 4 or 8
GPUs by splitting the batch. The free phase with the convolutional
architecture was done in half precision, which provided a $2\times$ speed up
compared to using full precision.

  *[BP]: backpropagation of error algorithm
  *[EP]: equilibrium propagation
  *[EBM]: energy-based model
  *[EBMs]: energy-based model
  *[IFT]: implicit function theorem
  *[AD]: automatic differentiation
  *[RBP]: recurrent backpropagation
  *[MLP]: multi-layer perceptron
  *[hEP]: holomorphic EP
  *[TP]: target propagation