# Minimizing Control for Credit Assignment with Strong Feedback

Alexander Meulemans  Matilde Tristany Farinha  Maria Cervera  João Sacramento
Benjamin F. Grewe

###### Abstract

The success of deep learning ignited interest in whether the brain learns
hierarchical representations using gradient-based learning. However, current
biologically plausible methods for gradient-based credit assignment in deep
neural networks need infinitesimally small feedback signals, which is
problematic in biologically realistic noisy environments and at odds with
experimental evidence in neuroscience showing that top-down feedback can
significantly influence neural activity. Building upon deep feedback control
(DFC), a recently proposed credit assignment method, we combine strong
feedback influences on neural activity with gradient-based learning and show
that this naturally leads to a novel view on neural network optimization.
Instead of gradually changing the network weights towards configurations with
low output loss, weight updates gradually minimize the amount of feedback
required from a controller that drives the network to the supervised output
label. Moreover, we show that the use of strong feedback in DFC allows
learning forward and feedback connections simultaneously, using learning rules
fully local in space and time. We complement our theoretical results with
experiments on standard computer-vision benchmarks, showing competitive
performance to backpropagation as well as robustness to noise. Overall, our
work presents a fundamentally novel view of learning as control minimization,
while sidestepping biologically unrealistic assumptions.

Machine Learning, ICML

##  1 Introduction

The error backpropagation (BP) method (Rumelhart et al., 1986; Werbos, 1982)
has emerged as the method of choice for training deep artificial neural
networks, due to its efficient computation of gradients, which is crucial for
learning in high-dimensional parameter spaces (LeCun et al., 2015). Although
biological networks in the neocortex form similar hierarchies (Yamins &
DiCarlo, 2016), it is not clear yet if they follow similar principles to
implement credit assignment (CA), i.e., determining how to change synaptic
connectivity strengths to get closer to some desired output. Several key
properties of the BP method are highly incompatible with cortical networks
(Crick, 1989; Lillicrap et al., 2020). For example, synaptic plasticity in
biological neurons is local in space and time and tightly coupled to neural
activity. By contrast, artificial neural networks trained with BP process
activity and weight update signals separately and in distinct phases while
requiring exactly the same synaptic weights in the forward and feedback
pathways.

To couple synaptic plasticity directly to neural activity, recent work
designed cortical network models that propagate CA signals by leveraging the
network dynamics (Scellier & Bengio, 2017; Whittington & Bogacz, 2017;
Sacramento et al., 2018; Meulemans et al., 2021) or using multiplexed neural
signals (Payeur et al., 2021). However, these methods require tightly
coordinated plasticity mechanisms, possibly through the use of distinct
training phases, and it remains unclear whether cortical networks can exhibit
these levels of coordination. Moreover, to approximate gradient-based
optimization these methods consider the weak-feedback limit where, to avoid
interference, feedback cannot significantly alter the feedforward processing
of the network. This weak-feedback assumption is problematic in biologically
realistic noisy environments, and is in stark contrast with biological
observations showing that top-down feedback can significantly affect forward
processing (Jordan & Keller, 2020; Keller et al., 2020).

To date, several studies used strong feedback control signals in combination
with error-based learning rules to drive the network activity to track a
reference output (Gilra & Gerstner, 2017; Alemi et al., 2018; Denève et al.,
2017; Bourdoukan & Deneve, 2015). Although these models use strong feedback
influences on neural activity, they have only been successfully applied to
train single-layer recurrent neural networks with fixed feedback weights and
it remains unclear whether they can be extended for training deep neural
networks. Furthermore, even though these methods are capable of reducing the
output loss during training, it remains an open question whether they extract
gradient information to enable efficient learning. Hence, gradient-based CA
theories for training deep neural networks with strong feedback influences and
without the need for tightly coordinated plasticity mechanisms are still
lacking.

Building upon Deep Feedback Control (DFC) (Meulemans et al., 2021), we
introduce Strong-DFC, a gradient-based approach to CA in deep neural networks
that makes use of strong feedback influences on neural activity. Like DFC, we
use a feedback controller to drive the network to a desired output target and
update the network weights with a learning rule fully local in space and time
which leverages the dynamic change in neural activity. Unlike DFC, we allow
for strong feedback influences by taking the supervised label as the desired
output target instead of weakly nudging the current network output towards the
supervised label. This new setting results in significant feedback influence
on the neural activities, as the initial network output is far from the
supervised label at the beginning of training.

Interestingly, because the output targets remain invariant throughout
training, Strong-DFC naturally leads to a fundamentally different perspective
on neural network optimization. As the controller drives the network to the
supervised output labels at each point in training, the goal of optimizing the
network parameters shifts from minimizing the supervised loss towards
minimizing the amount of feedback control needed to reach the supervised
output labels. To formalize this minimizing control perspective we introduce a
surrogate loss that quantifies the amount of control feedback, and show that,
under flexible conditions on the feedback connectivity, the Strong-DFC updates
follow its negative gradient. Intriguingly, Strong-DFC enables us to learn the
forward and feedback weights simultaneously, overcoming the two-phase
requirement of DFC, while inheriting the close connection to recently proposed
dendritic compartment models of cortical pyramidal neurons equipped with
voltage-dependent plasticity rules (Sacramento et al., 2018). Finally, we
complement our theoretical findings with experimental results on standard
computer vision benchmarks, showing that Strong-DFC provides principled CA in
deep neural networks using strong feedback influences without the need for
tightly coordinated plasticity mechanisms.

##  2 Background

We revisit DFC (Meulemans et al., 2021) to investigate CA with strong feedback
influence on neural activities. The control framework used by DFC was designed
for feedforward networks with continuous dynamics. More specifically, DFC
considers a multilayer network with the following dynamics (for $1\leq i\leq
L)$:

$\displaystyle\tau_{v}\dot{\mathbf{v}}_{i}(t)$
$\displaystyle=-\mathbf{v}_{i}(t)+W_{i}\phi\big{(}\mathbf{v}_{i-1}(t)\big{)}+Q_{i}\mathbf{u}(t),$
(1)

where $\mathbf{v}_{i}$ is a vector of pre-nonlinearity activations in layer
$i$, $\dot{\mathbf{v}}_{i}$ its time derivative, $W_{i}$ the forward weights,
$\phi$ a smooth nonlinear activation function, $Q_{i}$ the feedback weights,
and $\tau_{v}$ the neural time constant. We define
$\mathbf{r}_{i}=\phi(\mathbf{v}_{i})$ as the post-nonlinearity activity of
layer $i$. Note that without feedback influences ($\mathbf{u}=0$), the steady
state of the above dynamics corresponds to a feedforward neural network, as
the input $\mathbf{r}_{0}$ remains constant over time.

The feedback signal $\mathbf{u}$ is computed by a proportional integral
feedback controller, using the output error
$\mathbf{e}(t)=\mathbf{r}_{L}^{*}-\mathbf{r}_{L}(t)$, and serves to push the
network to a desired output target $\mathbf{r}_{L}^{*}$:

$\displaystyle\mathbf{u}(t)=\mathbf{u}^{\text{int}}(t)+k\mathbf{e}(t),$
$\displaystyle\tau_{u}\dot{\mathbf{u}}^{\text{int}}(t)=\mathbf{e}(t)-\alpha\mathbf{u}^{\text{int}}(t),$
(2)

where $k$ is the proportional control constant and $\alpha$ the leakage
constant limiting the magnitude of $\mathbf{u}^{\mathrm{int}}$. DFC weakly
nudges the network output in the direction of lower loss, since the output
target $\mathbf{r}_{L}^{*}$ is defined as:

$\displaystyle\mathbf{r}_{L}^{*}=\mathbf{r}_{L}^{-}-\lambda\left.\frac{\partial\mathcal{L}}{\partial\mathbf{r}_{L}}\right|_{\mathbf{r}_{L}=\mathbf{r}_{L}^{-}}^{T},$
(3)

with $\mathbf{r}_{L}^{-}$ the steady state of the output without the presence
of a controller, $\mathcal{L}$ the loss function, and $\lambda$ the nudging
strength. The theoretical results linking DFC to a variant of Gauss-Newton
optimization consider the weak-feedback regime of $\lambda\rightarrow 0$
(Meulemans et al., 2021). Note that we cannot enforce strong feedback in this
theoretical framework by keeping the small nudging factor $\lambda\rightarrow
0$ and making $Q_{i}$ arbitrarily large, as the limit considered by Meulemans
et al. (2021) relies on vanishing feedback influence on neural activity,
similar to other seminal work in dynamic CA methods that approximate gradient-
based optimization (Sacramento et al., 2018; Scellier & Bengio, 2017;
Whittington & Bogacz, 2017). DFC then updates the feedforward weights with the
following plasticity rule that is fully local in space and time:

$\displaystyle\tau_{W}\dot{W}_{i}(t)=\big{(}\phi(\mathbf{v}_{i}(t))-\phi(W_{i}\mathbf{r}_{i-1}(t))\big{)}\mathbf{r}_{i-1}(t)^{T}.$
(4)

Figure 1: Strong-DFC trains deep neural networks by minimizing control. (A)
Schematic of a feedforward network trained with Strong-DFC, where feedback
weights $Q$ are uncoupled from forward weights $W$, and used to carry a
control signal $\mathbf{u}$ from the controller to all hidden neurons. (B)
Neurons are multi-compartmental and allow integrating separate forward drive
$v_{i}^{\text{ff}}$ and feedback control $v_{i}^{\text{fb}}$ signals. (C)
Since the supervised target $\mathbf{r}_{L}^{\text{true}}$ is always reached
whenever the controller is on, the training loss $\mathcal{L}$ is equal to
zero throughout, and training can be better understood as the minimization of
a surrogate loss $\mathcal{H}$ describing the amount of control required. (D)
Early in training (top) the controller feedback strongly influences the neural
activity to drive the network output to the desired target, and training
minimizes the amount of control required (bottom). Figure A-B adapted from
Meulemans et al. (2021).

Using a multi-compartment model of a neuron (Fig. 1B), one can interpret this
learning rule as a delta rule that uses the difference between the outgoing
firing rate of the neuron $\phi(\mathbf{v}_{i})$ and the estimated firing rate
based only on the feedforward input of the neuron
$\phi(\mathbf{v}^{\mathrm{ff}}_{i})=\phi(W_{i}\mathbf{v}_{i-1})$ (Urbanczik &
Senn, 2014). Standard DFC learns its feedback weights $Q_{i}$ in a separate
‘sleep’ phase with a simple anti-Hebbian learning rule, making use of the
feedback compartment of the multi-compartment neuron,
$\mathbf{v}^{\mathrm{fb}}_{i}=Q_{i}\mathbf{u}$. Note that DFC assumes a
separation of timescales between the plasticity dynamics (4) and the network
dynamics (1-2), such that the weights $W_{i}$ and $Q_{i}$ can be considered
constant on the timescale of the network-controller interactions.

##  3 Training neural networks with strong feedback by minimizing control

Here, we introduce our minimizing control framework that allows for gradient-
based CA with strong feedback influences in deep neural networks. We use the
feedback controller (2) to strongly drive the network towards the supervised
output label at each point during training. Consequently, the focus of
optimizing the network parameters shifts from minimizing the training loss
$\mathcal{L}$ towards minimizing the amount of feedback control needed to
reach the supervised output label. We quantify the magnitude of the control
input with a surrogate loss $\mathcal{H}$, and show that under flexible
constraints on the feedback connectivity, the Strong-DFC updates perform
gradient descent on this new loss function. Finally, we show that by
minimizing $\mathcal{H}$, we implicitly minimize $\mathcal{L}$ as well.
Throughout this section, we assume the network is stable, on which we
elaborate at the end of the section. Although the theoretical results shown
here are tailored towards the proportional integral controller (2), these can
easily be extended to other controller types such as a pure proportional
controller (App. A.5).

#### Strong-DFC.

We adapt DFC to allow for strong feedback influence on neural activity by
setting the output target $\mathbf{r}_{L}^{*}$ to the supervised output label
$\mathbf{r}_{L}^{\mathrm{true}}$, instead of weakly nudging it towards values
of lower loss (3) (see Fig. 1). This new setting results in strong feedback
influence, as the network output without feedback is far from the supervised
labels at the beginning of training. Furthermore, as the existing theory on
DFC relies on the assumption that the feedback only infinitesimally changes
the network activity, it is no longer valid in this new setting of DFC with
strong feedback (Strong-DFC).

#### Learning as control minimization.

Here, we propose a novel view on neural network optimization that justifies
strong feedback influences on neural activity. In Strong-DFC, at every stage
of training the feedback controller drives the network to the supervised label
$\mathbf{r}_{L}^{\mathrm{true}}$, resulting in zero training loss at steady
state. Hence, instead of gradually changing the network weights towards
configurations with low output loss, the goal of optimizing the network
parameters is now to gradually reduce the amount of feedback required from the
controller to reach the supervised output label (see Fig. 1C-D). At the
individual neuron level, one can interpret this new perspective on
optimization as each neuron minimizing its feedback input by altering its
feedforward weights.

We formalize this notion of minimizing control by introducing a surrogate loss
function representing the magnitude of the feedback control input:

$\displaystyle\mathcal{H}=\frac{1}{2}\sum_{b}\|Q\mathbf{u}^{(b)}_{\mathrm{ss}}\|_{2}^{2},$
(5)

with $Q\triangleq[Q^{T}_{1},...,Q^{T}_{L}]^{T}$ the concatenated feedback
weights and $\mathbf{u}^{(b)}_{\mathrm{ss}}$ the steady-state control feedback
of a data point $b$. The control feedback $\mathbf{u}_{\mathrm{ss}}$ depends
on the feedforward weights $W_{i}$ through the network and controller dynamics
(Eq. 1-2). Using the implicit function theorem, in Theorem 1 we compute the
total derivative of $\mathcal{H}$ w.r.t. $W_{i}$ and show that, under flexible
constraints on the feedback weights, it results in a local update in space and
time that aligns with a simplified version of the Strong-DFC weight updates
(full proof in App. A.1).

######  Theorem 1.

Assuming $J_{\mathrm{ss}}Q+\tilde{\alpha}I$ is invertible, the total
derivative of $\mathcal{H}$ w.r.t. $W_{i}$ is given by:

$\displaystyle\frac{\mathrm{d}\mathcal{H}}{\mathrm{d}W_{i}^{T}}=-\sum_{b}J_{i,\mathrm{ss}}^{(b)T}\left(J_{\mathrm{ss}}^{(b)}Q+\tilde{\alpha}I\right)^{-T}Q^{T}Q\mathbf{u}_{\mathrm{ss}}^{(b)}\mathbf{r}_{i-1,\mathrm{ss}}^{(b)T},$
(6)

with
$J_{i,\mathrm{ss}}\triangleq\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{i}}$
and
$J_{\mathrm{ss}}\triangleq[\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{1}},...,\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{L}}]$,
both evaluated at steady state, and
$\tilde{\alpha}\triangleq\frac{\alpha}{1+k\alpha}$.

Furthermore, assuming Condition 1 holds, the network is stable,
$\alpha\rightarrow 0$, and $J_{\mathrm{ss}}Q+\tilde{\alpha}I$ is invertible,
then the following steady-state (ss) updates for the forward weights, averaged
over the batch samples,

$\displaystyle\Delta
W_{i,\mathrm{ss}}=\eta(\mathbf{v}_{i,\mathrm{ss}}-\mathbf{v}_{i,\mathrm{ss}}^{\mathrm{ff}})\mathbf{r}_{i-1,\mathrm{ss}}^{T}\,\,,$
(7)

with $\eta$ the stepsize, align with the negative gradient
$-\frac{\text{d}\mathcal{H}}{\text{d}W_{i}}$.

######  Condition 1.

The column space of the concatenated feedback weights $Q$ is equal to the row
space of
$J_{\mathrm{ss}}\triangleq[\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{1}},...,\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{L}}]\big{|}_{\mathbf{v}_{i}=\mathbf{v}_{i,\mathrm{ss}}}$,
the network Jacobian at steady state.

Theorem 1 links Strong-DFC to gradient descent on a surrogate loss function
$\mathcal{H}$. Although the steady-state updates (7) of Theorem 1 are not
exactly equal to the DFC updates (4), previous work showed that it is a good
approximation (see App. A.3-A.4 and Meulemans et al., 2021).

#### Flexible feedback connectivity constraint.

Condition 1 is a flexible connectivity constraint on the feedback weights that
represents the requirement that the feedback controller influences the network
activity in the most efficient manner (see App. A.2). More specifically, if
Condition 1 is not satisfied, the feedback input $Q\mathbf{u}_{\mathrm{ss}}$
is partially lost in the nullspace of the network and hence the same influence
on the network output could be reached with a smaller feedback input.
Interestingly, Condition 1 generalizes the feedback requirements of commonly
used learning rules for approximating gradient-based optimization methods. For
example, backpropagation and its approximations require $Q=J^{T}$ (Rumelhart
et al., 1986; Lansdell et al., 2020; Akrout et al., 2019) and common variants
of target propagation require the pseudoinverse $Q=J^{\dagger}$ (Meulemans et
al., 2020), both satisfying Condition 1. More generally, Condition 1 implies
that the feedback weights of Strong-DFC do not necessarily need to align with
the forward pathway or its inverse to provide gradient-based CA by following
the negative gradient of $\mathcal{H}$. Remarkably, the fundamentally distinct
theoretical framework of standard DFC gives rise to a similar condition on
$Q$, however with $J$ evaluated at the feedforward activities without feedback
present instead of the steady-state activities (Meulemans et al., 2021).

#### Minimizing $\mathcal{H}$ leads to minimizing $\mathcal{L}$.

Theorem 1 links Strong-DFC to gradient descent on the surrogate loss function
$\mathcal{H}$. Hence, to show that Strong-DFC minimizes the original loss
function $\mathcal{L}$, we need to connect $\mathcal{H}$ to $\mathcal{L}$.
Proposition 2 shows that if $\mathcal{H}$ is minimized to zero, $\mathcal{L}$
is minimized to zero as well, which is trivial to prove as $\mathcal{H}=0$
translates to all training samples reaching exactly their supervised output
label without help from the controller.

######  Proposition 2.

For all $\mathcal{L}$ that obey $\mathcal{L}(\mathbf{x},\mathbf{x})=0$, we
have that

$\displaystyle\mathcal{H}=\sum_{b}\|Q\mathbf{u}^{(b)}_{\mathrm{ss}}\|_{2}^{2}=0\iff\sum_{b}\mathcal{L}(\mathbf{r}_{L}^{-(b)},\mathbf{r}_{L}^{\mathrm{true,}(b)})=0,$

with $\mathbf{r}_{L}^{-}$ the steady state of the output without the presence
of a controller.

Recent work in the field of deep learning has shown that, for
overparameterized networks, it is possible to minimize the training loss
$\mathcal{L}$ to zero while reaching good generalization (Belkin et al.,
2019). We show empirically in Section 5 that Strong-DFC can minimize
$\mathcal{H}$, and hence $\mathcal{L}$ to zero, thereby validating the
minimizing control framework for optimizing neural networks.

#### Stability of Strong-DFC.

So far, we assumed that the network and controller dynamics reach their steady
state. However, this is not a given, as the feedback interaction could lead to
unstable loops causing the network to diverge. Adapting the results of
Meulemans et al. (2021), we investigate the local stability of a simplified
version of Strong-DFC that uses integral control and assumes a separation of
timescales $\tau_{v}\ll\tau_{u}$. At steady state, the Jacobian matrix of this
dynamical system is equal to $-(J_{\mathrm{ss}}Q+\alpha I)$, resulting in the
following stability condition on the feedback weights (see App. B).

######  Condition 2.

Given the network Jacobian evaluated at steady state $J_{\mathrm{ss}}$, the
real parts of the eigenvalues of $-J_{\mathrm{ss}}Q$ are all below $\alpha$.

In App. B, we compute the Jacobian matrix of the full dynamical system
representing Strong-DFC with proportional integral control and with no
assumptions on a separation of timescales. However, a local stability analysis
on this system does not reveal clearly interpretable stability conditions for
the feedback weights. Therefore, we monitor instead the stability of the full
system empirically.

##  4 A single-phase learning scheme for forward and feedback weights

Conditions 1 and 2 highlight the importance of the feedback weights for
learning and stability, respectively. As training causes the forward weights
and thus the network Jacobian $J_{\mathrm{ss}}$ to change, the set of feedback
connectivity patterns satisfying Conditions 1 and 2 shift during learning as
well. To keep the feedback weights compatible with the changing network, we
learn the feedback weights together with the forward weights.

Here, we introduce plasticity dynamics for the feedback weights that make it
possible to update both forward and feedback weights simultaneously in a
single phase, in contrast to standard DFC and other methods that need two
separate phases. Our new feedback plasticity dynamics leverage two essential
properties of Strong-DFC. (i) The strong feedback influences make it possible
to have noise in the dynamics while learning the forward weights, while
keeping a good signal-to-noise ratio of the learning signals. (ii) The
conditions on the feedback weights for learning and stability required for our
minimizing control framework consider the network activity at steady state.
Hence, the feedback weights should be learned at the same steady state used
for training the forward weights. Leveraging the noise intrinsically available
to neurons, we use a simple anti-Hebbian learning rule for the feedback
weights to drive them towards satisfying Conditions 1 and 2.

#### Noisy neural dynamics.

Our feedback weight plasticity dynamics make use of noise injected in the
neural dynamics, which we model as follows:

$\displaystyle\tau_{v}\dot{\mathbf{v}}_{i}(t)$
$\displaystyle=-\mathbf{v}_{i}(t)+W_{i}\phi\big{(}\mathbf{v}_{i-1}(t)\big{)}+Q_{i}\mathbf{u}(t)+\sigma\boldsymbol{\epsilon}_{i}(t),$
(8)

where the noise $\boldsymbol{\epsilon}$ represents exponentially filtered
white noise with time constant $\tau_{\epsilon}$, i.e., an Ornstein-Uhlenbeck
process, and $\sigma$ the standard deviation of the noise. Note that such
Ornstein-Uhlenbeck processes can, for example, be obtained in the limit of
many independent Poisson inputs with infinitesimal synaptic weights (Gerstner
et al., 2014).

#### An anti-Hebbian learning rule for $Q$.

The noise fluctuations that propagate through the network carry information
about the network Jacobian $J$ and we can extract this information by
correlating the noise fluctuations at the output with the noise fluctuations
at each layer. Inspired by the feedback weight learning of standard DFC
(Meulemans et al., 2021), we instantiate this correlation technique with the
following simple anti-Hebbian plasticity rule:

$\displaystyle\tau_{Q}\dot{Q}_{i}(t)$
$\displaystyle=-\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{L-i}\mathbf{v}^{\mathrm{fb}}_{i}(t)\mathbf{u}_{\mathrm{hp}}(t)^{T}-\beta
Q_{i}.$ (9)

The feedback compartment $\mathbf{v}^{\mathrm{fb}}$ contains a part of the
injected noise,
$\mathbf{v}^{\mathrm{fb}}_{i}=Q_{i}\mathbf{u}+\sigma\boldsymbol{\epsilon}_{i}^{\mathrm{fb}}$,
and the pre-synaptic plasticity signal $\mathbf{u}_{\mathrm{hp}}$ is a high-
pass filtered version of the incoming control signal $\mathbf{u}$, which
extracts the output noise fluctuations.111We high-pass filter $\mathbf{u}(t)$
by subtracting the exponentially averaged signal from it:
$\mathbf{u}_{\mathrm{hp}}(t)\triangleq\mathbf{u}(t)-\mathbf{u}_{\mathrm{lp}}(t)$,
with
$\mathbf{u}_{\mathrm{lp}}\triangleq\frac{1}{\tau_{f}}\int_{0}^{t}\exp\big{(}-\frac{1}{\tau_{f}}(t-t^{\prime})\big{)}\mathbf{u}(t^{\prime})\text{d}t^{\prime}$
and $\tau_{f}$ the filtering time constant. Our feedback weight plasticity
dynamics (9) introduce two important novelties with respect to DFC. (i) The
high-pass filtering of $\mathbf{u}$ makes it possible to learn the feedback
weights using the noise fluctuations of $\mathbf{u}$, while simultaneously
letting the controller push the network towards the output target for learning
the forward weights. (ii) The scaling factor $(1+\tau_{v}/\tau_{\epsilon})$
enables the network to learn an optimal feedback connectivity pattern, while
having a finite neural time constant $\tau_{v}$, which introduces delays in
the network. As the correlation signal decreases for increasing time delays
between the pre- and postsynaptic noise fluctuations, this needs to be
compensated by an appropriate scaling factor.

#### Feedback weights satisfy Conditions 1 and 2.

Theorem 3 shows that under simplifying conditions, the plasticity rule Eq. (9)
drives the feedback connectivity $Q$ to satisfy Conditions 1 and 2 (full proof
in App. C).

######  Theorem 3.

Assume a separation of timescales $\tau_{v},\tau_{\epsilon}\ll\
\tau_{u}\ll\tau_{f}\ll\tau_{Q}$, integral control ($k_{p}=0$), stable network
dynamics, and small noise perturbations $\sigma\ll\lVert\mathbf{v}_{i}\rVert$.
Then, for a fixed data sample, the feedback plasticity dynamics (9) let the
first moment of $Q$ converge approximately towards

$\displaystyle\mathbb{E}[Q_{\mathrm{ss}}]\approx J_{\mathrm{ss}}^{T}M,$ (10)

with $M$ a positive definite symmetric matrix, which satisfies Conditions 1
and 2.

To ensure that the network is stable at all stages, we pre-train the feedback
weights with a large $\alpha$ to guarantee network stability while aligning
the feedback weights with Eq. (10), after which we start the full network
training with a small $\alpha$. Theorem 3 considers training the feedback
weights on a single input sample. In reality, however, multiple input samples
are used for training. For linear networks, the network Jacobian
$J_{\mathrm{ss}}$ is independent of the input samples and Theorem 3 holds
exactly. For nonlinear networks, $J_{\mathrm{ss}}$ will change depending on
the input, hence, the feedback plasticity (9) will drive $Q$ to align with
some average of Eq. (10) for all input samples.

#### Debiasing the forward weight updates.

While we need noise correlations for training the feedback weights, these can
bias the forward weight updates (see App. C.4). We solve this challenge by
low-pass filtering the pre-synaptic plasticity signal for the forward weight
updates, such that the pre-synaptic noise fluctuations are removed:

$\displaystyle\tau_{W}\frac{\text{d}}{\text{d}t}W_{i}(t)$
$\displaystyle=\Big{(}\phi\big{(}\mathbf{v}_{i}(t)\big{)}-\phi\big{(}W_{i}\mathbf{r}_{i-1}(t)\big{)}\Big{)}\mathbf{r}_{\mathrm{lp},i-1}^{T},$
(11)

with
$\mathbf{r}_{\mathrm{lp,}i}\triangleq\frac{1}{\tau_{f}}\int_{0}^{t}\exp\big{(}{-}\frac{1}{\tau_{f}}(t-t^{\prime})\big{)}\mathbf{r}_{i}(t^{\prime})\text{d}t^{\prime}$
an exponential moving average of $\mathbf{r}_{i}$ with filtering time constant
$\tau_{f}$.

##  5 Experiments

Here, we validate our theoretical findings with empirical results. By studying
in detail a synthetic toy example and modest computer-vision benchmarks, we
show that the minimizing control framework succeeds in minimizing the training
loss $\mathcal{L}$, and that Strong-DFC equipped with single-phase learning
succeeds in approximating the gradient of $\mathcal{H}$, thereby validating
itself as a principled CA method. Finally, we show that Strong-DFC is
significantly more noise-robust compared to standard DFC, due to its strong
feedback influences on neural activity.222Source code for all experiments is
available at: https://github.com/mariacer/strong_dfc.

#### Strong-DFC performs CA by minimizing feedback control.

We start with validating our main theoretical claim that minimizing control is
a suitable new theoretical framework for learning. To this end, we first focus
on a low-dimensional teacher-student nonlinear function approximation problem,
which can be efficiently simulated.

In this problem, the goal of learning is to reduce the squared error
$\mathcal{L}=\mathbb{E}_{\mathbf{x}\sim\mathcal{N}(0,I)}[\|f_{\text{trgt}}(\mathbf{x})-\mathbf{v}_{L,\mathrm{ss}}(\mathbf{x})\|^{2}]$
between the target $f_{\text{trgt}}(\mathbf{x})$ provided by the teacher
network we wish to approximate and the steady-state last-layer activity of our
student network, which is a function of the input $\mathbf{x}$. The student
network is trained with Strong-DFC, and its input layer activity
$\mathbf{r}_{0}$ is clamped to the same random pattern
$\mathbf{x}\sim\mathcal{N}(0,I)$ used to probe the teacher function. To
confirm that Theorem 1 and Proposition 2 hold in practice, we study an
idealized setting. We disable neural noise ($\sigma=0$), simulate the neural
dynamics (8) to equilibrium, and then apply a single weight change to $W$
according to the steady-state Strong-DFC update (7). Moreover, to ensure that
Condition 1 is satisfied in this idealized setting at all points during
training, we manually set the feedback weights equal to the transposed network
Jacobian: $Q=\left.J^{T}\right|_{\mathbf{v}=\mathbf{v}(t)}$ (see Condition 1).
We use Euler’s forward method to simulate the dynamics and we initialize the
network state to the feedforward prediction for each datapoint.

Our experiments confirm that the surrogate loss $\mathcal{H}$ is reduced as
training progresses, eventually reaching a low value (of the order of
$10^{-7}$), in accordance with Theorem 1 (see Fig. 2A). Crucially, we also
observe that as $\mathcal{H}$ decreases, so does the actual loss of interest
$\mathcal{L}$. At the end of training, $\mathcal{L}$ is of the order
$10^{-3}$, which supports our claim that the minima of the two loss functions
coincide (cf. Proposition 2). Furthermore, by varying the width of the hidden
layers of the student (Fig. 2B), we show that the amount of control
$\mathcal{H}$ required at the end of training depends on how well the original
training loss $\mathcal{L}$ can be minimized, with more overparameterized
networks reaching lower loss $\mathcal{L}$ and requiring less help
$\mathcal{H}$. We note that such low $\mathcal{L}$ cannot be achieved by a
linear function approximator, nor by a shallow student that only learns its
output weights. This confirms that the ideal version of Strong-DFC performs
useful credit assignment. We conclude that learning by minimizing control
succeeds in solving our nonlinear function approximation problem.

Figure 2: Strong-DFC minimizes the original training loss and the amount of
control. (A) The training loss $\mathcal{L}$ as well as the surrogate loss
$\mathcal{H}$ indicating the amount of control received are minimized during
training with Strong-DFC and ideal feedback weights. The results are obtained
in a nonlinear student-teacher setting with a fixed teacher of size
30-10-10-10-5, and an overparameterized student of size 30-50-50-50-5, both
equipped with $\tanh$ nonlinearities. (B) As hinted by Proposition 2, the
amount of control $\mathcal{H}$ at the end of training depends on how well
$\mathcal{L}$ can be minimized, and is thus dependent on the capacity of the
network. Same setting as in (A), with the loss at convergence plotted in
function of the width of the hidden layers of the student that varies from 3
to 50. We compare the Strong-DFC results against backpropagation (BP) and a
student network of the same architecture but where only the last layer is
trained (BP-shallow).

We confirm the above results for the idealized setting on the MNIST dataset
(LeCun et al., 2010) (Table 1). Here we also observe that both the original
training loss $\mathcal{L}$ and the surrogate loss $\mathcal{H}$ reach low
values, which in turn translates into high testing accuracy. Overall, these
results confirm that the ideal version of Strong-DFC performs useful credit
assignment in more challenging problems.

Table 1: Strong-DFC (with ideal feedback weights) successfully minimizes the original loss $\mathcal{L}$ and the surrogate loss $\mathcal{H}$ in MNIST. Results (mean $\pm$ std, 5 seeds) using a multi-layer perceptron with three hidden layers of 256 neurons. Values are reported for the epoch with lowest training loss out of 40. | Strong-DFC (ideal) | BP
---|---|---
Train loss $\mathcal{L}$ | $8.09^{\pm 0.35}\cdot 10^{-2}$ | $3.38^{\pm 1.30}\cdot 10^{-5}$
Test loss $\mathcal{L}$ | $1.71^{\pm 0.10}\cdot 10^{-1}$ | $9.60^{\pm 0.04}\cdot 10^{-2}$
Train error | $0.07^{\pm 0.09}\%$ | $0.00^{\pm 0.00}\%$
Test error | $1.98^{\pm 0.11}\%$ | $1.78^{\pm 0.04}\%$
Train loss $\mathcal{H}$ | $2.64^{\pm 3.29}\cdot 10^{-5}$ | N/A

#### Feedback learning dynamically tracks feedforward learning.

Figure 3: Strong-DFC simulations for a nonlinear student-teacher regression
task. Standard Strong-DFC, where the feedback weights $Q$ are both pre-trained
before training the forward weights and trained simultaneously with the
forward weights according to Eq. 9 (red), is compared to a set of controls
where the feedback weights are: trained without the
$(1+\tau_{v}/\tau_{\epsilon})^{(L-i)}$ scaling term in the updates (orange),
only pre-trained (dark blue), or kept fixed to a random initialization during
the entire simulation (light blue). The blue vertical dotted lines indicate
the moment where training with random feedback weights leads to
instabilities.(A) Strong-DFC successfully reduces the loss $\mathcal{L}$
during training. (B) When simultaneously learning the feedback weights, the
forward weight updates of Strong-DFC align well with the ideal updates
following the gradient of $\mathcal{H}$, showing that the amount of control is
being minimized with our learning rules. (C) This is achieved through the use
of _strong_ feedback, which allows driving the network to the desired output
even at early stages of training when the ratio of feedback ($\|Qu\|_{2}^{2}$)
to forward ($\|Wr\|_{2}^{2}$) input is high. (D) The single-phase learning
scheme succeeds in driving the feedback weights to satisfying Condition 1,
enabling efficient training of the forward weights. We measure the ratio
between the norm of Q projected into the row space of $J_{\mathrm{ss}}$, and
$Q$ itself, with a value of $1.0$ representing perfect compliance (see App.
E). The results are obtained in a student-teacher setting with a fixed teacher
of size 30-10-10-10-5, and an overparameterized student of size 30-50-50-50-5
with $\tanh$ nonlinearities. Shaded regions indicate standard deviations
across five runs.

Having verified that Strong-DFC fulfills the minimizing control objective, and
that this solves our task, we move to a more challenging setting where
feedback weights are no longer manually fixed to ideal values, and where
neural activity is noisy. The learning problem is otherwise unchanged. In
these experiments, we now adjust the feedback weights $Q$ according to our
anti-Hebbian plasticity rule (9) concurrently as the forward weights evolve
through our always-on, debiased plasticity rule (11). For computational
efficiency, we accumulate our weight changes over time and apply them once the
neural steady-state is reached. We compare Strong-DFC to various ablations on
the feedback weight learning (Fig. 3), and focus on the initial stage of
learning as this provides the clearest differences. We observe that Strong-DFC
successfully decreases the training loss (Fig. 3A) through updates that
approximately align with the negative gradient of $\mathcal{H}$ (Fig. 3B).
Furthermore, our single-phase feedback learning rule succeeds in dynamically
tracking the ongoing feedforward learning process and in finding values of $Q$
which approximately respect Condition 1 (Fig. 3D). We note that throughout our
simulation the feedforward and feedback inputs to any given neuron are of the
same order of magnitude during early training, which confirms that the network
operates in the strong feedback regime (Fig. 3C). Importantly, training
performance, alignment with the gradient of $\mathcal{H}$ and compliance with
Condition 1 are all negatively impacted whenever the feedback weights $Q$ are
not trained during the simulation. Furthermore, without any pre-training,
updates quickly lead to numerical instabilities. This indicates that purely
random feedback mappings are not sufficient for activity-dependent learning
with strong feedback modulation, in contrast to random feedback mappings that
are solely used for plasticity (Lillicrap et al., 2016; Nøkland, 2016).

#### Strong-DFC exhibits robustness to noise and competitive performance in
standard computer vision benchmarks.

We now investigate if Strong-DFC can successfully perform credit assignment on
a deeper neural network model, that we apply to two standard image
classification benchmarks, MNIST and Fashion-MNIST (Xiao et al., 2017).
Moreover, we set the focus of our investigation on noisy neural states, where
$\sigma$ in (8) is non-negligible ($\sigma\sim 10^{-2}$), to determine whether
the strong feedback introduced in Strong-DFC enables learning in the high-
noise regime, and whether our plasticity rule (9) can use the fluctuations to
learn feedback connections that are useful for credit assignment on deeper
networks. In this section, we study a fully-connected three-hidden-layer
neural network (256 units per layer) with a final softmax classification layer
and choose the cross-entropy as our loss function $\mathcal{L}$, which can be
handled in our framework as detailed in Appendix A.6. When using Strong-DFC,
we learn feedback and feedforward weights concurrently with always-on
plasticity rules; plasticity is only interrupted when switching patterns and
re-initializing the network states. To ensure network stability at the
beginning of training, we pre-train the feedback weights using Eq. (9) and a
high controller leak rate $\alpha$.

Table 2: Computer vision results. Test errors (mean $\pm$ std, 5 seeds) obtained after 40 training epochs on a multi-layer perceptron with three hidden layers of 256 neurons each. For DFC, results are provided without noise (original) and with noise during the forward weight training. For Strong-DFC, we consider the (original) single-phase setting with noise and a two-phase setting without noise in the forward weights training phase. MNIST- | standard ($\%$) | fashion ($\%$)
---|---|---
BP | $1.83^{\pm 0.11}$ | $10.60^{\pm 0.44}$
Strong-DFC (with noise) | $2.19^{\pm 0.05}$ | $12.07^{\pm 0.16}$
DFC (with noise) | $15.15^{\pm 0.44}$ | $16.29^{\pm 0.41}$
Strong-DFC (no noise) | $1.98^{\pm 0.05}$ | $11.36^{\pm 0.17}$
DFC (no noise) | $2.09^{\pm 0.10}$ | $11.31^{\pm 0.14}$

In Table 2, we compare the test set prediction error rates of Strong-DFC to
DFC, as well as to standard backpropagation (BP), which we provide as a
reference baseline. We optimize the hyperparameters of each method
independently, for best performance on a validation set of 5000 datapoints. We
use BP to train a noiseless network, and DFC to train both a noiseless and a
noisy network, where we add noise to the dynamics according to Eq. (8). In the
noisy setting, we equip DFC with the debiased plasticity rule (11) for fair
comparison. As previously reported (Meulemans et al., 2021), in the noiseless
case, the test-set performance of DFC is close to that achieved by BP.
However, the situation changes in the large-$\sigma$ regime: the performance
of (weak feedback) DFC drastically drops, whereas Strong-DFC is undisturbed on
MNIST and not severely impaired on Fashion-MNIST (we further illustrate this
in a student-teacher regression setting in Appendix E.2). When we train the
forward weights with Strong-DFC in a separate phase without noise, its
performance is competitive to both DFC and BP. Moreover, we observed that
without feedback plasticity, we did not succeed in training our network when
starting from a random initial $Q$, as the network was unstable. Our feedback
plasticity thus succeeds in pushing the initial network state to an
appropriate one and in keeping it throughout learning.

##  6 Discussion

Modern deep neural networks are invariably learned by backpropagation-of-
error. While backpropagation has proven to be a highly effective gradient-
based credit assignment method for training deep artificial neural networks,
it is unsatisfactory as a model of learning in cortical networks in at least
three important aspects. First, it needs precisely-symmetric forward and
feedback connections. Second, backpropagated errors do not influence neural
activity. Third, it requires two precisely-clocked phases. Instead of
attempting to remedy these issues one by one, we propose a fundamentally new
approach for deep learning, that overcomes all of these by design, and provide
first results indicating that it is an effective method for credit assignment.

Our study shows that these three critical issues can be overcome at once
within the framework of gradient-based optimization, by casting learning as a
control minimization problem. For this, we first augment standard deep neural
networks with a feedback controller, which dynamically changes the hidden
neural activity to reach the desired output. Then, we derive a synaptic
plasticity rule with the objective of minimizing the influence of the feedback
controller on the neural activity. This principled derivation led us to
Strong-DFC, a variation of the recently proposed DFC method which inherits all
of its appealing properties as a biologically plausible deep learning
framework. Like DFC, our weight update rule is local in space and time for a
broad range of (not necessarily symmetric) feedback connectivity patterns and
neural feedback controller architectures, and it can be interpreted as a
biologically inspired dendritic voltage-dependent synaptic plasticity rule
(Urbanczik & Senn, 2014). Importantly, however, top-down feedback signals are
no longer required to be vanishingly small for Strong-DFC.

This strong feedback influence sets our minimizing control framework apart
from previous theories of how the brain might estimate objective function
gradients by measuring small changes in neural activity, generated by
‘nudging’ the network towards a slightly better state (Hinton & McClelland,
1988; O’Reilly, 1996; Xie & Seung, 2003; Scellier & Bengio, 2017; Sacramento
et al., 2018; Lillicrap et al., 2020; Meulemans et al., 2021). In biology,
such theories might face a fundamental obstacle during learning, since
cortical activity is notoriously noisy (Rusakov et al., 2020), making it hard
to estimate small activity changes as plasticity signals. Moreover, they are
at odds with a growing number of experimental neuroscience observations
reporting large neural activity changes in the presence of unexpected or novel
events, which are thought to drive synaptic plasticity (Keller & Mrsic-Flogel,
2018), and the existence of top-down feedback signals that are sufficient to
drive secondary receptive field responses in pyramidal neurons (Keller et al.,
2020). Our minimizing control framework shows how to perform principled
mathematical optimization of an objective function, while being consistent
with these findings.

A long-standing question in neuroscience is whether noise plays a
computational role in the brain – whether it is ‘a bug or a feature’ of
biological neural circuits (Rusakov et al., 2020). Our study suggests that
noise could play a key role in ensuring that gradient-based credit assignment
information is transmitted backwards through the network. Building on the
noise-tolerance of Strong-DFC, we leverage noise to concurrently and
seamlessly learn both feedforward and feedback weights so that Condition 1 is
constantly fulfilled. In contrast to standard two-phase learning of forward
and backward connections (Akrout et al., 2019; Lee et al., 2015; Meulemans et
al., 2021), our plasticity model does not require any pauses, distinct phases,
plasticity switches, or alternations between noise- and input-driven modes.
Our always-on plasticity rules are arguably more biologically plausible than
alternative two-phase learning algorithms, and simpler to implement in
neuromorphic hardware.

In our pyramidal cell model, top-down apical synapses transmit feedback
information which instructs the plasticity of bottom-up basal synapses, an
idea that can be traced back to the seminal work of Körding & König (2001).
Interpreted as such, our results suggest that the plasticity rules for top-
down and bottom-up connections might display major differences, which is
experimentally testable. Both rules include a voltage-dependent postsynaptic
factor, in agreement with phenomenological plasticity models (Clopath &
Gerstner, 2010). However, the presynaptic components of our apical and basal
plasticity rules are high- and low-pass filtered, respectively. Such
presynaptic activity filters are a recurring element in reward-driven synaptic
plasticity rules (Seung, 2003). Moreover, our apical plasticity rule is anti-
Hebbian and its postsynaptic term depends only on the local dendritic voltage,
not on the somatic spiking activity. This may be consistent with the fact
that, for pyramidal cells, backpropagating action potentials often fail to
reach the apical dendrite, questioning their role as the main drivers of
apical synaptic plasticity (Spruston, 2008; Gambino et al., 2014).

We highlight that, despite the close methodological connections between DFC
and Strong-DFC, our minimizing control framework is fundamentally different
from the underlying theory of DFC and many other recent biologically plausible
learning methods that link their weight updates to a new variant of Gauss-
Newton optimization with a minibatch size of one (Meulemans et al., 2021,
2020; Podlaski & Machens, 2020; Bengio, 2020). While the specific optimization
characteristics of these methods are not yet fully understood, our minimizing
control framework can be cast as gradient descent on a surrogate loss
function. Hence, the rich field studying the application of stochastic
gradient descent to neural networks (Bottou et al., 2018) can be applied on
our new loss to uncover the optimization characteristics of Strong-DFC.
Furthermore, our minimizing control approach draws from control theory, but it
is distinct from standard optimal control. In particular, it is well known
that the adjoint state method from optimal control is equivalent to
backpropagation on the supervised loss function, when applied to deep
feedforward neural networks (LeCun et al., 1988), and hence different from our
minimizing control framework that performs gradient descent on a surrogate
loss function.

In practice, the Strong-DFC updates do not exactly follow the negative
gradient of $\mathcal{H}$, due to (i) the limited number of iterations for
training the feedback weights, (ii) the limited capacity of the linear pathway
$Q$ to satisfy Condition 1 for every input sample, and (iii) the difference
between the steady-state weight update (7) and the Strong-DFC weight update
(4). Despite these challenges, Fig 3 shows that Strong-DFC approximates the
theory well. Future work can further improve Strong-DFC by investigating
alternative feedback mappings and network architectures. As Strong-DFC uses
continuous dynamics, it is costly to simulate on conventional deep learning
hardware, preventing it from being tested on large-scale machine learning
benchmarks. A promising solution is to implement Strong-DFC on analog and
neuromorphic hardware, where the dynamics can coincide with the physical
dynamics of the analog components. Consequently, this can make Strong-DFC an
attractive principled approach for CA on analog deep learning implementations,
commonly used in Edge-AI and other low-energy applications of deep learning
(Xiao et al., 2020; Misra & Saha, 2010).

By providing a novel optimization approach for deep neural networks that
differs from the usual direct minimization of the output loss, our minimizing
control framework allowed us to derive a novel gradient-based credit
assignment method that relies on strong feedback influence. Instantiating our
framework with Strong-DFC naturally leads to biologically desired
characteristics, such as single-phase learning, flexible feedback connectivity
requirements, and local learning rules, further underlining its promise for
understanding credit assignment in the brain.

## Acknowledgements

This work was supported by the Swiss National Science Foundation (B.F.G.
CRSII5-173721 and 315230_189251), ETH project funding (B.F.G. ETH-20 19-01)
and the Human Frontiers Science Program (RGY0072/2019). João Sacramento was
supported by an Ambizione grant (PZ00P3_186027) from the Swiss National
Science Foundation. We would like to thank Nicolas Zucchet for his valuable
feedback and engaging discussions on the optimization theory aspects of the
minimizing control framework; and Aditya Gilra for insightful discussions on
the learning through control framework.

## References

  * Akrout et al. (2019) Akrout, M., Wilson, C., Humphreys, P., Lillicrap, T., and Tweed, D. B.  Deep learning without weight transport.  In _Advances in Neural Information Processing Systems 32_ , pp. 974–982, 2019.
  * Alemi et al. (2018) Alemi, A., Machens, C., Denève, S., and Slotine, J.-J.  Learning arbitrary dynamics in efficient, balanced spiking networks using local plasticity rules.  _AAAI Conference on Artificial Intelligence (AAAI)_ , 2018.
  * Belkin et al. (2019) Belkin, M., Hsu, D., Ma, S., and Mandal, S.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  _Proceedings of the National Academy of Sciences_ , 116(32):15849–15854, 2019.
  * Bengio (2020) Bengio, Y.  Deriving differential target propagation from iterating approximate inverses.  _arXiv preprint arXiv:2007.15139_ , 2020.
  * Bottou et al. (2018) Bottou, L., Curtis, F. E., and Nocedal, J.  Optimization methods for large-scale machine learning.  _SIAM Review_ , 60(2):223–311, 2018.
  * Bourdoukan & Deneve (2015) Bourdoukan, R. and Deneve, S.  Enforcing balance allows local supervised learning in spiking recurrent networks.  _Advances in Neural Information Processing Systems_ , 28:982–990, 2015.
  * Clopath & Gerstner (2010) Clopath, C. and Gerstner, W.  Voltage and spike timing interact in stdp–a unified model.  _Frontiers in synaptic neuroscience_ , 2:25, 2010.
  * Crick (1989) Crick, F.  The recent excitement about neural networks.  _Nature_ , 337(6203):129–132, 1989.
  * Denève et al. (2017) Denève, S., Alemi, A., and Bourdoukan, R.  The brain as an efficient and robust adaptive learner.  _Neuron_ , 94(5):969–977, 2017.
  * Gambino et al. (2014) Gambino, F., Kehayas, V., Baptista, D., Tatti, R., Carleton, A., and Holtmaat, A.  Sensory-evoked LTP driven by dendritic plateau potentials in vivo.  _Nature_ , 515(7525):116–119, 2014.
  * Gerstner et al. (2014) Gerstner, W., Kistler, W. M., Naud, R., and Paninski, L.  _Neuronal dynamics: From single neurons to networks and models of cognition_.  Cambridge University Press, 2014.
  * Gilra & Gerstner (2017) Gilra, A. and Gerstner, W.  Predicting non-linear dynamics by stable local learning in a recurrent spiking neural network.  _Elife_ , 6:e28295, 2017.
  * Hinton & McClelland (1988) Hinton, G. E. and McClelland, J. L.  Learning representations by recirculation.  In _Advances in Neural Information Processing Systems_ , pp. 358–366, 1988.
  * Jordan & Keller (2020) Jordan, R. and Keller, G. B.  Opposing influence of top-down and bottom-up input on excitatory layer 2/3 neurons in mouse primary visual cortex.  _Neuron_ , 108(6):1194–1206, 2020.
  * Keller et al. (2020) Keller, A. J., Roth, M. M., and Scanziani, M.  Feedback generates a second receptive field in neurons of the visual cortex.  _Nature_ , 582(7813):545–549, 2020.
  * Keller & Mrsic-Flogel (2018) Keller, G. B. and Mrsic-Flogel, T. D.  Predictive processing: a canonical cortical computation.  _Neuron_ , 100(2):424–435, 2018.
  * Körding & König (2001) Körding, K. P. and König, P.  Supervised and unsupervised learning with two sites of synaptic integration.  _Journal of Computational Neuroscience_ , 11(3):207–215, 2001.
  * Lansdell et al. (2020) Lansdell, B. J., Prakash, P., and Kording, K. P.  Learning to solve the credit assignment problem.  In _International Conference on Learning Representations_ , 2020.
  * LeCun et al. (1988) LeCun, Y., Touresky, D., Hinton, G., and Sejnowski, T.  A theoretical framework for back-propagation.  In _Proceedings of the 1988 Connectionist Models Summer School_ , volume 1, pp. 21–28. Morgan Kaufmann, 1988.
  * LeCun et al. (2010) LeCun, Y., Cortes, C., and Burges, C.  MNIST handwritten digit database.  _AT &T Labs. Available: http://yann.lecun.com/exdb/mnist_, 2010.
  * LeCun et al. (2015) LeCun, Y., Bengio, Y., and Hinton, G.  Deep learning.  _Nature_ , 521(7553):436, 2015.
  * Lee et al. (2015) Lee, D.-H., Zhang, S., Fischer, A., and Bengio, Y.  Difference target propagation.  In _Joint european conference on machine learning and knowledge discovery in databases_ , pp. 498–515. Springer, 2015.
  * Lillicrap et al. (2016) Lillicrap, T. P., Cownden, D., Tweed, D. B., and Akerman, C. J.  Random synaptic feedback weights support error backpropagation for deep learning.  _Nature Communications_ , 7:13276, 2016.
  * Lillicrap et al. (2020) Lillicrap, T. P., Santoro, A., Marris, L., Akerman, C. J., and Hinton, G.  Backpropagation and the brain.  _Nature Reviews Neuroscience_ , pp. 1–12, 2020.
  * Lu & Shiou (2002) Lu, T.-T. and Shiou, S.-H.  Inverses of 2$\times$ 2 block matrices.  _Computers & Mathematics with Applications_, 43(1-2):119–129, 2002.
  * Meulemans et al. (2020) Meulemans, A., Carzaniga, F., Suykens, J., Sacramento, J. a., and Grewe, B. F.  A theoretical framework for target propagation.  _Advances in Neural Information Processing Systems_ , 33:20024–20036, 2020.
  * Meulemans et al. (2021) Meulemans, A., Farinha, M. T., Ordóñez, J. G., Aceituno, P. V., Sacramento, J., and Grewe, B. F.  Credit assignment in neural networks through deep feedback control.  _arXiv preprint arXiv:2106.07887_ , 2021.
  * Misra & Saha (2010) Misra, J. and Saha, I.  Artificial neural networks in hardware: A survey of two decades of progress.  _Neurocomputing_ , 74(1-3):239–255, 2010.
  * Nøkland (2016) Nøkland, A.  Direct feedback alignment provides learning in deep neural networks.  In _Advances in neural information processing systems_ , pp. 1037–1045, 2016.
  * O’Reilly (1996) O’Reilly, R. C.  Biologically plausible error-driven learning using local activation differences: The generalized recirculation algorithm.  _Neural Computation_ , 8(5):895–938, 1996.
  * Payeur et al. (2021) Payeur, A., Guerguiev, J., Zenke, F., Richards, B., and Naud, R.  Burst-dependent synaptic plasticity can coordinate learning in hierarchical circuits.  _Nature neuroscience_ , 24(5):1546, 2021.
  * Podlaski & Machens (2020) Podlaski, W. F. and Machens, C. K.  Biological credit assignment through dynamic inversion of feedforward networks.  _Advances in Neural Information Processing Systems 33_ , 2020.
  * Rumelhart et al. (1986) Rumelhart, D. E., Hinton, G. E., and Williams, R. J.  Learning representations by back-propagating errors.  _Nature_ , 323(6088):533, 1986.
  * Rusakov et al. (2020) Rusakov, D. A., Savtchenko, L. P., and Latham, P. E.  Noisy synaptic conductance: bug or a feature?  _Trends in Neurosciences_ , 43(6):363–372, 2020\.
  * Sacramento et al. (2018) Sacramento, J., Costa, R. P., Bengio, Y., and Senn, W.  Dendritic cortical microcircuits approximate the backpropagation algorithm.  In _Advances in Neural Information Processing Systems 31_ , pp. 8721–8732, 2018.
  * Särkkä & Solin (2019) Särkkä, S. and Solin, A.  _Applied stochastic differential equations_ , volume 10.  Cambridge University Press, 2019.
  * Scellier & Bengio (2017) Scellier, B. and Bengio, Y.  Equilibrium propagation: Bridging the gap between energy-based models and backpropagation.  _Frontiers in computational neuroscience_ , 11:24, 2017\.
  * Seung (2003) Seung, H. S.  Learning in spiking neural networks by reinforcement of stochastic synaptic transmission.  _Neuron_ , 40(6):1063–1073, 2003.
  * Spruston (2008) Spruston, N.  Pyramidal neurons: dendritic structure and synaptic integration.  _Nature Reviews Neuroscience_ , 9(3):206–221, 2008.
  * Urbanczik & Senn (2014) Urbanczik, R. and Senn, W.  Learning by the dendritic prediction of somatic spiking.  _Neuron_ , 81(3):521–528, 2014.
  * Werbos (1982) Werbos, P. J.  Applications of advances in nonlinear sensitivity analysis.  In _System modeling and optimization_ , pp. 762–770. Springer, 1982\.
  * Whittington & Bogacz (2017) Whittington, J. C. and Bogacz, R.  An approximation of the error backpropagation algorithm in a predictive coding network with local hebbian synaptic plasticity.  _Neural computation_ , 29(5):1229–1262, 2017\.
  * Xiao et al. (2017) Xiao, H., Rasul, K., and Vollgraf, R.  Fashion-MNIST: a novel image dataset for benchmarking machine learning algorithms.  _arXiv preprint arXiv:1708.07747_ , 2017.
  * Xiao et al. (2020) Xiao, T. P., Bennett, C. H., Feinberg, B., Agarwal, S., and Marinella, M. J.  Analog architectures for neural network acceleration based on non-volatile memory.  _Applied Physics Reviews_ , 7(3):031301, 2020\.
  * Xie & Seung (2003) Xie, X. and Seung, H. S.  Equivalence of backpropagation and contrastive Hebbian learning in a layered network.  _Neural Computation_ , 15(2):441–454, 2003.
  * Yamins & DiCarlo (2016) Yamins, D. L. and DiCarlo, J. J.  Using goal-driven deep learning models to understand sensory cortex.  _Nature neuroscience_ , 19(3):356–365, 2016.

##  Appendix A Proofs and discussion for Section 3

###  A.1 Proof Theorem 1

In this appendix section, we prove Theorem 1, which we copy below for
convenience.

######  Theorem S1.

Assuming $J_{\mathrm{ss}}Q+\tilde{\alpha}I$ is invertible, the total
derivative of $\mathcal{H}$ w.r.t. $W_{i}$ is given by:

$\displaystyle\frac{\mathrm{d}\mathcal{H}}{\mathrm{d}W_{i}}=-\sum_{b}J_{i,\mathrm{ss}}^{(b)T}\left(J_{\mathrm{ss}}^{(b)}Q+\tilde{\alpha}I\right)^{-T}Q^{T}Q\mathbf{u}_{\mathrm{ss}}^{(b)}\mathbf{r}_{i-1,\mathrm{ss}}^{(b)T},$
(12)

with
$J_{i,\mathrm{ss}}\triangleq\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{i}}$
and
$J_{\mathrm{ss}}\triangleq[\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{1}},...,\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{L}}]$,
both evaluated at steady state, and
$\tilde{\alpha}\triangleq\frac{\alpha}{1+k\alpha}$.

Furthermore, assuming Condition 1 holds, the network is stable,
$\alpha\rightarrow 0$, and $J_{\mathrm{ss}}Q+\tilde{\alpha}I$ is invertible,
then the following steady-state (ss) updates for the forward weights, averaged
over the batch samples,

$\displaystyle\Delta
W_{i,\mathrm{ss}}=\eta(\mathbf{v}_{i,\mathrm{ss}}-\mathbf{v}_{i,\mathrm{ss}}^{\mathrm{ff}})\mathbf{r}_{i-1,\mathrm{ss}}^{T}\,\,,$
(13)

with $\eta$ the stepsize, align with the negative gradient
$-\frac{\text{d}\mathcal{H}}{\text{d}W_{i}}$.

######  Condition S1.

The column space of $Q$ is equal to the row space of
$J_{\mathrm{ss}}\triangleq[\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{1}},...,\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{L}}]$
at steady state.

The intuition behind the proof is as follows. First, we compute in Lemma S2
the gradient of $\mathcal{H}$ using the Implicit Function Theorem, as this
provides us with the means to calculate how the control signal
$\mathbf{u}_{\mathrm{ss}}$ at steady steady state changes when $W_{i}$ is
changed. In Lemma S3, we show that $\frac{\text{d}\mathcal{H}}{\text{d}W_{i}}$
collapses to a local update if Condition 1 is satisfied. Finally, we prove
Theorem 1 by showing that this local update is equal to the steady-state
weight update in Eq. (7).

######  Lemma S2.

The gradient of $\mathcal{H}$ w.r.t. $\boldsymbol{\theta}$ for a single
sample, evaluated at steady state, is given by

$\displaystyle\frac{1}{2}\frac{\mathrm{d}\|Q\mathbf{u}_{\mathrm{ss}}\|^{2}_{2}}{\mathrm{d}\boldsymbol{\theta}}=-\mathbf{u}_{\mathrm{ss}}^{T}Q^{T}Q\left(J_{\mathrm{ss}}Q+\tilde{\alpha}I\right)^{-1}J_{\mathrm{ss}}R_{\mathrm{ss}}^{T},$
(14)

with $\boldsymbol{\theta}=[\vec{W}_{1}^{T},...\vec{W}_{L}^{T}]^{T}$ and
$R_{\mathrm{ss}}$ defined below.

$\displaystyle
R_{\mathrm{ss}}^{T}\triangleq\begin{bmatrix}\mathbf{r}_{0}^{T}\otimes I&0&0\\\
0&\ddots&0\\\ 0&0&\mathbf{r}_{L-1}^{T}\otimes I\end{bmatrix}$ (15)

with $\otimes$ the Kronecker product.

###### Proof.

Using the chain rule, we get the following expression for the gradient of
interest:

$\displaystyle\frac{1}{2}\frac{\mathrm{d}\|Q\mathbf{u}_{\mathrm{ss}}\|^{2}_{2}}{\mathrm{d}\boldsymbol{\theta}}=\mathbf{u}_{\mathrm{ss}}^{T}Q^{T}Q\frac{\mathrm{d}\mathbf{u}_{\mathrm{ss}}}{\mathrm{d}\boldsymbol{\theta}}$
(16)

Now we proceed to find
$\frac{\mathrm{d}\mathbf{u}_{\mathrm{ss}}}{\mathrm{d}\boldsymbol{\theta}}$. At
steady state, the network and controller dynamics result in the following
equilibrium equations:

$\displaystyle\mathbf{v}_{i,\mathrm{ss}}$
$\displaystyle=W_{i}\phi(\mathbf{v}_{i-1,\mathrm{ss}})+Q_{i}\mathbf{u}_{\mathrm{ss}},\quad
1\leq i\leq L$ (17) $\displaystyle\mathbf{u}_{\mathrm{ss}}$
$\displaystyle=k\mathbf{e}_{\mathrm{ss}}+\mathbf{u}^{\mathrm{int}}_{\mathrm{ss}}$
(18) $\displaystyle 0$
$\displaystyle=\mathbf{e}_{\mathrm{ss}}-\alpha\mathbf{u}^{\mathrm{int}}_{\mathrm{ss}}$
(19) $\displaystyle\mathbf{e}_{\mathrm{ss}}$
$\displaystyle=\mathbf{r}_{L}^{*}-\mathbf{r}_{L,\mathrm{ss}}$ (20)

The equilibrium solutions in Eq. (17) can be used to obtain an explicit
function of $\mathbf{r}_{L}$ w.r.t. $\boldsymbol{\theta}$ and
$\mathbf{u}_{\mathrm{ss}}$, which in turn can be used to obtain an explicit
function of $\mathbf{e}_{\mathrm{ss}}$ i.f.o. $\boldsymbol{\theta}$ and
$\mathbf{u}_{\mathrm{ss}}$. Combining all equations results in the following
implicit function relating $\mathbf{u}_{\mathrm{ss}}$ to
$\boldsymbol{\theta}$:

$\displaystyle
F(\mathbf{u}_{\mathrm{ss}},\boldsymbol{\theta})\triangleq\mathbf{e}_{\mathrm{ss}}(\mathbf{u}_{\mathrm{ss}},\boldsymbol{\theta})-\tilde{\alpha}\mathbf{u}_{\mathrm{ss}}=0,$
(21)

with $\tilde{\alpha}\triangleq\frac{\alpha}{1+\alpha k}$. If $\frac{\partial
F}{\partial\mathbf{u}_{\mathrm{ss}}}$ is invertible, we can use the Implicit
Function Theorem to obtain the total derivative of $\mathbf{u}_{\mathrm{ss}}$
w.r.t. $\boldsymbol{\theta}$:

$\displaystyle\frac{\mathrm{d}\mathbf{u}_{\mathrm{ss}}}{\mathrm{d}\boldsymbol{\theta}}=-\left(\frac{\partial
F}{\partial\mathbf{u}_{\mathrm{ss}}}\right)^{-1}\frac{\partial
F}{\partial\boldsymbol{\theta}}.$ (22)

Investigating each part separately, we get

$\displaystyle\frac{\partial F}{\partial\mathbf{u}_{\mathrm{ss}}}$
$\displaystyle=-\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{u}}-\tilde{\alpha}I$
(23)
$\displaystyle=-\sum_{i=1}^{L}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{i}}\frac{\partial\mathbf{v}_{i}}{\partial\mathbf{u}}-\tilde{\alpha}I$
(24) $\displaystyle=-J_{\mathrm{ss}}Q-\tilde{\alpha}I,$ (25)

with all derivatives evaluated at the steady-state activations, and

$\displaystyle\frac{\partial F}{\partial\boldsymbol{\theta}}$
$\displaystyle=\frac{\partial\mathbf{r}_{L}}{\partial\boldsymbol{\theta}}$
(26) $\displaystyle=-J_{\mathrm{ss}}R_{\mathrm{ss}}^{T},$ (27)

with $R_{\mathrm{ss}}^{T}$ defined in Eq. (15). Bringing everything together
results in

$\displaystyle\frac{\mathrm{d}\mathbf{u}_{\mathrm{ss}}}{\mathrm{d}\boldsymbol{\theta}}$
$\displaystyle=-\left(J_{\mathrm{ss}}Q+\tilde{\alpha}I\right)^{-1}J_{\mathrm{ss}}R_{\mathrm{ss}}^{T}$
(28)
$\displaystyle\frac{1}{2}\frac{\mathrm{d}\|Q\mathbf{u}_{\mathrm{ss}}\|^{2}_{2}}{\mathrm{d}\boldsymbol{\theta}}$
$\displaystyle=-\mathbf{u}_{\mathrm{ss}}^{T}Q^{T}Q\left(J_{\mathrm{ss}}Q+\tilde{\alpha}I\right)^{-1}J_{\mathrm{ss}}R_{\mathrm{ss}}^{T}.$
(29)

∎

######  Lemma S3.

When Condition 1 is satisfied and $J_{\mathrm{ss}}Q$ is invertible,
$\frac{\mathrm{d}\|Q\mathbf{u}_{\mathrm{ss}}\|^{2}_{2}}{\mathrm{d}W_{i}}$
results in a local update rule in the limit of $\tilde{\alpha}\rightarrow 0$:

$\displaystyle\lim_{\tilde{\alpha}\rightarrow
0}\frac{1}{2}\frac{\mathrm{d}\|Q\mathbf{u}_{\mathrm{ss}}\|^{2}_{2}}{\mathrm{d}W_{i}}=-Q_{i}\mathbf{u}_{\mathrm{ss}}\mathbf{r}_{i-1,\mathrm{ss}}^{T}.$
(30)

###### Proof.

Using a similar proof as Lemma S2 in Meulemans et al. (2021), we have that

$\displaystyle\lim_{\tilde{\alpha}\rightarrow
0}\left(J_{\mathrm{ss}}Q+\tilde{\alpha}I\right)^{-1}J_{\mathrm{ss}}=Q^{\dagger}$
(31)

iff $\mathrm{Row}(J_{\mathrm{ss}})=\mathrm{Col}(Q)$, with $Q^{\dagger}$ the
Moore-Penrose pseudoinverse of $Q$. Using the fact that
$Q^{T}QQ^{\dagger}=Q^{T}$, we have that

$\displaystyle\lim_{\tilde{\alpha}\rightarrow
0}\frac{1}{2}\frac{\mathrm{d}\|Q\mathbf{u}_{\mathrm{ss}}\|^{2}_{2}}{\mathrm{d}\boldsymbol{\theta}}$
$\displaystyle=-\mathbf{u}_{\mathrm{ss}}^{T}Q^{T}R_{\mathrm{ss}}^{T}$ (32)

iff $\mathrm{Row}(J_{\mathrm{ss}})=\mathrm{Col}(Q)$ and $J_{\mathrm{ss}}$ is
of full row rank. Going from the vectorized notation for $\boldsymbol{\theta}$
towards the matrix notation for $W_{i}$ gives us

$\displaystyle\lim_{\tilde{\alpha}\rightarrow
0}\frac{1}{2}\frac{\mathrm{d}\|Q\mathbf{u}_{\mathrm{ss}}\|^{2}_{2}}{\mathrm{d}W_{i}^{T}}=-Q\mathbf{u}_{\mathrm{ss}}\mathbf{r}_{i-1,\mathrm{ss}}^{T}.$
(33)

∎

Now, we are ready to prove Theorem 1.

###### Proof.

From the linearity of the derivative operator, we have that

$\displaystyle\frac{\mathrm{d}\mathcal{H}}{\mathrm{d}W_{i}^{T}}=\sum_{b}\frac{1}{2}\frac{\mathrm{d}\|Q\mathbf{u}_{\mathrm{ss}}^{(b)}\|^{2}_{2}}{\mathrm{d}W_{i}}.$
(34)

Rewriting Lemma S2 from its vectorized form towards matrix form, and filling
it in the above equation, proves the first part of Theorem 1:

$\displaystyle\frac{\mathrm{d}\mathcal{H}}{\mathrm{d}W_{i}^{T}}=-\sum_{b}J_{i,\mathrm{ss}}^{(b)T}\left(J_{\mathrm{ss}}^{(b)}Q+\tilde{\alpha}I\right)^{-T}Q^{T}Q\mathbf{u}_{\mathrm{ss}}^{(b)}\mathbf{r}_{i-1,\mathrm{ss}}^{(b)T},$
(35)

For the second part of the Theorem, the conditions of Lemma S3 are satisfied
by the stated assumptions of this theorem. Hence, we have that:

$\displaystyle\lim_{\tilde{\alpha}\rightarrow
0}\frac{\mathrm{d}\mathcal{H}}{\mathrm{d}W_{i}^{T}}=-\sum_{b}Q_{i}\mathbf{u}_{\mathrm{ss}}^{(b)}\mathbf{r}_{i-1,\mathrm{ss}}^{(b)T}.$
(36)

Finally, from the steady-state solution for the network dynamics, we have that
$\mathbf{v}_{i,\mathrm{ss}}-\mathbf{v}^{\mathrm{fb}}_{i,\mathrm{ss}}=Q_{i}\mathbf{u}_{\mathrm{ss}}$,
thereby concluding the proof. ∎

###  A.2 Interpreting Condition 1

Here, we investigate in more depth why Condition 1 is needed in Theorem 1. The
main line of argumentation goes as follows: if Condition 1 is not satisfied,
the steady-state weight update given by Eq. (7) has components inside the
nullspace of the network, and hence, can only partially influence the network
output. The gradient of $\|Q\mathbf{u}_{\mathrm{ss}}\|_{2}^{2}$ is its
steepest ascent direction, thus all the components of this gradient should
influence the network output and consequently, the amount of needed control
feedback, because otherwise there would exist a ’steeper’ direction.
Therefore, if the weight update given by Eq. (7) contains nullspace
components, it cannot be fully aligned with the gradient.

Now, let us specify the nullspace of the network. As the dimension of the
concatenated forward weights $\boldsymbol{\theta}$ is higher than the
dimension of the network output, there exist certain weight updates
$\Delta\boldsymbol{\theta}$ that do not result in any change to the output
$\mathbf{r}_{L}$ of the network, as they lie inside the nullspace of the
network. For small weight updates, the nullspace of the network at steady
state is connected to the network Jacobian, as we have that
$\Delta\mathbf{r}_{L}\approx
J_{\mathrm{ss}}R_{\mathrm{ss}}^{T}\Delta\boldsymbol{\theta}$, with
$R_{\mathrm{ss}}$ as defined in Lemma S2. From this formulation, one can see
that the nullspace of the network is equal to the nullspace of
$J_{\mathrm{ss}}R_{\mathrm{ss}}^{T}$.

When vectorizing the steady-state weight updates (Eq. 7), we see that they lie
inside the rowspace of $Q^{T}R_{\mathrm{ss}}^{T}$. If we have that the column
space of $Q$ (hence, the rowspace of $Q^{T}$) is equal to the row space of
$J_{\mathrm{ss}}$, the steady-state weight updates will always lie inside the
rowspace of $J_{\mathrm{ss}}R_{\mathrm{ss}}^{T}$, and hence, can never lie
inside the nullspace of the network.

The negative gradient of $\|Q\mathbf{u}_{\mathrm{ss}}\|_{2}^{2}$ is its
steepest descent direction, thus all the components of the gradient should
influence the network output and consequently the amount of needed control
input. Therefore, if the weight update Eq. (7) contains nullspace components,
it cannot be fully aligned with the gradient. Consequently, Condition 1 is a
necessary condition for the steady-state update (Eq. 7) to align with the
negative gradient. The sufficiency of Condition 1 follows from the proof of
Theorem 1.

To illustrate the above arguments further, let us investigate in more detail
what happens when we update the forward network parameters
$\boldsymbol{\theta}$. Assuming a small weight update
$\Delta\boldsymbol{\theta}$, we can do a first-order Taylor expansion of the
steady-state network output after the weight update:

$\displaystyle\mathbf{r}_{L,\mathrm{ss}}^{(m+1)}\approx\mathbf{r}_{L,\mathrm{ss}}^{(m)}+J_{\boldsymbol{\theta}}\Delta\boldsymbol{\theta}+J_{\mathrm{ss}}Q\Delta\mathbf{u}_{\mathrm{ss}},$
(37)

with
$J_{\boldsymbol{\theta}}=\frac{\partial\mathbf{r}_{L}}{\partial\boldsymbol{\theta}}=J_{\mathrm{ss}}R_{\mathrm{ss}}^{T}$
and $m$ the training iteration. The weight update changes the amount of
control needed to reach the output target, which we capture by introducing
$\Delta\mathbf{u}_{\mathrm{ss}}$. As the feedback controller drives the
network to the supervised label during each training iteration (assuming
$\alpha\rightarrow 0$), we have that
$\mathbf{r}_{L,\mathrm{ss}}^{(m+1)}=\mathbf{r}_{L,\mathrm{ss}}^{(m)}=\mathbf{r}_{L}^{\mathrm{true}}$
and hence

$\displaystyle
J_{\mathrm{ss}}R_{\mathrm{ss}}^{T}\Delta\boldsymbol{\theta}\approx-
J_{\mathrm{ss}}Q\Delta\mathbf{u}_{\mathrm{ss}}.$ (38)

The above equation relates the update $\Delta\boldsymbol{\theta}$ to the
change in steady-state control input $Q\Delta\mathbf{u}_{\mathrm{ss}}$. In
order to maximally reduce $\mathcal{H}$ and hence $\lVert
Q\mathbf{u}_{\mathrm{ss}}\rVert^{2}_{2}$, the magnitude of the LHS should be
as big as possible. As a result, we want $\Delta\boldsymbol{\theta}$ to lie
fully inside the rowspace of $J_{\mathrm{ss}}R_{\mathrm{ss}}^{T}$, because
otherwise, $\Delta\boldsymbol{\theta}$ will partially be lost in its
nullspace, reducing the magnitude of the LHS. As
$\Delta\boldsymbol{\theta}=R_{\mathrm{ss}}Q\mathbf{u}_{\mathrm{ss}}$, we have
that the columnspace of $R_{\mathrm{ss}}Q$ needs to lie inside the rowspace of
$J_{\mathrm{ss}}R_{\mathrm{ss}}^{T}$. This condition is equivalent to the
columnspace of $Q$ lying inside the rowspace of $J_{\mathrm{ss}}$. Note that
the above arguments should be interpreted on an intuitive level, a detailed
mathematical consideration of the above arguments is more nuanced, but leads
to the same high-level intuition.

###  A.3 Effect of the nonlinearity $\phi$ on the weight updates

The main effect of including the nonlinearity $\phi$ into the Strong-DFC
weight updates (Eq. 4), in contrast to the steady-state updates (Eq. 7) of
Theorem 1, is ensuring that saturated neurons refrain from updating their
synaptic weights for the current data sample, as
$\phi(\mathbf{v}_{i})-\phi(\mathbf{v}^{\mathrm{ff}}_{i})\approx 0$ in this
case.

Meulemans et al. (2021) showed that including $\phi$ in the learning rule is a
useful heuristic for improving the alignment of the DFC weight updates with
the theoretically ideal update (which in our case would be the negative
gradient of $\mathcal{H}$), when the column space condition on $Q$ is not
perfectly satisfied. Consequently, including $\phi$ also improves the
performance in this case. When the column space condition on $Q$ is perfectly
satisfied, including $\phi$ in the learning rule makes no difference.

###  A.4 Steady-state weight updates vs. continuous weight updates

Our theoretical results assume that the weight update is performed at steady
state, giving rise to the update defined in Eq. (7). However, in practice,
plasticity is always on, giving rise to the continuous plasticity dynamics
defined in Eq. (4). Following Meulemans et al. (2021), we show that when the
steady state is quickly reached, the steady-state update (7) approximates the
resulting continuous update (4) well.

Assuming a separation of time scales between the network dynamics and the
synaptic plasticity ($\tau_{v},\tau_{u}\ll\tau_{W}$), the accumulated weight
update resulting from the plasticity dynamics (4) is given by

$\displaystyle\int_{0}^{t_{1}}\text{d}W_{i}=\int_{0}^{t_{\mathrm{ss}}}\text{d}W_{i}+\frac{t_{1}-t_{\mathrm{ss}}}{\tau_{W}}\left(\phi(\mathbf{v}_{i,\mathrm{ss}})-\phi(\mathbf{v}^{\mathrm{ff}}_{i,\mathrm{ss}})\right)\mathbf{r}_{i-1,\mathrm{ss}}^{T}.$
(39)

Now, if the dynamics quickly settle towards their steady state, relative to
the time interval the input is presented, i.e., $t_{\mathrm{ss}}\ll t_{1}$,
the weight update is dominated by the second term in the RHS of the above
equation, which corresponds to the steady-state weight update with $\phi$
included. Furthermore, the controller quickly pushes the activity close to its
steady state, with possibly some oscillations around it, so the first term
$\int_{0}^{t_{\mathrm{ss}}}\text{d}W_{i}$ will approximate the steady state
well.

In conclusion, the Strong-DFC weight update resulting from the plasticity
dynamics described in Eq. (4) will approximate the steady-state weight update
$\left(\phi(\mathbf{v}_{i,\mathrm{ss}})-\phi(\mathbf{v}^{\mathrm{ff}}_{i,\mathrm{ss}})\right)\mathbf{r}_{i-1,\mathrm{ss}}^{T}$.
Combined with the arguments of Appendix A.3, we see that the Strong-DFC weight
update approximates the steady-state update of Theorem 1 well.

###  A.5 The minimizing control framework is compatible with various
controller types

In Section 3, we introduced our theory on the minimizing control framework
tailored towards a proportional intergral (PI) controller used by DFC.
However, our theory can easily be extended to other controller types that
reach the correct steady state.

In the proof of Theorem 1, the only assumption on the controller is its steady
state equation (21), given by Eq. (21):

$\displaystyle
F(\mathbf{u}_{\mathrm{ss}},\boldsymbol{\theta})\triangleq\mathbf{e}_{\mathrm{ss}}(\mathbf{u}_{\mathrm{ss}},\boldsymbol{\theta})-\tilde{\alpha}\mathbf{u}_{\mathrm{ss}}=0.$

Therefore, Theorem 1 applies to all stable controllers that give rise to this
same steady-state equation, $F(\mathbf{u}_{\mathrm{ss}},\boldsymbol{\theta})$.

For example, pure proportional control, $\mathbf{u}(t)=k\mathbf{e}(t)$ gives
rise to the same steady-state equation
$F(\mathbf{u}_{\mathrm{ss}},\boldsymbol{\theta})$, with
$\tilde{\alpha}=\frac{1}{k}$. For pure integral control, we have the same
steady-state equation $F(\mathbf{u}_{\mathrm{ss}},\boldsymbol{\theta})$ but
with $\tilde{\alpha}=\alpha$. If we add derivative control to any of the above
controller types, we reach the same steady-state equation as the derivative of
the error vanishes at steady state.

In conclusion, Theorem 1 is compatible with all possible combinations of the
proportional-integral-derivative (PID) controller framework, except for the
pure derivative control. Furthermore, all other controller types that satisfy
the steady-state equation $F(\mathbf{u}_{\mathrm{ss}},\boldsymbol{\theta})$
are also compatible with Theorem 1.

###  A.6 Adapting the minimizing control framework for classification with
softmax and cross-entropy loss

For classification with neural networks, one uses conventionally a one-hot
label $\mathbf{r}_{L}^{\mathrm{true}}$ in combination with a cross-entropy
loss and a softmax layer as output of the network. In our minimizing control
framework, the controller drives the network to perfectly match the output
target $\mathbf{r}_{L}^{\mathrm{true}}$. However, a softmax layer only results
in a one-hot vector if it has an infinite input. To overcome this challenge,
we use soft targets which have $a$ on its label entry and
$\frac{1-a}{n_{L}-1}$ on the other entries, with $a$ close to 1 (e.g., 0.99).

A nice property of the combination of the cross-entropy loss with the softmax
output is that the curvature of the exponentials in the softmax cancel out the
curvature of the logarithms in the cross-entropy loss, leading to clean
gradients that do not saturate easily. However, if our controller uses as
control error $\mathbf{e}(t)=\mathbf{r}_{L}^{*}-\mathbf{r}_{L}(t)$, with
$\mathbf{r}_{L}^{*}$ the soft target and $\mathbf{r}_{L}(t)$ the softmax
output, it cannot make use of the cross-entropy loss to provide clean
gradients. To overcome this challenge, we use a linear output layer,
$\mathbf{r}_{L}$, and absorb the softmax into the cross-entropy loss, leading
to the combined loss:

$\displaystyle\mathcal{L}^{\mathrm{combined}}=-\sum_{b=1}^{B}\mathbf{r}_{L}^{\mathrm{true},(b)T}\log\big{(}\mathrm{Softmax}(\mathbf{r}_{L}^{(b)})\big{)}.$
(40)

Next, we slightly adapt the minimizing control framework to include this
combined loss in a principled way. Currently, the controller dynamics is given
by

$\displaystyle\mathbf{u}(t)$
$\displaystyle=k\mathbf{e}(t)+\mathbf{u}^{\mathrm{int}}(t),\quad\tau_{u}\frac{\text{d}}{\text{d}t}\mathbf{u}^{\mathrm{int}}(t)=\mathbf{e}(t)-\alpha\mathbf{u}^{\mathrm{int}}(t)$
(41) $\displaystyle\mathbf{e}(t)$
$\displaystyle=\mathbf{r}_{L}^{*}-\mathbf{r}_{L}(t).$ (42)

Note that $\mathbf{e}(t)$ can be interpreted as
$-\frac{\partial\mathcal{L}}{\partial\mathbf{r}_{L}}\rvert_{\mathbf{r}_{L}=\mathbf{r}_{L}(t)}$,
with
$\mathcal{L}=\frac{1}{2}\lVert\mathbf{r}_{L}^{*}-\mathbf{r}_{L}(t)\rVert_{2}^{2}$.
Now, if we use another loss $\mathcal{L}$ (i.e.,
$\mathcal{L}^{\mathrm{combined}}$), we can generalize the control error
$\mathbf{e}(t)$ to

$\displaystyle\mathbf{e}(t)=-\left.\frac{\partial\mathcal{L}}{\partial\mathbf{r}_{L}}\right\rvert_{\mathbf{r}_{L}=\mathbf{r}_{L}(t)}^{T}.$
(43)

For the cross entropy loss with the softmax, this boils down to

$\displaystyle\mathbf{e}(t)=-\left.\frac{\partial\mathcal{L}}{\partial\mathbf{r}_{L}}\right\rvert_{\mathbf{r}_{L}=\mathbf{r}_{L}(t)}^{T}=\mathbf{p}^{*}-\mathrm{Softmax}(\mathbf{r}_{L}(t)),$
(44)

with $\mathbf{p}^{*}$ the soft target.

##  Appendix B Stability analysis

In this section, we build further upon the stability analysis results of
Meulemans et al. (2021), which are applicable to Strong-DFC with minimal
adaptations.

As a first step, we investigate a pure integral controller ($k=0$) and assume
a separation of timescales $\tau_{v}\ll\tau_{u}$, i.e., we can replace the
network dynamics by a deterministic function and only retain the controller
dynamics. This provides us with the following controller dynamics:

$\displaystyle\tau_{u}\dot{\mathbf{u}}(t)=\mathbf{e}\left(\mathbf{u}(t)\right)-\alpha\mathbf{u}(t),$
(45)

where the error, $\mathbf{e}$, can be seen as a deterministic function of the
control input, $\mathbf{u}(t)$ (for simplicity, we omit the dependence on the
network parameters and input, as they are considered fixed during the
controller dynamics). To investigate the local stability at steady state, we
consider the eigenvalues of the Jacobian matrix of these dynamics evaluated at
steady state. Using Eq. (23), we have that:

$\displaystyle\frac{\partial}{\partial\mathbf{u}}\left.\left(\mathbf{e}\left(\mathbf{u}\right)-\alpha\mathbf{u}\right)\right|_{\mathbf{u}=\mathbf{u}_{\mathrm{ss}}}=-\left(J_{\mathrm{ss}}Q+\alpha
I\right).$ (46)

Hence, the condition for local stability around the steady state is that all
eigenvalues of $-\left(J_{\mathrm{ss}}Q+\alpha I\right)$ should be negative,
resulting in Condition 2.

When considering the proportional integral controller without a separation of
timescales, the local stability condition becomes harder to interpret.
Following a similar derivation as proposed by Meulemans et al. (2021), we
reach the following Jacobian of the network-controller dynamics:

$\displaystyle
A_{\mathrm{PI}}=-\begin{bmatrix}\frac{1}{\tau_{v}}(I-W\phi^{\prime})&-\frac{1}{\tau_{v}}Q\\\
\frac{1}{\tilde{\tau}_{u}}S\phi^{\prime}-\frac{k}{\tau_{v}}S\phi^{\prime}(I-W\phi^{\prime})&\frac{k}{\tau_{v}}S\phi^{\prime}Q\frac{\tilde{\alpha}}{\tilde{\tau}_{u}}I,\end{bmatrix}$
(47)

with $\tilde{\tau}_{u}\triangleq\frac{\tau_{u}}{1+\alpha k}$,
$\tilde{\alpha}\triangleq\frac{\alpha}{1+\alpha k}$,
$\phi^{\prime}\triangleq\left.\frac{\partial\phi(\mathbf{v})}{\partial\mathbf{v}}\right|_{\mathbf{v}=\mathbf{v}_{\mathrm{ss}}}$,
and

$\displaystyle W$ $\displaystyle\triangleq\begin{bmatrix}0&0&0&\ldots&0\\\
W_{2}&0&0&\ldots&0\\\ 0&W_{3}&0&\ldots&0\\\ 0&0&\ddots&\ldots&0\\\
0&0&\ldots&W_{L}&0\end{bmatrix}$ (48) $\displaystyle S$
$\displaystyle\triangleq\begin{bmatrix}0&\ldots&0&I\end{bmatrix}$ (49)

The condition for local stability is now that all eigenvalues of
$A_{\mathrm{PI}}$ are negative. However, this condition provides no
straightforward interpretations for defining conditions on the feedback
weights, $Q$, such that local stability is achieved. One can apply
Gershgoring’s circle theorem for finding sufficient conditions on
$J_{\mathrm{ss}}$ and Q to ensure local stability, however, the resulting
conditions are too conservative and hence not usable for designing feedback
weight learning rules.

##  Appendix C Proofs and discussion for Section 4

###  C.1 Stochastic dynamics for the single-phase learning setting

In this section, we provide a detailed description of the stochastic dynamics
used for training the forward and feedback weights in a single phase. The
network and controller obey the following dynamics:

$\displaystyle\tau_{v}\frac{\text{d}}{\text{d}t}\mathbf{v}_{i}(t)$
$\displaystyle=-\mathbf{v}_{i}(t)+W_{i}\phi\big{(}\mathbf{v}_{i-1}(t)\big{)}+Q_{i}\mathbf{u}(t)+\sigma\boldsymbol{\epsilon}_{i}(t)$
(50)
$\displaystyle\tau_{u}\frac{\text{d}}{\text{d}t}\mathbf{u}^{\mathrm{int}}(t)$
$\displaystyle=\mathbf{r}_{L}^{*}-\mathbf{r}_{L}(t)-\alpha\mathbf{u}^{\mathrm{int}}(t)$
(51) $\displaystyle\mathbf{u}(t)$
$\displaystyle=k\big{(}\mathbf{r}_{L}^{*}-\mathbf{r}_{L}(t)\big{)}+\mathbf{u}^{\mathrm{int}}(t)$
(52)

with $\sigma^{2}$ the noise variance and $\boldsymbol{\epsilon}(t)$
exponentially filtered white noise, with time constant $\tau_{\epsilon}$:

$\displaystyle\boldsymbol{\epsilon}_{i}(t)=\frac{1}{\tau_{\epsilon}}\int_{-\infty}^{t}\exp\big{(}-\frac{1}{\tau_{\epsilon}}(t-t^{\prime})\big{)}\boldsymbol{\xi}_{i}(t^{\prime})\text{d}t^{\prime},$
(54)

where $\boldsymbol{\xi}$ is white noise. The time constant $\tau_{\epsilon}$
determines the auto-correlation of $\boldsymbol{\epsilon}$. For simplicity, we
take that the noise $\boldsymbol{\epsilon}$ enters solely in the feedback
compartment:

$\displaystyle\mathbf{v}^{\mathrm{fb}}_{i}(t)=Q_{i}\mathbf{u}(t)+\sigma\boldsymbol{\epsilon}_{i}(t).$
(55)

However, all the arguments and theorems of the following sections also apply
for noise entering all compartments.

Now, we define the following error-based learning rule for the forward weights
and anti-Hebbian learning rule for the feedback weights:

$\displaystyle\tau_{W}\frac{\text{d}}{\text{d}t}W_{i}(t)$
$\displaystyle=\Big{(}\phi\big{(}\mathbf{v}_{i}(t)\big{)}-\phi\big{(}W_{i}\mathbf{r}_{i-1}(t)\big{)}\Big{)}\bar{\mathbf{r}}_{i-1}^{T}$
(56) $\displaystyle\tau_{Q}\frac{\text{d}}{\text{d}t}Q_{i}(t)$
$\displaystyle=-\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{L-i}\mathbf{v}^{\mathrm{fb}}_{i}(t)\tilde{\mathbf{u}}(t)^{T}-\beta
Q_{i}$ (57) $\displaystyle\tilde{\mathbf{u}}(t)$
$\displaystyle=\mathbf{u}(t)-\bar{\mathbf{u}}(t)$ (58)
$\displaystyle\bar{\mathbf{u}}(t)$
$\displaystyle=\frac{1}{\tau_{f}}\int_{0}^{t}\exp\big{(}\frac{1}{\tau_{f}}(t-t^{\prime})\big{)}\mathbf{u}(t^{\prime})\text{d}t^{\prime}$
(59) $\displaystyle\bar{\mathbf{r}}_{i}(t)$
$\displaystyle=\frac{1}{\tau_{f}}\int_{0}^{t}\exp\big{(}\frac{1}{\tau_{f}}(t-t^{\prime})\big{)}\mathbf{r}_{i}(t^{\prime})\text{d}t^{\prime}$
(60)

with $\tau_{f}$ the filtering time constant used to high-pass filter
$\mathbf{u}$ and low-pass filter $\mathbf{r}_{i}$.

###  C.2 Intuition behind the feedback weight learning rule

Here we expand on the intuition behind our feedback weight learning rule. The
central technique, inspired by the weight mirroring method (Akrout et al.,
2019) and DFC (Meulemans et al., 2021), is to inject white noise in the
network that carries information about the network Jacobian,
$J_{\mathrm{ss}}$, towards the output. If we now measure the correlation of
the output noise w.r.t. the injected noise in the layers, we can recover the
information on $J_{\mathrm{ss}}$. Below, we summarize the intuitive
explanation on a normal feedforward neural network, provided by Meulemans et
al. (2021). We end with explaining how to go from this intuitive explanation
to our single-phase feedback weight plasticity, and how it differs from the
original feedback weight plasticity of DFC.

Now, let us consider a normal feedforward neural network
$\mathbf{r}_{i}=\phi(\mathbf{v}_{i})=\phi(W_{i}\mathbf{r}_{i-1})$. We perturb
each layer with white noise $\sigma\boldsymbol{\xi}$, which is propagated
forward through the network:

$\displaystyle\tilde{\mathbf{v}}_{i}$
$\displaystyle=W_{i}\phi(\tilde{\mathbf{v}}_{i-1})+\sigma\boldsymbol{\xi}_{i}$
(61) $\displaystyle\tilde{\mathbf{r}}_{0}$ $\displaystyle=\mathbf{r}_{0}$ (62)

Assuming $\sigma$ is small, we can approximate the network output as:

$\displaystyle\tilde{\mathbf{r}}_{L}\approx\mathbf{r}_{L}+\sigma
J\boldsymbol{\xi},$ (63)

with
$J=\begin{bmatrix}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{1}},&\ldots&,\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{L}}\end{bmatrix}$
and $\boldsymbol{\xi}$ the concatenation of all $\boldsymbol{\xi}_{i}$. Now,
we can use an output error $\mathbf{e}=\mathbf{r}_{L}-\tilde{\mathbf{r}}_{L}$
to extract the noise, and use this in a feedback learning rule $\Delta
Q=-\boldsymbol{\xi}\mathbf{e}^{T}-\beta Q$. In expectation, this update rule
results in

$\displaystyle\mathbb{E}[\Delta Q]=\sigma^{2}J^{T}-\beta Q.$ (64)

Hence, we see that this simple update rule drives the feedback weights to
align with $J^{T}$.

Although the above explanation grasps the intuition behind the feedback
learning rule (9), there are some important differences when we consider the
continuous single-phase plasticity dynamics.

  1. 1.

The pre-synaptic plasticity signal of the feedback weight update (9) is the
high-pass filtered control signal $\mathbf{u}_{\mathrm{hp}}$ instead of the
output noise we used before. When training the feedback and forward weights
simultaneously, $\mathbf{u}$ consists of two parts: (i) the control signal
used to drive the network in expectation towards the output target
$\mathbf{r}_{L}^{*}$, and (ii), the output noise fluctuations integrated over
time by the controller. As we are interested in the noise fluctuations, we
need to extract them, which we do by using the high-pass filtered version of
$\mathbf{u}$.

  2. 2.

The postsynaptic plasticity signal is the activity in the feedback compartment
$\mathbf{v}^{\mathrm{fb}}_{i}=Q_{i}\mathbf{u}+\boldsymbol{\epsilon}_{i}$,
instead of purely the noise. As $Q_{i}\mathbf{u}$ will correlate with
$\mathbf{u}_{\mathrm{hp}}$, and $\mathbf{u}_{\mathrm{hp}}$ contains the
integrated noise fluctuations over time, instead of the direct noise
fluctuations, the feedback weights $Q$ will not align exactly with $J^{T}$,
but with $J^{T}M$ instead, where $M$ is a positive definite matrix.

  3. 3.

To incorporate noise in the continuous dynamics of Strong-DFC, we need to use
stochastic differential equations (SDEs), instead of discrete noise
perturbations. Unlike discrete feedforward neural networks, the network
dynamics do not propagate the layer activity instantaneously towards the
output, but have a delay proportional to $\tau_{v}$. As white noise has no
correlation in time, naively correlating the output noise fluctuations with
the current layer noise would result in a signal of zero due to the time
delays. In order to have correlations over time, we use exponentially filtered
white noise $\boldsymbol{\epsilon}$ with time constant $\tau_{\epsilon}$,
instead of pure white noise $\boldsymbol{\xi}$. As the time delay between the
output noise fluctuations and the layer noise increases for upstream layers,
we need to compensate for the decreasing correlation signal by scaling the
plasticity updates with
$\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{L-i}$, as is done in Eq.
(9).

Although our single-phase plasticity dynamics (9) are closely related to the
two-phase plasticity dynamics of DFC introduced by Meulemans et al. (2021),
and these authors discussed briefly the possibility of having single-phase
plasticity dynamics, our plasticity dynamics introduce some important
advances, both in theory and practice.

  1. 1.

Our minimizing control framework explicitly motivates learning the forward and
feedback weights simultaneously, as the feedback weights are trained to
satisfy Condition 1, which is evaluated at the steady state of the network for
training the forward weights. By contrast, the theory for DFC in the weak-
nudging regime requires the feedback weights to be trained at the network
steady state without influence from the controller, and therefore different
from the network steady state used for training the forward weights (Meulemans
et al., 2021).

  2. 2.

In their theoretical results on the feedback learning rule, Meulemans et al.
(2021) require that the controller leakage $\alpha$ is big. However, this
cannot be combined with learning the feedback and feedforward weights
simultaneously, as for training the feedforward weights, we need
$\alpha\rightarrow 0$. In our new theoretical results, we overcome this
important challenge.

  3. 3.

In their theoretical results on the feedback learning rule, Meulemans et al.
(2021) assume instant network dynamics ($\tau_{v}\rightarrow 0$). However, in
practice, a finite network time constant is used, giving rise to delays in the
propagation of noise fluctuations to the output, leading to weaker correlation
signals for upstream layers. Our new theory embraces finite time constants and
investigates in detail the required relation between the various time
constants $\tau_{v}$, $\tau_{\epsilon}$, $\tau_{u}$, and $\tau_{f}$. As a
result, we introduce the scaling term
$\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{L-i}$ which improves the
feedback learning significantly.

###  C.3 Proof of Theorem 3

In the following, we prove Theorem 3, which we state in its full form below.

######  Theorem S4.

Assume a separation of timescales $\tau_{v},\tau_{\epsilon}\ll\
\tau_{u}\ll\tau_{f}\ll\tau_{Q}$, integral control ($k=0$), stable network
dynamics and small noise perturbations $\sigma\ll\lVert\mathbf{v}_{i}\rVert$.
Then, for a fixed data sample, the feedback plasticity dynamics (9) let the
first moment of $Q$ converge approximately towards the following configuration
that satisfies Conditions 1 and 2:

$\displaystyle\mathbb{E}[Q_{\mathrm{ss}}]\approx J_{\mathrm{ss}}^{T}M,$ (65)

with $M$ a positive definite symmetric matrix.

###### Proof.

We consider the noisy dynamics defined in Eq. (50-52), with $k=0$. Let us
define $\bar{\mathbf{v}}_{\mathrm{ss}}$ and $\bar{\mathbf{u}}_{\mathrm{ss}}$
as the steady state of the mean trajectory of the noisy dynamics, and
$\tilde{\mathbf{v}}\triangleq\mathbf{v}-\bar{\mathbf{v}}_{\mathrm{ss}}$ and
$\tilde{\mathbf{u}}\triangleq\mathbf{u}-\bar{\mathbf{u}}_{\mathrm{ss}}$ as the
(noisy) perturbations around the mean trajectories. As
$\sigma\ll\lVert\mathbf{v}_{i}\rVert$ and the dynamics are stable (hence,
contracting), we have that
$\lVert\tilde{\mathbf{v}}\rVert\ll\lVert\mathbf{v}\rVert$ and
$\lVert\tilde{\mathbf{u}}\rVert\ll\lVert\mathbf{u}\rVert$, allowing us to do a
first-order Taylor expansion of the dynamics, which for simplicity we assume
to be exact:

$\displaystyle\tau_{v}\frac{\text{d}}{\text{d}t}\tilde{\mathbf{v}}(t)$
$\displaystyle=-\bar{\mathbf{v}}_{\mathrm{ss}}+W\phi\big{(}\bar{\mathbf{v}}_{\mathrm{ss}}\big{)}+U\mathbf{r}_{0}+Q\bar{\mathbf{u}}_{\mathrm{ss}}(t)-\tilde{\mathbf{v}}(t)+W\phi^{\prime}\tilde{\mathbf{v}}(t)+Q\tilde{\mathbf{u}}(t)+\sigma\boldsymbol{\epsilon}(t)$
(66)
$\displaystyle=-\tilde{\mathbf{v}}(t)+W\phi^{\prime}\tilde{\mathbf{v}}(t)+Q\tilde{\mathbf{u}}(t)+\sigma\boldsymbol{\epsilon}(t)$
(67) $\displaystyle\tau_{u}\frac{\text{d}}{\text{d}t}\tilde{\mathbf{u}}(t)$
$\displaystyle=\mathbf{r}_{L}^{*}-S\phi(\bar{\mathbf{v}}_{\mathrm{ss}})-\alpha\bar{\mathbf{u}}_{\mathrm{ss}}-S\phi^{\prime}\tilde{\mathbf{v}}(t)-\alpha\tilde{\mathbf{u}}(t)$
(68)
$\displaystyle=-S\phi^{\prime}\tilde{\mathbf{v}}(t)-\alpha\tilde{\mathbf{u}}(t)$
(69)

where we define
$\phi^{\prime}\triangleq\left.\frac{\partial\phi(\mathbf{v})}{\partial\mathbf{v}}\right|_{\mathbf{v}=\bar{\mathbf{v}}_{\mathrm{ss}}}$,

$\displaystyle W$ $\displaystyle\triangleq\begin{bmatrix}0&0&0&\ldots&0\\\
W_{2}&0&0&\ldots&0\\\ 0&W_{3}&0&\ldots&0\\\ 0&0&\ddots&\ldots&0\\\
0&0&\ldots&W_{L}&0\end{bmatrix}$ (70) $\displaystyle U$
$\displaystyle\triangleq\begin{bmatrix}W_{1}^{T}&0&\ldots&0\end{bmatrix}^{T}$
(71) $\displaystyle S$
$\displaystyle\triangleq\begin{bmatrix}0&\ldots&0&I\end{bmatrix}$ (72)

and use that
$-\bar{\mathbf{v}}_{\mathrm{ss}}+W\phi\big{(}\bar{\mathbf{v}}_{\mathrm{ss}}\big{)}+U\mathbf{r}_{0}+Q\bar{\mathbf{u}}_{\mathrm{ss}}(t)=0$
and
$\mathbf{r}_{L}^{*}-S\phi(\bar{\mathbf{v}}_{\mathrm{ss}})-\alpha\bar{\mathbf{u}}_{\mathrm{ss}}=0$
at the steady state of the mean trajectory. Similar to $\Delta\mathbf{v}_{i}$,
we define
$\Delta\tilde{\mathbf{v}}_{i}\triangleq\tilde{\mathbf{v}}_{i}-W_{i}\phi^{\prime}_{i-1}\tilde{\mathbf{v}}_{i-1}$
with
$\phi^{\prime}_{i}\triangleq\left.\frac{\partial\phi(\mathbf{v}_{i})}{\partial\mathbf{v}_{i}}\right|_{\mathbf{v}=\bar{\mathbf{v}}_{\mathrm{ss},i}}$.
As $\Delta\tilde{\mathbf{v}}_{1}=\tilde{\mathbf{v}}_{1}$ and
$\tilde{\mathbf{v}}_{i}=\Delta\tilde{\mathbf{v}}_{i}+W_{i}\phi^{\prime}_{i-1}$,
we have that

$\displaystyle
S\phi^{\prime}\tilde{\mathbf{v}}=\phi^{\prime}_{L}\tilde{\mathbf{v}}_{L}=J_{\mathrm{ss}}\Delta\tilde{\mathbf{v}},$
(73)

with
$J_{\mathrm{ss}}\triangleq[\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{1}},...,\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{L}}]$
evaluated at $\bar{\mathbf{v}}_{\mathrm{ss}}$. Using
$\frac{\text{d}}{\text{d}t}\Delta\tilde{\mathbf{v}}=(I-W\phi^{\prime})\frac{\text{d}}{\text{d}t}\tilde{\mathbf{v}}$
and
$\boldsymbol{\gamma}(t)\triangleq[\Delta\tilde{\mathbf{v}}(t)^{T},\tilde{\mathbf{u}}(t)^{T}]^{T}$
to simplify notation we get

$\displaystyle\frac{\text{d}}{\text{d}t}\boldsymbol{\gamma}(t)=-\underbrace{\begin{bmatrix}\frac{1}{\tau_{v}}(I-W\phi^{\prime})&-\frac{1}{\tau_{v}}(I-W\phi^{\prime})Q\\\
\frac{1}{\tau_{u}}J_{\mathrm{ss}}&\frac{\alpha}{\tau_{u}}I\end{bmatrix}}_{\triangleq
A}\boldsymbol{\gamma}(t)+\underbrace{\begin{bmatrix}\frac{\sigma}{\tau_{v}}(I-W\phi^{\prime})\\\
0\end{bmatrix}}_{\triangleq B}\boldsymbol{\epsilon}(t)$ (74)

We solve this set of linear time-invariant differential equations using the
variation of constants method (Särkkä & Solin, 2019), taking
$\tilde{\mathbf{v}}(t_{0})=0$ and $\tilde{\mathbf{u}}(t_{0})=0$ (i.e., we
consider that the network and controller dynamics are already converged at
$t_{0}$).

$\displaystyle\boldsymbol{\gamma}(t)=\int_{t_{0}}^{t}\exp\big{(}-A(t-t^{\prime})\big{)}B\boldsymbol{\epsilon}(t^{\prime})\text{d}t^{\prime}.$
(75)

Now that we solved the differential equations for $\tilde{\mathbf{v}}$ and
$\tilde{\mathbf{u}}$, we turn our attention to the feedback weight plasticity
rule (9). We start by investigating
$\mathbb{E}[\mathbf{v}^{\mathrm{fb}}(t)\tilde{\mathbf{u}}(t)^{T}]$. As the
noise $\boldsymbol{\epsilon}$ enters through the apical compartment, we have
that:

$\displaystyle\mathbf{v}^{\mathrm{fb}}(t)$
$\displaystyle=Q\mathbf{u}(t)+\sigma\boldsymbol{\epsilon}(t)=Q\bar{\mathbf{u}}_{\mathrm{ss}}+Q\tilde{\mathbf{u}}(t)+\sigma\boldsymbol{\epsilon}(t)$
(76)
$\displaystyle\mathbb{E}[\mathbf{v}^{\mathrm{fb}}(t)\tilde{\mathbf{u}}(t)^{T}]$
$\displaystyle=Q\bar{\mathbf{u}}_{\mathrm{ss}}\mathbb{E}[\tilde{\mathbf{u}}(t)^{T}]+Q\mathbb{E}[\tilde{\mathbf{u}}(t)\tilde{\mathbf{u}}(t)^{T}]+\sigma\mathbb{E}[\boldsymbol{\epsilon}(t)\tilde{\mathbf{u}}(t)^{T}]$
(77)

The first term of Eq. (77) is equal to zero, as $\boldsymbol{\epsilon}$ and
consequently $\tilde{\mathbf{u}}$ are zero-mean. We now turn our attention to
the third term. Using $S_{u}\triangleq\begin{bmatrix}0&I\end{bmatrix}$ and
$\tilde{\mathbf{u}}=S_{u}\boldsymbol{\gamma}$, we have that:

$\displaystyle\mathbb{E}[\tilde{\mathbf{u}}(t)\boldsymbol{\epsilon}(t)^{T}]$
$\displaystyle=S_{u}\mathbb{E}[\boldsymbol{\gamma}(t)\boldsymbol{\epsilon}(t)^{T}]=S_{u}\int_{t_{0}}^{t}\exp\big{(}{-}A(t-t^{\prime})\big{)}B\mathbb{E}[\boldsymbol{\epsilon}(t^{\prime})\boldsymbol{\epsilon}(t)^{T}]\text{d}t^{\prime}$
(78)
$\displaystyle=\frac{1}{2\tau_{\epsilon}}S_{u}\int_{t_{0}}^{t}\exp\big{(}{-}(A+\frac{1}{\tau_{\epsilon}}I)(t-t^{\prime})\big{)}B\text{d}t^{\prime},$
(79)

for which we used
$\mathbb{E}[\boldsymbol{\epsilon}(t^{\prime})\boldsymbol{\epsilon}(t)^{T}]=\frac{1}{2\tau_{\epsilon}}\exp\big{(}\frac{1}{\tau_{\epsilon}}|t-t^{\prime}|\big{)}$.
Using $\int_{0}^{T}\exp(At)\text{d}t=\big{(}\exp(AT)-I\big{)}A^{-1}$ and
assuming stability and $t-t_{0}\gg\tau_{u}$, which is justified by the
condition $\tau_{Q}\gg\tau_{u}$, we can solve the previous integral:

$\displaystyle\mathbb{E}[\tilde{\mathbf{u}}(t)\boldsymbol{\epsilon}(t)^{T}]=\frac{1}{2\tau_{\epsilon}}S_{u}\big{(}A+\frac{1}{\tau_{\epsilon}}I\big{)}^{-1}B.$
(80)

Using the solution for the inverse of $2\times 2$ block matrices of Lu & Shiou
(2002), we get

$\displaystyle\mathbb{E}[\tilde{\mathbf{u}}(t)\boldsymbol{\epsilon}(t)^{T}]=$
(81)
$\displaystyle{-}\frac{1}{2\tau_{\epsilon}}\left(\frac{1}{\tau_{\epsilon}}+\frac{\alpha}{\tau_{u}}\right)^{-1}\frac{1}{\tau_{u}}J_{\mathrm{ss}}\left(\frac{1}{\tau_{\epsilon}}I+\frac{1}{\tau_{v}}\left(I-W\phi^{\prime}\right)+\left(\frac{1}{\tau_{\epsilon}}+\frac{\alpha}{\tau_{u}}\right)^{-1}\frac{1}{\tau_{v}}\left(I-W\phi^{\prime}\right)Q\frac{1}{\tau_{u}}J_{\mathrm{ss}}\right)^{-1}\frac{\sigma}{\tau_{v}}\left(I-W\phi^{\prime}\right)$
(82)

Now, using the condition $\tau_{v},\tau_{\epsilon}\ll\tau_{u}$, we can
approximate the above with:

$\displaystyle\mathbb{E}[\tilde{\mathbf{u}}(t)\boldsymbol{\epsilon}(t)^{T}]$
$\displaystyle\approx{-}\frac{\sigma}{2\tau_{u}}J_{\mathrm{ss}}\left(\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)I-W\phi^{\prime}\right)^{-1}\left(I-W\phi^{\prime}\right)$
(83)

Using Lemma S5, we have that:

$\displaystyle\mathbb{E}[\tilde{\mathbf{u}}(t)\boldsymbol{\epsilon}(t)^{T}]$
$\displaystyle\approx{-}\frac{\sigma}{2\tau_{u}}\left.\begin{bmatrix}\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{-L}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{1}}&\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{-(L-1)}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{2}}&\ldots&\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{-1}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{L}}\end{bmatrix}\right|_{\mathbf{v}=\mathbf{v}_{\mathrm{ss}}}$
(84)

Now we turn our attention towards the second term of Eq. (77). As
$\mathbb{E}[\tilde{\mathbf{u}}(t)\tilde{\mathbf{u}}(t)^{T}]$ is a covariance
matrix, it is positive semi-definite.

Combining Eq. (9) with Eq. (81) and the above insights, we get

$\displaystyle\tau_{Q}\frac{\text{d}}{\text{d}t}\mathbb{E}[Q(t)]\triangleq\tau_{Q}\frac{\text{d}}{\text{d}t}\bar{Q}(t)=\frac{\sigma^{2}}{2\tau_{u}}J_{\mathrm{ss}}^{T}-\bar{Q}(\tilde{M}+\beta
I),$ (85)

with $\tilde{M}$ a positive semi-definite matrix. We used the fact that the
layer-wise scaling with
$\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{L-i}$ in Eq. (9) cancels out
the scaling factors $\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{-(L-i)}$
in Eq. (81). Using $M\triangleq\tilde{M}+\beta I$, the stable point of the
above dynamics is

$\displaystyle\bar{Q}_{\mathrm{ss}}=\frac{\sigma^{2}}{\tau_{u}}J_{\mathrm{ss}}^{T}M^{-1}.$
(86)

As $M$ is positive definite, the dynamics are stable and the above steady
state is reached. Furthermore, as $J_{\mathrm{ss}}^{T}M^{-}1$ is a permutation
of the columns of $J_{\mathrm{ss}}^{T}$, it satisfies Condition 1; and as
$M^{-1}$ is positive definite, it satisfies Condition 2.

Now the only thing left to show is that
$\mathbf{u}_{\mathrm{hp}}\triangleq\mathbf{u}-\mathbf{u}_{\mathrm{lp}}\overset{?}{=}\tilde{\mathbf{u}}$,
with $\mathbf{u}_{\mathrm{lp}}$ the low-pass filtered version of $\mathbf{u}$:

$\displaystyle\mathbf{u}_{\mathrm{lp}}(t)=\frac{1}{\tau_{f}}\int_{t_{0}}^{t}\exp\left({-}\frac{1}{\tau_{f}}(t-t^{\prime})\right)\mathbf{u}(t^{\prime})\text{d}t^{\prime}.$
(87)

As we have that $\tau_{f}\gg\tau_{u},\tau_{\epsilon},\tau_{v}$, the low-pass
filtering will remove all fluctuations and we have that
$\mathbf{u}_{\mathrm{lp}}=\bar{\mathbf{u}}_{\mathrm{ss}}$. Consequently, we
have that
$\mathbf{u}_{\mathrm{hp}}=\mathbf{u}-\bar{\mathbf{u}}_{\mathrm{ss}}=\tilde{\mathbf{u}}$,
thereby concluding the proof. ∎

######  Lemma S5.

For some scalar $a\neq 0$, we have that

$\displaystyle
J\left(aI-W\phi^{\prime}\right)^{-1}\left(I-W\phi^{\prime}\right)=\left.\begin{bmatrix}a^{-L}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{1}}&a^{-(L-1)}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{2}}&\ldots&a^{-1}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{L}}\end{bmatrix}\right|_{\mathbf{v}=\mathbf{v}^{\prime}}$
(88)

with $J$ and $\phi^{\prime}$ evaluated at some activity pattern
$\mathbf{v}^{\prime}$.

###### Proof.

First, we note that $\left(aI-W\phi^{\prime}\right)$ is of full rank for
$a\neq 0$, as it is a triangular matrix with $a$ on its diagonal. The inverse
of $\left(aI-W\phi^{\prime}\right)$ is given by:

$\displaystyle\left(aI-W\phi^{\prime}\right)^{-1}=\begin{bmatrix}a^{-1}I&0&0&\ldots&0\\\
a^{-2}J_{2,1}&a^{-1}I&0&\ldots&0\\\
a^{-3}J_{2,1}&a^{-2}J_{3,2}&a^{-1}I&\ldots&0\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
a^{-L}J_{L,1}&\ldots&\ldots&a^{-2}J_{L,L-1}&a^{-1}I\end{bmatrix}$ (89)

with $J_{i,j}\triangleq\frac{\partial\mathbf{v}_{i}}{\partial\mathbf{v}_{j}}$.
One can easily verify this by computing
$\left(aI-W\phi^{\prime}\right)\left(aI-W\phi^{\prime}\right)^{-1}=I$ upon
noting that $J_{i,i-1}=W_{i}\phi^{\prime}_{i-1}$. Next, we have that:

$\displaystyle\left(aI-W\phi^{\prime}\right)^{-1}\left(I-W\phi^{\prime}\right)=\begin{bmatrix}a^{-1}I&0&0&\ldots&0\\\
(a^{-2}-a^{-1})J_{2,1}&a^{-1}I&0&\ldots&0\\\
(a^{-3}-a^{-2})J_{2,1}&(a^{-2}-a^{-1})J_{3,2}&a^{-1}I&\ldots&0\\\
\vdots&\ddots&\ddots&\ddots&\vdots\\\
(a^{-L}-a^{-(L-1)})J_{L,1}&\ldots&\ldots&(a^{-2}-a^{-1})J_{L,L-1}&a^{-1}I\end{bmatrix}$
(90)

Using
$J=[\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{1}},\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{2}},\ldots,\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{L}}]$,
we have

$\displaystyle
J\left(aI-W\phi^{\prime}\right)^{-1}\left(I-W\phi^{\prime}\right)=\begin{bmatrix}a^{-L}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{1}}&a^{-(L-1)}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{2}}&\ldots&a^{-1}\frac{\partial\mathbf{r}_{L}}{\partial\mathbf{v}_{L}},\end{bmatrix}$
(91)

thereby concluding the proof. ∎

###  C.4 Noise can bias the feedforward weight updates

As observed by Meulemans et al. (2021), when we inject noise into the dynamics
while training the forward weights, the resulting noise correlations can bias
the forward weight updates. This is a general issue for methods using error
feedback for learning while having realistic noise dynamics, as discussed in
more detail in App. C.4 in Meulemans et al. (2021). Here, we investigate how
these issues manifest themselves in Strong-DFC and how we can solve them.

The forward weight update is given by Eq. (4). If the noise fluctuations in
the postsynaptic term $\phi(\mathbf{v}_{i}(t))-\phi(W_{i}\mathbf{r}_{i-1}(t))$
are correlated to the noise fluctuations in the presynaptic term
$\mathbf{r}_{i-1}(t)$, the expected forward weight update contains this extra
correlation term and is hence biased.

The difference $\phi(\mathbf{v}_{i}(t))-\phi(W_{i}\mathbf{r}_{i-1}(t))$ is
mainly determined by the feedback input $Q\mathbf{u}(t)$. Using the notation
of Section C.3, we divide the control signal of the converged dynamics in two
parts: $\mathbf{u}(t)=\mathbf{u}_{\mathrm{ss}}+\tilde{\mathbf{u}}(t)$, with
$\mathbf{u}_{\mathrm{ss}}$ the part of the control signal useful for learning
the forward weights, and $\tilde{\mathbf{u}}$ the noise fluctuations that are
used for learning the feedback weights. These noise fluctuations are a
function of the noise fluctuations of all network layers (e.g., see Eq. 75).
Hence, the noise fluctuations of
$\phi(\mathbf{v}_{i}(t))-\phi(W_{i}\mathbf{r}_{i-1}(t))$ will contain a part
of the noise fluctuations of $\mathbf{r}_{i-1}(t)$, therefore causing the pre-
and post-synaptic term in the plasticity rule to be correlated.

A simple approach to remove the noise correlation in the plasticity rule for
the forward weights, is either to remove the noise fluctuations in the pre- or
post-synaptic term, or both. We choose to remove the noise fluctuations in the
pre-synaptic term, by using a low-pass filtered version of $\mathbf{r}_{i-1}$:

$\displaystyle\bar{\mathbf{r}}_{i}(t)$
$\displaystyle=\frac{1}{\tau_{f}}\int_{0}^{t}\exp\big{(}\frac{1}{\tau_{f}}(t-t^{\prime})\big{)}\mathbf{r}_{i}(t^{\prime})\text{d}t^{\prime},$
(92)

with $\tau_{f}$ the filtering time constant. This low-pass filtering can be
interpreted as a leaky biophysical mechanism that accumulates the pre-synaptic
input used for the plasticity of the synapse. As the pre-synaptic term of the
forward plasticity dynamics no longer contains the noise fluctuations, the
updates will no longer be biased.

##  Appendix D Simulation and algorithms

###  D.1 Simulating the dynamics

We use the Euler-Maruyama method to simulate the stochastic differential
equations (SDE), which simulates a (non)linear SDE

$\displaystyle\frac{\text{d}}{\text{d}t}\mathbf{x}(t)=\mathbf{f}(\mathbf{x},t)+\mathbf{g}(\mathbf{x},t)\boldsymbol{\xi}(t)$
(93)

as follows:

$\displaystyle\mathbf{x}[m+1]=\mathbf{x}[m]+\mathbf{f}(\mathbf{x}[m],m\Delta
t)\Delta t+\mathbf{g}(\mathbf{x}[m],m\Delta t)\sqrt{\Delta
t}\Delta\boldsymbol{\beta}_{m},$ (94)

with $\Delta\boldsymbol{\beta}$ a Gaussian random variable with distribution
$\mathcal{N}(\boldsymbol{0},I)$. Note that the noise term has stepsize
$\sqrt{\Delta t}$ instead of $\Delta t$, as Brownian motion has a variance of
$t$ and not $t^{2}$.

Applying the Euler-Maruyama method to the SDEs of Section C.1 defining the
neural and controller dynamics gives us the following discrete update
equations:

$\displaystyle\boldsymbol{\epsilon}_{i}[m+1]$
$\displaystyle=\boldsymbol{\epsilon}_{i}[m]+\frac{1}{\tau_{\epsilon}}\big{(}-\boldsymbol{\epsilon}_{i}[m]\Delta
t+\sqrt{\Delta t}\Delta\boldsymbol{\beta}_{i,m}\big{)}$ (95)
$\displaystyle\mathbf{v}^{\mathrm{fb}}_{i}[m]$
$\displaystyle=Q_{i}\mathbf{u}[m]+\sigma\boldsymbol{\epsilon}_{i}[m]$ (96)
$\displaystyle\mathbf{u}^{\mathrm{int}}[m+1]$
$\displaystyle=\mathbf{u}^{\mathrm{int}}[m]+(\mathbf{e}[m]-\alpha\mathbf{u}^{\mathrm{int}}[m])\frac{\Delta
t}{\tau_{u}}$ (97) $\displaystyle\mathbf{u}[m]$
$\displaystyle=k\mathbf{e}[m]+\mathbf{u}^{\mathrm{int}}[m]$ (98)
$\displaystyle\mathbf{v}_{i}[m+1]$
$\displaystyle=\mathbf{v}_{i}[m]+(-\mathbf{v}_{i}[m]+W_{i}\phi(\mathbf{v}_{i-1}[m])+\mathbf{v}^{\mathrm{fb}}_{i}[m])\frac{\Delta
t}{\tau_{v}}$ (99)

And for the weight updates:

$\displaystyle\bar{\mathbf{u}}[m+1]$
$\displaystyle=\bar{\mathbf{u}}[m]+\frac{1}{\tau_{f}}(-\bar{\mathbf{u}}[m]+\mathbf{u}[m])\Delta
t$ (100) $\displaystyle\bar{\mathbf{r}}_{i}[m+1]$
$\displaystyle=\bar{\mathbf{r}}_{i}[m]+\frac{1}{\tau_{f}}(-\bar{\mathbf{r}}_{i}[m]+\mathbf{r}_{i}[m])\Delta
t$ (101) $\displaystyle\Delta W_{i}$ $\displaystyle=\frac{\Delta
t}{\tau_{W}}\sum_{m}\big{(}\mathbf{r}_{i}[m]-\phi(W_{i}\mathbf{r}_{i-1}[m])\big{)}\bar{\mathbf{r}}_{i-1}[m]^{T}$
(102) $\displaystyle\Delta Q_{i}$ $\displaystyle=\frac{\Delta
t}{\tau_{Q}}\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{L-i}\sum_{m}-\mathbf{v}^{\mathrm{fb}}_{i}[m](\mathbf{u}[m]-\bar{\mathbf{u}}[m])^{T}-M_{\max}\frac{\Delta
t}{\tau_{Q}}\beta Q_{i}$ (103)

with $M_{\max}$ the total number of simulation steps.

###  D.2 Pseudocode for the Strong-DFC algorithm

Algorithm 1 provides the pseudocode for the Strong-DFC algorithm, making use
of the discretized dynamics provided in the previous section. We use the
Euler-Maruyama discretization of the SDEs with some slight modifications.

First, in order to have direct control over the hyperparameter
$\tilde{\alpha}$, we rewrite the controller dynamics as follows:

$\displaystyle\mathbf{u}(t)=\mathbf{u}^{\text{int}}(t)+k\mathbf{e}(t),\quad\tau_{u}\frac{\text{d}}{\text{d}t}\mathbf{u}^{\text{int}}(t)=\mathbf{e}(t)-\tilde{\alpha}\mathbf{u}(t).$
(104)

Note that the underlying dynamics for $\mathbf{u}(t)$ are exactly the same as
in Eq. (2).

In order to better incorporate the layered structure of the network into the
discrete network dynamics, we update a layer at timestep $m+1$ with the
previous layer at timestep $m+1$ instead of $m$:
$\mathbf{v}^{\mathrm{ff}}_{i}[m+1]=W_{i}\phi(\mathbf{v}_{i-1}[m+1])+\mathbf{b}_{i}$,
with $\mathbf{b}_{i}$ the bias parameters of layer $i$. For small stepsizes
$\Delta t$ this has almost no effect, but for bigger stepsizes it has the
benefit of better reflecting the layered structure.

In order to better incorporate the controller loop into the discrete dynamics,
we use the updated control signal $\mathbf{u}[m+1]$ for updating the feedback
compartment: $\mathbf{v}^{\mathrm{fb}}_{i}[m+1]=Q_{i}\mathbf{u}[m+1]$, such
that the control error $\mathbf{e}[m]$ of the previous timestep is used for
control instead of the error of two timesteps ago. For small stepsizes $\Delta
t$ this has almost no effect, but for bigger stepsizes it has the benefit of
better reflecting the control interaction.

Following the exponential filtering approach for the low-pass filtering of
$\mathbf{r}$ and $\mathbf{u}$, we use the current value $[m+1]$ of
$\mathbf{r}$ and $\mathbf{u}$ instead of the previous one prescribed by the
forward Euler method. This reduces the time delay introduced by the low pass
filtering.

Algorithm 1 Pseudocode for the Strong-DFC algorithm on a single input sample

# Initialize layer activations to the feedforward activations and parameter
update buffers to zero

for $i$ in range($1,L$) do

$\mathbf{v}_{i}[1]=\mathbf{v}_{i}^{-}$

$\mathbf{r}_{i}[1]=\mathbf{r}_{i}^{-}$

$\Delta W_{i}=0$

$\Delta\mathbf{b}_{i}=0$

$\Delta Q_{i}=0$

end for

$\mathbf{u}^{\mathrm{int}}[1]=0$

for $m$ in range(1,$M_{\max}$) do

# Update controller:

if softmax classification with cross-entropy loss then

$\mathbf{e}[m]=\mathbf{p}^{*}-\mathrm{softmax}(\mathbf{r}_{L}[m])$ {with
$\mathbf{p}^{*}$ the soft target and $\mathbf{r}_{L}$ a linear output layer}

else if regression with L2 loss then

$\mathbf{e}[m]=\mathbf{r}_{L}^{*}-\mathbf{r}_{L}[m]$

end if

$\mathbf{u}^{\mathrm{int}}[m+1]=\mathbf{u}^{\mathrm{int}}[m]+\frac{\Delta
t}{\tau_{u}}(\mathbf{e}[m]-\tilde{\alpha}\mathbf{u}[m])$

$\mathbf{u}[m+1]=\mathbf{u}^{\mathrm{int}}[m+1]+k\mathbf{e}[m]$

# Low-pass filter $\mathbf{u}$

if $m==1$ then

$\bar{\mathbf{u}}[m+1]=\mathbf{u}[m+1]$ {initialize $\bar{\mathbf{u}}$}

else

$\bar{\mathbf{u}}[m+1]=\bar{\mathbf{u}}[m]+\frac{1}{\tau_{f}}(-\bar{\mathbf{u}}[m]+\mathbf{u}[m+1])\Delta
t$

end if

# Update network:

for i in range(1,L) do

# Sample noise to update the exponentially filtered noise variable

$\Delta\boldsymbol{\beta}_{i,m}\sim\mathcal{N}(0,I)$

$\boldsymbol{\epsilon}_{i}[m+1]=\boldsymbol{\epsilon}_{i}[m]+\frac{1}{\tau_{\epsilon}}\big{(}-\boldsymbol{\epsilon}_{i}[m]\Delta
t+\sqrt{\Delta t}\Delta\boldsymbol{\beta}_{i,m}\big{)}$

# Inject noise in the feedback compartment and update the activations

$\mathbf{v}^{\mathrm{fb}}_{i}[m+1]=Q_{i}\mathbf{u}[m+1]+\sigma\boldsymbol{\epsilon}_{i}[m+1]$

$\mathbf{v}^{\mathrm{ff}}_{i}[m+1]=W_{i}\phi(\mathbf{v}_{i-1}[m+1])+\mathbf{b}_{i}$

$\mathbf{v}_{i}[m+1]=\mathbf{v}_{i}[m]+\frac{\Delta
t}{\tau_{v}}(-\mathbf{v}_{i}[m]+\mathbf{v}^{\mathrm{ff}}_{i}[m+1]+\mathbf{v}^{\mathrm{fb}}_{i}[m+1])$

$\mathbf{r}_{i}[m+1]=\phi(\mathbf{v}_{i}[m+1])$

# low-pass filter $\mathbf{r}_{i}$

$\bar{\mathbf{r}}_{i}[m+1]=\bar{\mathbf{r}}_{i}[m]+\frac{1}{\tau_{f}}(-\bar{\mathbf{r}}_{i}[m]+\mathbf{r}_{i}[m+1])\Delta
t$

# Buffer parameter updates:

$\Delta W_{i}=\Delta
W_{i}+\big{(}\mathbf{r}_{i}[m+1]-\phi(W_{i}\mathbf{r}_{i-1}[m+1])\big{)}\bar{\mathbf{r}}_{i-1}[m+1]^{T}$

$\Delta\mathbf{b}_{i}=\Delta\mathbf{b}_{i}+\mathbf{r}_{i}[m+1]-\phi(\mathbf{v}^{\mathrm{ff}}_{i}[m+1])$

$\Delta Q_{i}=\Delta
Q_{i}-\left(1+\frac{\tau_{v}}{\tau_{\epsilon}}\right)^{L-i}\mathbf{v}^{\mathrm{fb}}_{i}[m](\mathbf{u}[m+1]-\bar{\mathbf{u}}[m+1])^{T}-\beta
Q_{i}$

end for

end for

Update parameters with $\Delta W_{i}/M_{\max}$,
$\Delta\mathbf{b}_{i}/M_{\max}$ and $\Delta Q_{i}/M_{\max}$ using an optimizer
of choice

##  Appendix E Experiments

###  E.1 Alignment measures

Below we describe how the measures of Fig. 3 were computed.

#### Alignment of Strong-DFC updates with $\mathcal{H}$

Fig. 3B quantifies alignment between actual Strong-DFC updates and the
gradient of $\mathcal{H}$. For this, we vectorize the update matrices and
compute the angle between them in degrees.

#### Relative strength of feedback control

Fig. 3C quantifies the magnitude of the feedback input compared to the forward
drive, and is used to illustrate that our framework allows feedback to have
comparable magnitudes to the forward input. It is calculated as follows:

$\displaystyle\mathrm{ratio}_{\text{fb}/\text{ff}}=\frac{\|Q\mathbf{u}\|_{F}}{\|W\mathbf{r}\|_{F}},$
(105)

where $F$ indicates the Frobenius norm.

#### Condition 1

Fig. 3D quantifies the degree to which Condition 1 is satisfied. It is
computed through the Frobenius norm of the projection of $Q$ with the norm of
$Q$:

$\displaystyle\mathrm{ratio}_{\mathrm{Con1}}=\frac{\|P_{J_{ss}^{T}}Q\|_{F}}{\|Q\|_{F}}.$
(106)

where the projection $P_{J_{ss}^{T}}Q$ is given by:

$\displaystyle P_{J_{ss}^{T}}Q=J_{ss}^{T}(J_{ss}J_{ss}^{T})^{-1}J_{ss}Q.$
(107)

For an intuition of this metric, please refer to Appendix F.1 in Meulemans et
al. (2021).

###  E.2 Noise robustness measures

Below we further illustrate the robustness of Strong-DFC weight updates under
noisy conditions. Recall that this is a direct consequence of the main novelty
of Strong-DFC compared to standard DFC: the use of the true output target
$\mathbf{r}_{L}^{*}=\mathbf{r}_{L}^{\mathrm{true}}$ throughout training. For
this series of experiments (Figure S1) we vary the noise magnitude $\sigma$
added to the neural dynamics during training of the forward weights for DFC
and Strong-DFC. For a fair comparison between both methods, we low-pass filter
the presynaptic plasticity signal for debiasing the updates (11). In a
student-teacher regression setting with a nonlinear network of size 20-40-5,
we found that the performance of Strong-DFC is more robust to a wide range of
noise magnitudes compared to standard DFC. For non-noisy dynamics, we observed
that standard DFC can descend the loss landscape faster compared to Strong-DFC
in this setting, leading to a lower loss, as the Strong-DFC training has not
yet fully converged after 500 epochs. For the results of Table 2 in the main
manuscript, we used $\sigma=0.005$ for both methods.

Figure S1:  Noise robustness. Training loss after 500 epochs (training not
fully converged for Strong-DFC) for different noise magnitudes. The DFC
results for high noise values are missing since its training was unstable.

##  Appendix F Additional Results

###  F.1 Training Strong-DFC without re-visiting samples

To investigate whether Strong-DFC requires consistency or re-visiting training
samples multiple times to perform well, we performed the following experiment
using the student-teacher setting. Using a teacher of size 30-10-10-10-5 and a
student of size 30-50-50-50-5, we trained networks for 300 epochs, each
containing 500 training samples. In the standard setting, the 500 training
samples were identical across epochs. In the modified setting, each training
sample was independently drawn from the training distribution, so that no
sample was seen twice. This was achieved by feeding the teacher with a new set
of randomly-generated input samples at each epoch. Therefore, training
consisted of a total of 150000 samples, but in the modified setting all these
samples were unique. The results obtained are presented in Table S1.

Table S1: Effect of sample reuse when training Strong-DFC. Training setting | Train loss | Test loss
---|---|---
Standard | $1.945^{\pm 0.001}$ | $2.224^{\pm 0.002}$
Without sample reuse | $1.964^{\pm 0.004}$ | $2.226^{\pm 0.006}$

The results obtained are extremely similar, indicating that Strong-DFC does
not necessarily rely on re-visiting training samples many times to perform
well. However, it is important to note that Strong-DFC needs i.i.d. samples
and will not be immune to forgetting and interference issues when exposed to
correlated sequences of training samples; solving these issues will likely
require continual learning techniques, which is out of the scope for this
work.