supervised predictive coding updates with standard teacher clamping to
gradient updates on the least-control objective (3) thereby removing the need
for the weak-nudging limit.

Supervised predictive coding uses a probabilistic latent variable model to
compute the likelihood $p(y\mid x;\theta)$ of the output $y$ given the sensory
input $x$. More specifically, it uses a hierarchical Gaussian model, a
particular type of directed acyclic graphical model, with the following
factorized distribution

$\begin{split}p(z_{1},\ldots,z_{L-1},y\mid x;\theta)&=p(y\mid
z_{L-1})\prod_{l=1}^{L-1}p(\phi_{l}\mid z_{l-1})\\\
p(z_{l}\mid\phi_{l-1})&=\mathcal{N}(z_{l};~{}W_{l}\sigma(z_{l-1}),\Sigma_{l}))\\\
p(y\mid
z_{L-1})&=\mathcal{N}(y;~{}W_{L}\sigma(z_{L-1},\Sigma_{y}))\end{split}$ (94)

with latent variables $z=\\{z_{l}\\}_{l=1}^{L-1}$, parameters
$\theta=\\{W_{l}\\}_{l=1}^{L}$, $\mathcal{N}(v;\mu,\Sigma)$ a Gaussian
distribution over a variable $v$ with mean $\mu$ and correlation matrix
$\Sigma$, and where we take $\sigma(z_{0})=x$. The goal of training this model
is to change the model parameters $W_{l}$ to maximize the marginal likelihood
for input-output pairs $(x,y)$:

$p(y\mid x;\theta)=\int p(z,y\mid x;\theta)\mathrm{d}z$ (95)

As optimizing this marginal directly is intractable, predictive coding uses
the variational expectation maximization method neal_view_1998,
tzikas_variational_2008, bogacz_tutorial_2017, friston_free_2006. This method
uses a variational distribution $q(z;\phi)$ with parameters $\phi$ to
approximate the posterior $p(z\mid x,y)$ and then defines the free energy
$\mathcal{F}$ as an upper bound on the negative log-likelihood:

$\begin{split}\mathcal{F}(\phi,\theta)&=-\mathbb{E}_{q}[\log p(y,z\mid
x;\theta)]+\mathbb{E}_{q}[\log q(z;\phi)]\\\ &=-\log p(y\mid
x;\theta)+D_{\mathrm{KL}}(q(z;\phi)\|p(z\mid x,y;\theta))\geq-\log p(y\mid
x;\theta),\end{split}$ (96)

with $D_{\mathrm{KL}}$ the KL divergence which is always positive. Now the
variational expectation maximization method proceeds in two alternating
phases. In the expectation phase, it minimizes the free energy $\mathcal{F}$
w.r.t. $\phi$ for the current parameter setting of $\theta$, i.e. it finds a
good approximation $q(z;\phi)$ for the posterior. In the maximization phase,
it minimizes the free energy w.r.t. $\theta$ for the updated parameter setting
of $\phi$.

To be able to perform the two phases using simple neural dynamics and local
update rules, predictive coding uses the Dirac delta distribution
$q(z;\phi)=\delta(z-\phi)$ as the variational distribution
bogacz_tutorial_2017. Now, the entropy term $\mathbb{E}_{q}[\log p(y,z\mid
x;\theta)]$ is independent of $\phi$ and hence we ignore it in the expectation
and maximization step.222Note however that this entropy term is infinite for
the delta distribution, hence for making the derivation rigorous, certain
limits need to be taken. For example, friston_free_2006 uses a Laplace
approximation instead of the delta distribution. Taken together, this result
in the following free energy (without entropy term):

$\begin{split}\mathcal{F}(\phi,\theta)&=-\log p(y,\phi\mid x;\theta)\\\
&=\rho\frac{1}{2}\lVert
y-W\sigma(\phi_{L-1})\rVert^{2}+\frac{1}{2}\sum_{l=1}^{L-1}\lVert\phi_{l}-W_{l}\sigma(\phi_{l-1})\rVert^{2}+C\end{split}$
(97)

with $C$ a constant independent of $\phi$ and $\theta$, $x=\sigma(\phi_{0})$,
and where we took for simplicity $\Sigma_{l}=\mathrm{Id}$ and
$\Sigma_{y}=\rho^{-1}\mathrm{Id}$. For the expectation step, predictive coding
interprets the parameters $\phi=\\{\phi_{l}\\}_{l=1}^{L-1}$ as a prediction
neuron population and uses the gradient flow on $\mathcal{F}$ w.r.t. $\phi$ as
neural dynamics to find a $\phi_{*}$ that minimizes $\mathcal{F}$. For the
maximization step, predictive coding performs a gradient step on
$\mathcal{F}(\phi_{*},\theta)$ w.r.t. $\theta$, evaluated at the steady-state
neural activity $\phi_{*}$. Thus,

$\begin{split}\dot{\phi}&=-\frac{\partial\mathcal{F}}{\partial\phi}(\phi,\theta)^{\top},\\\
\Delta\theta&=-\frac{\partial\mathcal{F}}{\partial\theta}(\phi_{*},\theta)^{\top},\end{split}$
(98)

with the change $\Delta\theta$ applied after the neural dynamics has
converged, i.e., when $\dot{\phi}=0$.

####  S4.2 Equivalence between the free and augmented energies

To link supervised predictive coding to the least-control principle, we
further manipulate the free energy of Eq. 97 to relate it to the augmented
energy $F=\frac{1}{2}\lVert f(\phi,\theta)\rVert^{2}+\beta L$ introduced in
Section 2.4. For this, we introduce a ‘ghost layer’ $z_{L}$ in between
$z_{L-1}$ and $y$, that exists solely for the purpose of analyzing the free
energy. The joint probability distribution is now given by

$\begin{split}p(z,y\mid x;\theta)&=p(y\mid z_{L})\prod_{l=1}^{L}p(z_{l}\mid
z_{l-1})\\\ p(z_{l}\mid
z_{l-1})&=\mathcal{N}(z_{l};~{}W_{l}\sigma(z_{l-1}),\mathrm{Id})\\\ p(y\mid
z_{L})&=\mathcal{N}(y;~{}z_{L},\beta^{-1}\mathrm{Id})\end{split}$ (99)

We recover the same marginal likelihood $p(y\mid x;\theta)$ as before if we
have that $\rho^{-1}=\beta^{-1}+1$, assuming that $\rho^{-1}>1$. Using again
the dirac delta distribution as variational distribution, we get the following
augmented free energy

$\begin{split}\mathcal{F}(\phi,\theta,\beta)&=\log p(\phi,y\mid
x;\theta)=\frac{1}{2}\sum_{l=1}^{L}\lVert\phi_{l}-W_{l}\sigma(\phi_{l-1})\rVert^{2}+\beta\frac{1}{2}\lVert
y-\phi_{L}\rVert^{2}\\\
&=\frac{1}{2}\lVert\phi-W\sigma(\phi)-Ux\rVert^{2}+\beta L(D\phi)\\\
W&=\begin{bmatrix}0&0&0&0\\\ W_{2}&0&0&0\\\ 0&\ddots&0&0\\\
0&0&W_{L}&0\end{bmatrix},~{}~{}~{}U=\begin{bmatrix}W_{1}\\\
0\end{bmatrix},~{}~{}~{}D=\begin{bmatrix}0&\mathrm{Id}\end{bmatrix},\end{split}$
(100)

with $\phi$ the concatenation of $\\{\phi_{l}\\}_{l=1}^{L}$, $L$ the squared
error loss and where we ignored the constant term $C$. Hence, we see that the
free energy used in predictive coding is equivalent to the augmented energy
$F=\frac{1}{2}\lVert f(\phi,\theta)\rVert^{2}+\beta L$ introduced in Section
2.4. Consequently, we can use the least-control principle to characterize
predictive coding in the limit of $\beta\to\infty$.

####  S4.3 Learning with teacher clamping results in gradients on the least-
control objective

Whittington and Bogacz [36] connected supervised predictive coding to
gradient-based optimization of the loss $L$, when the variance of the output
layer $y$ goes to infinity. In this limit, there is only an infinitesimal
effect of the teaching signal on the rest of the network, arising by clamping
$y$ towards $y^{\mathrm{true}}$, as the model ‘does not trust’ the output
layer $y$ by construction. This weak nudging limit is captured in the
augmented energy of Eq. (100) as the limit of $\beta\to 0$, and can be linked
to gradient-based optimization of the loss $L$ for more general energy
functions beyond predictive coding [24].

Equipped with the least-control principle, we can now relate predictive coding
to gradient-based optimization on the least-control objective of Eq. 3,
without the need for infinite variance in the output layer. In Sections 2.4
and S2.5, we show that when the neural dynamics minimize the augmented energy
$F(\phi,\theta,\beta)$ for $\beta\to\infty$, the parameter update
$\Delta\theta$ of Eq. 98 follows the gradient on the least-control objective
(3). This perfect control limit of $\beta\to\infty$ corresponds to the
conventional teacher clamping setting, where the output $y$ has a finite
variance (without loss of generality, $\rho=1$ in this case in our derivation,
as $\rho^{-1}=\beta^{-1}+1$).

####  S4.4 The optimal control interpretation of predictive coding leads to
flexible dynamics.

Now that we connected predictive coding to the least-control principle, we can
further investigate predictive coding through the lens of optimal control. In
predictive coding, the neural dynamics follow the gradient flow on the
augmented free energy $\mathcal{F}$ of Eq. (100):

$\begin{split}\dot{\phi}&=-\frac{\partial\mathcal{F}}{\partial\phi}(\phi,\theta,\beta)^{\top}\\\
&=-\phi+W\sigma(\phi)+Ux+\sigma^{\prime}(\phi)W^{\top}\big{(}\phi-W\sigma(\phi)-Ux\big{)}-\beta
D^{\top}\nabla_{y}L(D\phi)\end{split}$ (101)

with $\sigma^{\prime}(\phi)=\partial_{\phi}\sigma(\phi)$. By comparing the
above equation to $\dot{\phi}=f(\phi,\theta)+\psi$, we have the optimal
control
$\psi=\sigma^{\prime}(\phi)W^{\top}\big{(}\phi-W\sigma(\phi)-Ux\big{)}-\beta
D^{\top}\nabla_{y}L(D\phi)$ For $\beta\to\infty$, the last term $\beta
D^{\top}\nabla_{y}L(D\phi)$ corresponds to an infinitely fast proportional
control that clamps the output $y=\phi_{L}$ to the teacher
$y^{\mathrm{true}}$. The term
$\sigma^{\prime}(\phi)W^{\top}\big{(}\phi-W\sigma(\phi)-Ux\big{)}$ then
optimally propagates the control at the output level towards the rest of the
network. In predictive coding, the term $\phi-W\sigma(\phi)-Ux$ is usually
interpreted as an error neuron population $\epsilon$. Hence, in the optimal
control view, the error neuron population $\epsilon$ optimally controls the
prediction neuron population $\phi$ via
$\psi=\sigma^{\prime}(\phi)W^{\top}\epsilon-\beta D^{\top}\nabla_{y}L(D\phi)$,
to reach a controlled network state that exactly matches the output target
$y^{\mathrm{true}}$, while having the smallest possible error signals
$\lVert\epsilon\rVert^{2}=\lVert\phi-W\sigma(\phi)-Ux\rVert^{2}$. Note that at
equilibrium, we have that $\psi=-f(\phi,\theta)=\epsilon$. This optimal
control view is further corroborated by isolating the error neurons $\epsilon$
at the equilibrium of the neural dynamics (101), leading to
$\epsilon_{*}=\psi_{*}=-(I-W\sigma^{\prime}(\phi_{*}))^{-\top}D^{\top}\beta\nabla_{y}L(D\phi_{*})$,
which satisfies the column space condition of Theorem 2 with
$u_{*}=-\beta\nabla_{y}L(D\phi_{*})$, confirming that the error neurons steady
state $\psi_{*}$ is an optimal control.

Importantly, using the least-control principle, we can go beyond the gradient-
flow dynamics on $F$ and generalize the predictive coding framework to more
flexible neural dynamics. First, instead of clamping the output of the network
towards the teacher value with an infinitely fast proportional controller, one
can use more general output controllers that satisfy the equilibrium condition
$\alpha u_{*}+\nabla_{y}L(D\phi)=0$ for $\alpha\to 0$, such as an integral
controller. Next, as we identified the error neurons
$\epsilon=\phi-W\sigma(\phi)-Ux$ as being an optimal control
$\epsilon_{*}=\psi_{\scalebox{0.42}{$\bigstar$}}$ at equilibrium, we can use
any dynamics satisfying Theorem 2 to compute these optimal error neurons and
resulting neural dynamics $\dot{\phi}$. For example, we can use the inversion
dynamcis of Eq. 8 and 14 to dynamically compute the error neuron activity.
Strikingly, at the steady state, the error neurons computed with these
inversion dynamics are indistinguishable from the error neurons computed with
the energy-based dynamics of Eq. 101, even though the underlying dynamics and
hence circuit implementation are completely different. This ‘circuit
invariance’ under the optimality condition of Theorem 2 opens new routes
towards finding cortical circuits implementing predictive coding, which is a
topic of high relevance in neuroscience [58]. Furthermore, the least-control
principle allows unifying predictive coding with other existing theories for
learning in the brain [30, 17] by uncovering possible equivalences of the
underlying credit assignment techniques.

####  S4.5 Beyond feedforward neural networks.

The underlying theory of predictive coding, which starts from an acyclic
graphical model, is tailored towards feedforward neural networks. However,
starting from the augmented energy $F$ of Eq. (100), we can generalize
predictive coding to equilibrium recurrent neural networks with arbitrary
synaptic connectivity matrices $W$ and $U$. Although the link with
probabilistic latent variable models is not straightforwardly extendable to
equilibrium RNNs, the predictive coding interpretation based on error and
prediction neurons remains valid. To the best of our knowledge, our
experiments testing this new setting of predictive coding for equilibrium RNNs
are the first of their kind.

###  Appendix S5 Remarks on the least-control principle when using multiple
data points

Without loss of generality, we consider a single data point in the formulation
of the least-control principle in Section 2, for clarity of presentation.
Here, we discuss in more detail how multiple data points can be incorporated
in the least-control principle.

####  S5.1 Loss defined over finite sample of data points

In many learning problems, as in supervised learning, the loss $L$ is defined
over a finite set of data samples:

$\displaystyle L(\phi)=\sum_{b=1}^{B}L^{b}(\phi).$ (102)

As the losses associated to different samples are different, the corresponding
teaching signal will impact the state differently. We therefore consider
$\phi$ to be a concatenation of all the sample-specific states
$\\{\phi^{b}\\}_{b=1}^{B}$. The loss is then equal to

$L(\phi)=\sum_{b=1}^{B}L^{b}(\phi^{b}).$ (103)

Similarly, we define $f(\phi,\theta)$ as a concatenation of all
$f(\phi^{b},\theta)$, and $\psi$ as a concatenation of all the controls
$\psi^{b}$. We can then apply the least-control principle on this concatenated
quantities.

As the least-control objective of Eq. 3 can be rewritten as the sum
$\sum_{b}\lVert{\psi_{\scalebox{0.42}{$\bigstar$}}}^{b}\rVert^{2}$, and has
constraints that do not interact in between different datapoints $b$, its
gradient can also be rewritten as a sum over the datapoints:

$\displaystyle\frac{\mathrm{d}\mathcal{H}}{\mathrm{d}\theta}(\theta)=\sum_{b=1}^{B}\frac{\mathrm{d}\mathcal{H}^{b}}{\mathrm{d}\theta}(\theta)$
(104)

with $\mathcal{H}^{b}(\theta)$ the least-control objective for a single data
point $b$. It follows that one can use standard stochastic or mini-batch
optimization methods to minimize the least-control objective
$\mathcal{H}(\theta)$.

####  S5.2 Loss defined by expectation over an infinite inputspace

When there exists an infinite number of data samples (e.g. a continuous input
space), the loss can be defined as an expectation over this infinity of data
samples:

$\displaystyle L(\phi)=\mathbb{E}_{b}[L_{b}(\phi)]$ (105)

One approach to incorporate this case in the least-control principle is to
first sample a finite amount of data samples, and then apply the arguments of
Section S5.1.

Another approach is to define the least-control objective as an expectation
over the infinity of samples:

$\begin{split}\mathcal{H}(\theta)&=\mathbb{E}_{b}[\mathcal{H}^{b}(\theta)]\\\
\mathcal{H}^{b}(\theta)&=\min_{\phi^{b},\psi^{b}}\lVert\phi^{b}\rVert^{2}\quad\text{s.t.}\enspace
f(\phi^{b},\theta)+\psi^{b}=0,~{}~{}\nabla_{y}L(h(\phi^{b})=0\end{split}$
(106)

Then, under standard regularity conditions for exchanging the gradient and
expectation operator, we can sample gradients of this least-control objective:

$\frac{\mathrm{d}\mathcal{H}}{\mathrm{d}\theta}(\theta)=\mathbb{E}\left[\frac{\mathrm{d}\mathcal{H}^{b}}{\mathrm{d}\theta}(\theta)\right]$
(107)

with $\mathrm{d}_{\theta}\mathcal{H}^{b}(\theta)$ given by Theorem 1.

###  Appendix S6 The least-control principle for equilibrium RNNs

Here, we provide additional details on the equilibrium recurrent neural
network models, the derivation of the local learning rules and the various
controller designs.

####  S6.1 Equilibrium recurrent neural network model specifications

We use an equilibrium RNN with the following free dynamics:

$\displaystyle\dot{\phi}=f(\phi,\theta)=-\phi+W\sigma(\phi)+Ux,$ (108)

with $\phi$ the neural activities of the recurrent layer, W the recurrent
synaptic weight matrix, U the input weight matrix and $\sigma$ the activation
function. We evaluate the performance of this RNN at equilibrium, by selecting
a set of output neurons and measuring its loss $L(D\phi^{*})$, with
$D=[0~{}\mathrm{Id}]$ a fixed decoder matrix.

Learning a decoder matrix. In practice, we can achieve better performance if
we learn a decoder matrix to map the recurrent activity $\phi$ to the output
$y$, instead of taking a fixed decoder $D=[0~{}\mathrm{Id}]$. To prevent the
loss $L(D\phi)$ from depending on the learned parameters $\theta$, we augment
the recurrent layer $\phi$ with an extra set of output neurons, and include a
learned decoder matrix $\tilde{D}$ inside the augmented recurrent weight
matrix $\bar{W}$, leading to the following free dynamics:

$\begin{split}\dot{\bar{\phi}}&=f(\bar{\phi},\theta)=-\bar{\phi}+\bar{W}\sigma(\bar{\phi})+\bar{U}x,\\\
\bar{W}&=\begin{bmatrix}W&0\\\
\tilde{D}&0\end{bmatrix},~{}~{}~{}\bar{U}=\begin{bmatrix}U\\\
0\end{bmatrix},\end{split}$ (109)

with $\bar{\phi}$ the concatenation of $\phi$ with a set of output neurons and
$\theta=\\{W,U,\tilde{D}\\}$ the set of learned parameters. To map the
augmented recurrent layer $\bar{\phi}$ to the output $y$, we again use a fixed
decoder $D=[0~{}\mathrm{Id}]$ that selects the ‘output neurons’ in
$\bar{\phi}$. At equilibrium, we now have
$y=D\bar{\phi}=\tilde{D}\sigma(\phi)$. Hence, we see that at equilibrium, the
augmented system of Eq. 109 with a fixed decoder $D$ and hence a loss
$L(D\bar{\phi})$ independent of $\theta$ is equivalent to the original system
of Eq. 108 with a learned decoder $D$ and loss $L(D\sigma(\phi))$. We can now
use the least-control principle on the augmented system of Eq. 109 to learn
the weight matrices $\theta=\\{W,U,\tilde{D}\\}$, which we show in the next
section.

Including biases. In practice, we equip the neurons with a bias parameter $b$,
leading to the following free dynamics

$\displaystyle\dot{\phi}=f(\phi,\theta)=-\phi+W\sigma(\phi)+Ux+b,$ (110)

with the set of learned parameters $\theta=\\{W,U,b\\}$. This can be extended
to the augmented system of Eq. 109 by adding the augmented bias parameters
$\bar{b}$, which is a concatenation of the original bias $b$ and the decoder
bias $\tilde{b}$. At equilibrium, we now have
$y=D\bar{\phi}=\tilde{D}\sigma(\phi)+\tilde{b}$. For ease of notation, we omit
the bias terms in most of the derivations in this paper.

Link with feedforward neural networks. When the recurrent weight matrix has a
lower block-diagonal structure, the equilibrium of the RNN corresponds to a
conventional deep feedforward neural network. More specifically, taking

$\displaystyle W=\begin{bmatrix}0&0&0&0\\\ W_{2}&0&0&0\\\ 0&\ddots&0&0\\\
0&0&W_{L}&0\end{bmatrix},~{}~{}~{}U=\begin{bmatrix}W_{1}\\\
0\end{bmatrix},~{}~{}~{}D=\begin{bmatrix}0&\mathrm{Id}\end{bmatrix},$ (111)

gives the following feedforward mappings at equilibrium

$\displaystyle\phi_{l}=W_{l}\sigma(\phi_{l-1})+b_{l},~{}~{}1\leq l\leq L$
(112)

where we structure $\phi$ into $L$ layers $\phi_{l}$, and take
$\sigma(\phi_{0})=x$ and $y=D\phi=\phi_{L}$. Hence, we see that the
equilibrium RNN model also includes the conventional deep feedforward neural
network. Similarly, by further structuring $W$ with blocks of Toeplitz
matrices, the equilibrium point corresponds with a convolutional feedforward
neural network.

Fixed point iterations. For computational efficiency, we use fixed point
iterations for finding the equilibrium of the free dynamics (108), which we
need for the recurrent backpropagation baseline, and for evaluating the test
performance of the trained equilibrium RNN.

$\displaystyle\phi_{\mathrm{next}}=W\sigma(\phi)+Ux$ (113)

Figure S2: Comparison of the different types of controllers we use on a
2-hidden-layer feedforward neural network in a supervised learning setting.
The light blue arrows represent the weights of the network that are used for
inference, the light green ones the feedback weights (the $Q$ mapping of
Eq.6), the dark green arrows identity mappings and the dark blue ones
transformations through the activation function. Dashed arrows indicates the
use of the activation function derivative. The grey rectangle shows how a
state neuron $\phi$ and its corresponding control neuron $\psi$ are connected,
along with how they are connected to the layers above and below. For example,
a light blue bottom-up arrow on the left indicates that the state neurons
receive bottom-up input from the state neurons of the layer below through the
inference weights. Alternatively, a light green top-down arrow from the right
to the left shows that the control signal from the layer above influence the
state neurons through the feedback weights. (A) Linear feedback (LF). The
output control signal $u$ is directly propagated to the hidden layers through
the linear mapping $\psi=Qu$. (B) Dynamic inversion (DI). The control signal
is propagated through the hidden layers by leveraging layer-wise feedback
mechanisms $\dot{\psi}=-\psi+\sigma^{\prime}(\phi)S\psi+D^{\top}u$ and
influences the state variable jointly with bottom up input:
$\dot{\phi}=-\phi+W\sigma(\phi)+Ux+\psi$. (C) Energy-based dynamics (EBD). The
control signal is here equal to the difference between the current state and
the bottom-up input ($\psi=\phi-W\sigma(\phi)-Ux$) and is indirectly
propagated to the remote layers through the state variables:
$\dot{\phi}=-\psi+\sigma^{\prime}(\phi)S\psi+D^{\top}u$.

####  S6.2 Deriving a local learning rule for equilibrium RNNs

The least-control principle prescribes to augment the RNN dynamics of Eq. 108
with an optimal controller $\psi$ that leads the neural dynamics to a
controlled equilibrium $\phi_{*}$, where the output neurons $y_{*}=D\phi_{*}$
minimize the loss $L$, while being of minimum norm
$\lVert\psi_{\scalebox{0.42}{$\bigstar$}}\rVert$. This optimal control
contains a useful learning signal, which we leverage in Theorem 1 to derive
local first-order parameter updates
$\Delta\theta=-\mathrm{d}_{\theta}\mathcal{H}(\theta)={\psi_{\scalebox{0.42}{$\bigstar$}}}^{\top}\partial_{\theta}f(\phi_{\scalebox{0.42}{$\bigstar$}},\theta)$,
with $(\psi_{\scalebox{0.42}{$\bigstar$}},\phi_{\scalebox{0.42}{$\bigstar$}})$
the minimizers of the least-control problem of Eq. 3, i.e. the optimal control
and optimally controlled state, respectively. In the next sections, we discuss
in detail how to compute
$(\psi_{\scalebox{0.42}{$\bigstar$}},\phi_{\scalebox{0.42}{$\bigstar$}})$ for
the equilibrium RNN using neural dynamics. Here, we apply the first-order
update $\Delta\theta$ to the equilibrium RNN to arrive at the local synaptic
updates of Eq. 12.

Derivation of the local update rules. To avoid using tensors in
$\partial_{\theta}f(\phi_{\scalebox{0.42}{$\bigstar$}},\theta)$, we vectorize
the weight matrices $W$ and $U$ and use Kronecker products to rewrite the free
dynamics of Eq. 108:333For this, we use the following property:
$\mathrm{vec}(ABC)=(C^{\top}\otimes A)\mathrm{vec}(B)$, with $\mathrm{vec}$
the vectorization operator that concatenates the columns of a matrix and
$\otimes$ the Kronecker product.

$\displaystyle\dot{\phi}=-\phi+(\sigma(\phi)^{\top}\otimes\mathrm{Id})\vec{W}+(x^{\top}\otimes\mathrm{Id})\vec{U}+b$
(114)

with $\otimes$ the Kronecker product, $\mathrm{Id}$ the identity matrix and
$\vec{W}$ and $\vec{U}$ the vectorized form (concatenated columns) of $W$ and
$U$, respectively. Taking $\theta$ the concatenation of $\vec{W}$, $\vec{U}$
and $b$, we get

$\displaystyle\Delta\theta={\psi_{\scalebox{0.42}{$\bigstar$}}}^{\top}\frac{\partial
f}{\partial\theta}(\phi_{\scalebox{0.42}{$\bigstar$}},\theta)={\psi_{\scalebox{0.42}{$\bigstar$}}}^{\top}\begin{bmatrix}(\sigma(\phi)^{\top}\otimes\mathrm{Id})&(x^{\top}\otimes\mathrm{Id})&\mathrm{Id}\end{bmatrix}$
(115)

We can rewrite the resulting update in matrix form again, leading to the
following parameter updates for the RNN:

$\displaystyle\Delta
W=\psi_{\scalebox{0.42}{$\bigstar$}}\sigma(\phi_{\scalebox{0.42}{$\bigstar$}})^{\top},~{}~{}~{}\Delta
U=\psi_{\scalebox{0.42}{$\bigstar$}}x^{\top},~{}~{}~{}\Delta
b=\psi_{\scalebox{0.42}{$\bigstar$}}.$ (116)

Now assuming we have a controller $\psi$ that leads the dynamics to the
optimally controlled state
$(\psi_{*},\phi_{*})=(\psi_{\scalebox{0.42}{$\bigstar$}},\phi_{\scalebox{0.42}{$\bigstar$}})$,
we arrive at the updates of Eq. 12. Note that we can use these updates in any
gradient-based optimizer. Applying a similar derivation to the augmented RNN
dynamics of Eq. 109, we arrive at the following updates for
$\theta=\\{W,U,b,\tilde{D},\tilde{b}\\}$:

$\Delta
W=\psi_{\scalebox{0.42}{$\bigstar$}}\sigma(\phi_{\scalebox{0.42}{$\bigstar$}})^{\top},~{}~{}~{}\Delta
U=\psi_{\scalebox{0.42}{$\bigstar$}}x^{\top}~{}~{}~{},\Delta\tilde{D}=\tilde{\psi}_{\scalebox{0.42}{$\bigstar$}}\sigma(\phi_{\scalebox{0.42}{$\bigstar$}})^{\top},~{}~{}~{}\Delta
b=\psi_{\scalebox{0.42}{$\bigstar$}},~{}~{}~{}\Delta\tilde{b}=\tilde{\psi}_{\scalebox{0.42}{$\bigstar$}}$
(117)

with $\tilde{\psi}$ the control applied on the set of output neurons within
$\bar{\phi}$, and $\psi$ the control applied to the set of recurrent neurons
$\phi$ (i.e. $\bar{\phi}$ without the output neurons).

Controller dynamics. Theorem 2 shows that we can compute an optimal control
$\psi_{\scalebox{0.42}{$\bigstar$}}$ by using any controller $\psi$ that leads
the system to a controlled equilibrium

$\displaystyle 0=f(\phi_{*},\theta)+Q(\phi_{*},\theta)u_{*},~{}~{}~{}0=-\alpha
u_{*}-\frac{\partial L}{\partial y}(h(\phi_{*}))^{\top}$ (118)

in the limit of $\alpha\to 0$ and with $Q(\phi_{*},\theta)$ satisfying the
column space condition

$\displaystyle\mathrm{col}\left[Q(\phi_{*},\theta)\right]=\mathrm{row}\left[\frac{\partial
h}{\partial\phi}(\phi_{*})\left(\frac{\partial
f}{\partial\phi}(\phi_{*},\theta)\right)^{-1}\right].$ (119)

Under these conditions, we have that $Q(\phi_{*},\theta)u_{*}$ is an optimal
control $\psi_{\scalebox{0.42}{$\bigstar$}}$. Applied to the equilibrium RNN
setting, this result in the following column space condition.

$\displaystyle\mathrm{col}\left[Q(\phi_{*},\theta)\right]=\mathrm{row}\left[D\big{(}\mathrm{Id}-W\sigma^{\prime}(\phi_{*})\big{)}^{-1}\right],$
(120)

with $\sigma^{\prime}(\phi):=\partial_{\phi}\sigma(\phi)$. In the next
sections, we propose various control dynamics $\dot{\psi}$ that satisfy these
conditions at equilibrium and hence lead to an optimal control
$\psi_{\scalebox{0.42}{$\bigstar$}}$.

####  S6.3 A controller with direct linear feedback

We start with the simplest controller, that broadcasts the output controller
$u$ to the recurrent neurons $\phi$ with a direct linear mapping $Q$,
resulting in the following dynamics:

$\displaystyle\tau\dot{\phi}=-\phi+W\sigma(\phi)+Ux+Qu,~{}~{}~{}\tau_{u}\dot{u}=-\frac{\partial
L}{\partial y}(D\phi)^{\top}-\alpha u,$ (121)

with $\tau$ and $\tau_{u}$ the time constants of the neural dynamics and
output controller dynamics, respectively. We visualize an example of such a
direct linear feedback controller on a feedforward neural network in Figure
S2.A. We assume without further discussion that the output controller $u$
which integrates the output error $\partial_{y}L$ is implemented by some
neural integrator circuit, a topic which has been extensively studied in ref.
goldman_neural_2010.

#####  S6.3.1 A single-phase learning scheme for updating the feedback weights

Theorem 2 guarantees that this dynamics converges to the optimal control
$\psi_{\scalebox{0.42}{$\bigstar$}}$ if the columnspace condition of Eq. 120
is satisfied. As $Q$ is now a linear mapping, this condition cannot be exactly
satisfied for all samples in the dataset, as the right-hand side of Eq. 120 is
data-dependent, whereas $Q$ is not. Still, we learn the feedback weights $Q$
simultaneously with the other weights $W$ and $U$, to approximately fulfill
this column space condition, without requiring separate phases. For this, we
add noise to the dynamics and we adapt a recently proposed feedback learning
rule for feedforward neural networks [30, 51] towards our RNN setup, leading
to the following system dynamcis and simple anti-Hebbian plasticity dynamics:

$\displaystyle\tau\dot{\phi}$
$\displaystyle=-\phi+W\sigma(\phi)+Ux+Qu+s\epsilon$ (122)
$\displaystyle\tau_{Q}\dot{Q}$ $\displaystyle=-s\epsilon
u_{\mathrm{hp}}^{\top}-\gamma_{Q}Q$ (123)

Here, $\gamma_{Q}$ is a small weight-decay term, $\epsilon$ is noise from an
Ornstein-Uhlenbeck process, i.e. exponentially filtered white noise with
standard deviation $s$ and timeconstant $\tau_{\epsilon}$. As post-synaptic
plasticity signal, we use a high-pass filtered version $u_{\mathrm{hp}}$ of
the output control $u$, to extract the noise fluctuations. We model the
plasticity with continuous dynamics instead of a single update at equilibrium,
as the noise $\epsilon$ changes continuously. Building upon the work of
Meulemans et al. [30, 51], we show in Proposition S10 that the feedback
learning dynamics of Eq. 123 drives the feedback weights to a setting that
satisfies the column space condition of Eq. 120, when trained on a single data
point.

######  Proposition S10.

Assuming (i) a separation of timescales
$\tau\ll\tau_{\epsilon}\ll\tau_{u}\ll\tau_{Q}$, (ii) a small noise variance
$s^{2}\ll\lVert\phi\rVert^{2}$, (iii) a high-pass filtered control signal
$u_{\mathrm{hp}}=u-u_{\mathrm{ss}}$ with $u_{\mathrm{ss}}$ the steady state of
the first moment of the output controller trajectory $u$, and (iv) stable
dynamics, then, for a single fixed input $x$, the feedback learning dynamics
of Eq. 123 let the first moment $\mathbb{E}[Q]$ converge towards a setting
satisfying the column space condition of Theorem 2:

$\displaystyle\mathbb{E}\left[Q_{\mathrm{ss}}^{\top}\right]\propto\frac{\partial^{2}L}{\partial
y^{2}}(D\phi_{\mathrm{ss}})D(\mathrm{Id}-W\sigma^{\prime}(\phi_{\mathrm{ss}}))^{-1}$
(124)

with $\phi_{\mathrm{ss}}$ the steady state of the first moment of the state
trajectory $\phi$, and $\sigma^{\prime}(\phi)=\partial_{\phi}\sigma(\phi)$.

###### Proof.

The general idea behind the feedback learning rule of Eq. 123 is to correlate
the noise fluctuations of the output controller $u$ with the noise $\epsilon$
inside the neurons that caused these fluctuations. As the injected noise
propagates through the whole network, and the neural dynamics $\dot{\phi}$
evolve on a faster timescale compared to the noise $\epsilon$, the noise
fluctuations in the output contain useful information on the network
transformation $D(\mathrm{Id}-W\sigma^{\prime})^{-1}$ which can be extracted
by correlating these fluctuations with the noise $\epsilon$. We refer the
interested reader to Appendix C of Meulemans et al. [30] for a detailed
discussion of this type of feedback learning rule for feedforward neural
networks.

First, we define the neural fluctuations
$\tilde{\phi}:=\phi-\phi_{\mathrm{ss}}$ and output controller fluctuations
$\tilde{u}=u-u_{\mathrm{ss}}$, with $\phi_{\mathrm{ss}}$ and $u$ the steady
state of the first moment of $\phi$ and $u$, respectively. Using the
assumption that the noise variance $s^{2}$ is much smaller than the neural
activity $\lVert\phi\rVert^{2}$ and that we have stable (hence contracting)
dynamics, we perform a first-order Taylor approximation around
$\phi_{\mathrm{ss}}$ and $u$, which for simplicity we assume to be exact,
leading to the following dynamics:

$\displaystyle\tau\dot{\tilde{\phi}}$
$\displaystyle=-\tilde{\phi}+W\sigma^{\prime}(\phi_{\mathrm{ss}})\tilde{\phi}+Q\tilde{u}+\epsilon,$
(125) $\displaystyle\tau_{u}\dot{\tilde{u}}$
$\displaystyle=\frac{\partial^{2}L}{\partial
y^{2}}(D\phi_{\mathrm{ss}})D\tilde{\phi}-\alpha\tilde{u}$ (126)

where we used that
$-\phi_{\mathrm{ss}}+W\sigma(\phi_{\mathrm{ss}})+Ux+Qu_{\mathrm{ss}}=0$ and
$-\partial_{y}L(D\phi_{\mathrm{ss}})-\alpha u_{\mathrm{ss}}=0$ by definition.
To compress notation, we introduce the variable $\nu$ as a concatenation of
$\tilde{\phi}$ and $\tilde{u}$.

$\begin{split}\dot{\nu}&=-A\nu+B\epsilon\\\
A&:=\begin{bmatrix}\frac{1}{\tau}(\mathrm{Id}-W\sigma^{\prime}(\phi_{\mathrm{ss}}))&-\frac{1}{\tau}Q\\\
\frac{1}{\tau_{u}}\partial^{2}_{y}L(D\phi_{\mathrm{ss}})D&\frac{\alpha}{\tau_{u}}\mathrm{Id}\end{bmatrix}~{}~{}~{}B:=\begin{bmatrix}\frac{s}{\tau}\mathrm{Id}\\\
0\end{bmatrix}\end{split}$ (127)

We solve this linear time-invariant stochastic differential equation using the
method of variation of constants sarkka_applied_2019, while assuming that
$\nu=0$ at time $t_{0}$, i.e. that the dynamics already have converged to the
steady state at $t_{0}$.

$\nu(t)=\int_{t_{0}}^{t}\exp\left(-A(t-t^{\prime})\right)B\epsilon(t^{\prime})\mathrm{d}t^{\prime}.$
(128)

Now we turn to the first moment of the feedback weight learning (123), using
$\tilde{u}=C\nu$ with $C=[0~{}\mathrm{Id}]$:

$\mathbb{E}\left[\tilde{u}(t)\epsilon(t)^{\top}\right]=C\mathbb{E}\left[\nu(t)\epsilon(t)^{\top}\right]=\frac{1}{2\tau_{\epsilon}}C\int_{t_{0}}^{t}\exp\left(-\left(A+\frac{1}{\tau_{\epsilon}}\mathrm{Id}\right)(t-t^{\prime})\right)B\mathrm{d}t^{\prime}$
(129)

for which we used that
$\mathbb{E}\left[\epsilon(t)\epsilon(t^{\prime})^{\top}\right]=\frac{1}{2\tau_{\epsilon}}\exp(\frac{1}{\tau_{\epsilon}}|t-t^{\prime}|)$
for an Ornstein-Uhlenbeck process with time constant $\tau_{\epsilon}$. Now
assuming that $t\gg t_{0}$ and that we have stable dynamics (i.e. $A$ has
strictly positive eigenvalues), we can solve the integral of the matrix
exponential:

$\mathbb{E}\left[\tilde{u}(t)\epsilon(t)^{\top}\right]=\frac{1}{2\tau_{\epsilon}}C\left(A+\frac{1}{\tau_{\epsilon}}\mathrm{Id}\right)^{-1}B$
(130)

We solve the inverse of this $2\times 2$ block matrix analytically (following
e.g. lu_inverses_2002), leading to

$\begin{split}\mathbb{E}\left[\tilde{u}(t)\epsilon(t)^{\top}\right]&=\frac{1}{2\tau_{\epsilon}}\frac{s}{\tau}\left[-\Delta^{-1}\frac{1}{\tau_{u}}\frac{\partial^{2}L}{\partial
y^{2}}(D\phi_{\mathrm{ss}})D\tau(\mathrm{Id}-W\sigma^{\prime}(\phi_{\mathrm{ss}}))^{-1}\right]\\\
\Delta&:=\left(\frac{\alpha}{\tau_{u}}+\frac{1}{\tau_{\epsilon}}\right)\mathrm{Id}+\frac{1}{\tau_{u}}\frac{\partial^{2}L}{\partial
y^{2}}(D\phi_{\mathrm{ss}})D\left(\frac{1}{\tau}(\mathrm{Id}-W\sigma^{\prime}(\phi_{\mathrm{ss}}))\right)^{-1}\frac{1}{\tau}Q\end{split}$
(131)

Using the separation of timescales $\tau\ll\tau_{\epsilon}\ll\tau_{u}$, this
simplifies to

$\begin{split}\mathbb{E}\left[\tilde{u}(t)\epsilon(t)^{\top}\right]&\approx-\frac{s}{2\tau_{u}}\frac{\partial^{2}L}{\partial
y^{2}}(D\phi_{\mathrm{ss}})D(\mathrm{Id}-W\sigma^{\prime}(\phi_{\mathrm{ss}}))^{-1}\end{split}$
(132)

Remembering that $\tilde{u}=u_{\mathrm{hp}}$, we get the following plasticity
dynamics using Eq. 123:

$\tau_{Q}\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}\left[Q(t)^{\top}\right]:=\tau_{Q}\frac{\mathrm{d}}{\mathrm{d}t}\bar{Q}(t)=\frac{s^{2}}{2\tau_{u}}\frac{\partial^{2}L}{\partial
y^{2}}(D\phi_{\mathrm{ss}})D(\mathrm{Id}-W\sigma^{\prime}(\phi_{\mathrm{ss}}))^{-1}-\gamma_{Q}\bar{Q}(t)$
(133)

with $\bar{Q}(t)$ the first moment of $Q(t)$. For $\gamma>0$, these dynamics
are stable and lead to

$\displaystyle\mathbb{E}\left[Q_{\mathrm{ss}}^{\top}\right]=\frac{s^{2}}{2\tau_{u}\gamma_{Q}}\frac{\partial^{2}L}{\partial
y^{2}}(D\phi_{\mathrm{ss}})D(\mathrm{Id}-W\sigma^{\prime}(\phi_{\mathrm{ss}}))^{-1}$
(134)

thereby concluding the proof. ∎

Isolating the postsynaptic noise $\epsilon$ with a multicompartment neuron
model. The feedback learning rule of Eq. 123 and the corresponding Proposition
S10 represent the feedback learning idea in its simplest form. However, for
this learning rule, the feedback synapses need some mechanism to isolate the
postsynaptic noise $\epsilon$ from the (noisy) neural activity. Meulemans et
al. [30, 51] solve this by considering a multicompartment model of the neuron,
inspired by recent dendritic compartment models of the cortical pyramidal
neuron [17, 41, 14]. Here, a feedback (apical) compartment integrates the
feedback input $Qu$, a basal compartment integrates the recurrent and
feedforward input $W\sigma(\phi)+Ux$ and a central (somatic) compartment
combines the two other compartments together and transmits the output firing
rate $\sigma(\phi)$ to the other neurons. In this model, we can assume that a
part of the noise $\epsilon$ enters through the feedback compartment, and
change the feedback learning rule to

$\displaystyle\tau_{Q}\dot{Q}$
$\displaystyle=-\big{(}Qu+s\epsilon^{\mathrm{fb}}\big{)}u_{\mathrm{hp}}^{\top}-\gamma_{Q}Q,$
(135)

where $\big{(}Qu+s\epsilon^{\mathrm{fb}}\big{)}$ is the feedback compartment
activity, and hence locally available for the synaptic updates. As
$\mathbb{E}[Quu_{\mathrm{hp}}^{\top}]=Q\mathbb{E}[u_{\mathrm{hp}}u_{\mathrm{hp}}^{\top}]=Q\Sigma_{u}$
with $\Sigma_{u}$ the positive definite auto-correlation matrix of
$u_{\mathrm{hp}}$, the result of Proposition S10 can be adapted to
$\mathbb{E}\left[Q_{\mathrm{ss}}^{\top}\right]\propto
M\partial_{y}^{2}L(D\phi)_{\mathrm{ss}})D(\mathrm{Id}-W\sigma^{\prime}(\phi_{\mathrm{ss}}))^{-1}$,
with $M=(\Sigma_{u}+\gamma\mathrm{Id})^{-1}$ a positive definite matrix. This
new fixed point of $Q$ also satisfies the column space condition of Eq. 120.

Improving the parameter updates. When the column space condition is not
perfectly satisfied, as is the case with this direct linear feedback
controller when multiple data points are used, the empirical performance can
be improved by changing the parameter updates towards

$\displaystyle\Delta
W=\sigma^{\prime}(\phi_{*})Qu_{*}\sigma(\phi_{*})^{\top},\quad\Delta
U=\sigma^{\prime}(\phi_{*})Qu_{*}x^{\top}$ (136)

This learning rule uses the local derivative of the nonlinearity $\sigma$ as a
heuristic to prevent saturated neurons from updating their weights, which is
known to improve performance in feedforward networks [30, 51].

Alignment results. We empirically tested the controller with direct linear
feedback and trained feedback weights $Q$ on the MNIST digit recognition task.
Figure S3 shows that the updates of Eq. 136 are approximately aligned to the
gradients $-\nabla_{\theta}\cal{H}(\theta)$ on the least-control objective,
indicating that the feedback weight learning for $Q$ is successful and the
controller with direct linear feedback approximates the optimal control.

Figure S3: The linear feedback controller and inversion dynamics lead to
approximate optimal control. Cosine similarity between the true gradient
$-\nabla_{W}\mathcal{H}(\theta)$ and approximate updates $\Delta W$ of the
least-control principle for training an equilibrium RNN on MNIST: a controller
with trained linear feedback connections (LCP-LF) and a controller using
dynamic inversion with trained feedback weights via Kolen-Pollack learning
(LCP-DI KP).

#####  S6.3.2 Simulation details

For computational efficiency, we run fixed-point iterations for finding the
equilibrium of the controlled neural dynamics shown in Eq. 121, which in this
case is equivalent to simulating the dynamics with the forward Euler method
with a simulation stepsize equal to the timeconstant:

$\begin{split}\phi_{\mathrm{next}}&=W\sigma(\phi)+Ux+Qu\\\
u_{\mathrm{next}}&=(1-\alpha)u-\frac{\partial L}{\partial
y}(D\phi)\end{split}$ (137)

We use the fixed points $\phi_{*}$ and $u_{*}$ to update the parameters
$\theta=\\{W,U,b,\tilde{D},\tilde{b}\\}$ according to Eq. 117 with
$\psi_{*}=Qu_{*}$. For training the feedback weights, we simulate the
stochastic differential equations of Eq. 122 with the Euler-Maruyama method
sarkka_applied_2019, starting from the fixed points computed in Eq. 137:

$\begin{split}\epsilon_{\mathrm{next}}&=\epsilon+\frac{\Delta
t}{\tau_{\epsilon}}\left(-\epsilon+\frac{1}{\sqrt{\Delta
t}}\Delta\xi\right)\\\ \phi_{\mathrm{next}}&=\phi+\frac{\Delta
t}{\tau}\left(-\phi+W\sigma(\phi)+Ux+Qu+s\epsilon\right)\\\
u_{\mathrm{next}}&=u+\frac{\Delta t}{\tau_{u}}\left(-\frac{\partial
L}{\partial y}(D\phi)-\alpha u\right)\end{split}$ (138)

with $\Delta t$ the simulation stepsize and
$\Delta\xi\sim\mathcal{N}(0,\mathrm{Id})$ white Gaussian noise that is drawn
independently for each simulation step. The factor $\frac{1}{\sqrt{\Delta t}}$
in front of $\Delta\xi$ is due specific properties of Brownian motion in the
simulation of stochastic differential equations. For computational efficiency,
we accumulate the weight updates for $Q$ over the simulation interval, and
apply them at the end:

$\Delta
Q=-\frac{1}{\tau_{Q}}\Big{[}\frac{1}{s^{2}}\sum_{m}(Qu[m]+\epsilon[m])(u[m]-u_{*})^{\top}-\gamma
Q\Big{]}$ (139)

where $m$ is the simulation step and we use $u_{\mathrm{hp}}=u-u_{*}$, with
$u_{*}$ the fixed point of Eq. (137), and where we scale by the noise variance
$s^{2}$ to remove the dependence of the update magnitude on this
hyperparameter. For all parameter updates, we use the Adam optimizer.

####  S6.4 Using general-purpose dynamics for computing an optimal control

In the main text (Section 3.1.2) we discussed the application of our two
general-purpose optimal control methods to an equilibrium RNN. Here, we
provide an expanded discussion on the neural network architectures each method
yields. While we leave the development of more detailed biological models for
future work, we comment on how these controlled RNNs may be implemented in
cortical circuits, sketching a number of alternative circuit designs. We end
this section with providing additional simulation and experimental details.

#####  S6.4.1 Dendritic error compartments for dynamic inversion

In Section 3.1.2 we applied our general-purpose dynamic inversion method (8)
to steer equilibrium RNNs towards optimally-controlled states. We restate here
the resulting optimal control circuit (14) and discuss afterwards two possible
network implementations, where feedback error signals are represented in
dendritic compartments. Assuming without loss of generality that the output
controller is of the leaky-integral kind, the resulting network dynamics (14)
is:

$\begin{split}\dot{\phi}&=-\phi+W\sigma(\phi)+Ux+\psi,\\\
\dot{\psi}&=-\psi+\sigma^{\prime}(\phi)S\psi+D^{\top}u\\\ \dot{u}&=-\alpha
u-\nabla_{y}L(y).\end{split}$ (140)

where $y=D\phi$ is the output neuron state, a subset of the full neuron state
$\phi$ determined by a fixed selector matrix $D$. Figure S2.B shows an
instantiation of this circuit on a feedforward neural network.

In Eq. 140, the control variable $\psi$ follows its own dynamics $\dot{\psi}$,
which serves the purpose of inverting $\partial_{\phi}f(\phi,\theta)^{\top}$.
This can be seen by inspecting the controlled steady state:

$\psi_{*}=(\mathrm{Id}-\sigma^{\prime}(\phi_{*})S)^{-1}D^{\top}u_{*}.$ (141)

If the inverted factor appearing in the equation above would correspond to
$\partial_{\phi}f(\phi_{*},\theta)^{-\top}$, the network would be at an
optimally-controlled state, where the weight updates following from our theory
assume the Hebbian form

$\begin{split}\Delta W&=\psi_{*}\sigma(\phi_{*})^{\top}\\\ \Delta
U&=\psi_{*}x^{\top}.\end{split}$ (142)

For this to be the case, we need to satisfy the symmetry requirement
$S=W^{\top}$. This corresponds to a form of weight transport, assuming that
$S$ and $W$ are two distinct sets of synaptic connections.

We present now two alternative circuit implementations which allow maintaining
$S\approx W^{\top}$ using local learning rules, so that the equations
governing the change $\Delta W$, $\Delta U$ and $\Delta S$ of all synaptic
weights only depend on quantities presumed to be available pre- and post-
synaptically and at the same time. In both cases, we build upon previous
theories of error-driven cortical learning and assume that our control signal
$\psi$ corresponds to the apical dendritic activity of a cortical pyramidal
cell, and $\phi$ to its somatic activity. According to this view, apical
dendrites encode neuron-specific credit assignment signals $\psi_{*}$ which
are responsible for determining which neurons should potentiate or depress
their connections, and by how much. We therefore extend previous models of
apically-controlled synaptic plasticity [14, 41, 20, 23, 18], for which there
is some experimental evidence magee_synaptic_2020, to the case of recurrently-
connected networks.

Multiplexed codes. Interpreting our apical dendrite dynamics $\dot{\psi}$ in
literal terms appears to require dendrodendritic synapses which bypass the
soma, due to the apical-to-apical term $\psi\leftarrow S\psi$ occurring in
(140). While such connections are known to exist in the brain
shepherd_dendrodendritic_2009, they are rare. We thus highlight an alternative
implementation, which builds on the interesting proposal that pyramidal cells
might be capable of transmitting credit assignment signals alongside their
intended outputs, by means of a multiplexed code [23, 14]. The basic premise
of this theory is that certain synaptic connections are sensitive to high-
frequency bursts of spikes only, which can be triggered by increasing apical
activity, whereas other connections are tuned to filter out such events. Each
connection type would then differentially reflect either apical or somatic
activity, respectively; the required specific tuning to burst events or
otherwise can be achieved for example by resorting to different profiles of
short-term plasticity for each connection type. Here, we leverage on this idea
and propose to replace the dendrodendritic (direct apical-to-apical)
connections appearing in (140) by common (axodendritic) apical-targetting
synapses mediated through the soma (i.e., soma-to-apical $\psi\leftarrow\phi$
synapses), tuned to be sensitive to burst events. This design enables keeping
the circuit approximately symmetric ($S\approx W^{\top}$) through local burst-
dependent plasticity rules, $\Delta S=\sigma(\phi_{*})\psi_{*}^{\top}-\gamma
S$ and $\Delta W=\psi_{*}\sigma(\phi_{*})^{\top}-\gamma W$, a form of Kolen-
Pollack learning [54, 42].

Interestingly, the resulting network model differs in a conceptually important
– and experimentally testable – way from the multiplexing model originally
proposed by Payeur et al. [23]: in our model, apical control signals should
steer (thus, influence) the network state $\phi$ towards a controlled
equilibrium, and therefore should _not_ be filtered out from the somatic
voltage dynamics $\dot{\phi}$, while the model of Payeur et al. [23] fully
multiplexes credit assignment information and the actual computations carried
out by the network. Critically, our major functional requirement is that
apical-targetting connections be capable of selectively transmitting burst
events (e.g., through short-term facilitating plasticity), and we make no
further claims about the short-term plasticity of the remaining connections,
leaving open the possibility of using this mechanism for other computational
purposes unrelated to multiplexing.

Dendritic errors via lateral inhibition. Our network dynamics (140) is also
compatible with an alternative proposal for how apical dendrites may come to
encode credit assignment information, then used to steer synaptic plasticity
within a neuron [17, 39]. In that model, dendritic errors emerge thanks to an
auxiliary population of neurons, modeled after a particular class of
inhibitory interneurons which preferentially target apical dendrites. These
interneurons are set up such that they learn to reproduce the activity of
pyramidal neurons. Both pyramidal neurons and interneurons feedback to apical
dendrites, but with opposite sign. Sacramento et al. [17] showed that when the
interneurons perfectly track pyramidal neurons (sometimes referred to as a
‘tight’ inhibitory-excitatory balance hennequin_inhibitory_2017) and an output
error is generated, the remaining steady-state apical activity corresponds to
a backpropagated error.

We briefly sketch here a proof-of-concept mapping of our dynamics (140) onto
this model, while leaving more realistic implementations to future work. We
first replicate our recurrently-connected neurons with a population of
interneurons, whose state we denote by $\delta$, and let both types of neurons
feedback onto apical dendrites. Now, we follow the supplementary analysis of
ref. [17] and assume that both populations have exactly the same connectivity
and weights, a condition to ensure successful tracking and error propagation,
which can be reached with the local weight updates for the interneurons
introduced in ref. [17]. Moreover, for simplicity, here we assume our neurons
are linear and ignore factors involving the derivative of the transfer
function $\partial_{\phi}\sigma$, although it is possible to take such factors
into account with a more detailed construction. Taking
$\dot{\psi}=-\psi+S\phi-S\delta+D^{\top}u$ as our apical dendrite dynamics and
inspecting its steady-state yields
$\psi_{*}=(\mathrm{Id}-S)^{-1}D^{\top}u_{*}$, as the interneuron activity
perfectly tracks the pyramidal neuron activity $W\phi+Ux$, but does not
include the apical activity $\psi$. This is the optimally-controlled
equilibrium we sought. As with the multiplexed code described above,
approximate symmetry $S\approx W^{\top}$ can then be maintained with local
rules via Kolen-Pollack learning.

#####  S6.4.2 Error neurons for energy-based dynamics

In the main text and in Sections S2.5 and S4, we have seen that the optimal
control problem arising under our principle can be generically approached
using a (constrained) energy-minimizing dynamics. We restate here the simple
gradient-flow dynamics (15) presented in the main text:

$\begin{split}\dot{\phi}&=-\psi+\sigma^{\prime}(\phi)S\psi+D^{\top}u,\\\
\psi&=\phi-W\sigma(\phi)-Ux,\end{split}$ (143)

where $u$ can be any output controller which enforces the constraint that the
loss $L$ is minimized with respect to the output neural activity $y$. The
circuit corresponding to those dynamics is shown on Figure S2.C. We now
briefly discuss how this system can be implemented using a predictive coding
architecture featuring two distinct subpopulations of neurons, generalizing
the model of Whittington and Bogacz [36] to the RNN setting. In this
architecture, $\psi$ is interpreted as the activity of a population of
prediction error neurons, while $\phi$ corresponds to the activity of
prediction neurons. Error neurons $\psi$ act as comparators, which measure
deviations of the actual state of prediction neurons $\phi$ from the expected
state $\mu:=W\sigma(\phi)-Ux$; in this model, such deviations dynamically
adjust the predictions $\phi$ towards an optimally-controlled steady-state in
which the output loss is minimal (c.f. Section S4). The critical architectural
difference to the models of Whittington and Bogacz [36] and Rao and Ballard
[56] is that in our case the expected state $\mu$ does not necessarily stem
from different network layers. A hierarchically-structured connectivity matrix
$W$ is a special case of our model; we allow for top-down, lateral and bottom-
up contributions to the expected states $\mu$.

#####  S6.4.3 Training deep equilibrium models with the least-control
principle.

We use deep equilibrium models that consist out of three main modules: (i) an
encoder $g(x,\theta)$ that maps the input $x$ to the implicit layer $\phi$,
(ii) an implicit layer defined by $f(\phi,\theta)=0$ where $g(x,\theta)$ is
incorporated inside $f(\phi,\theta)$, and (iii) a decoder $h(\phi,\theta)$
that maps the implicit layer towards the output $y$.

When $h(\phi,\theta)$ is independent of $\theta$, we can use the inversion
dynamics of Eq. 8 or the energy-based dynamics of Eq. 10 to compute an optimal
control $\psi_{\scalebox{0.42}{$\bigstar$}}$ and use the first-order parameter
updates of Theorem 1. To harvest the full power of deep equilibrium models, it
is beneficial to train the decoder $h(\phi,\theta)$ as well. However, this
poses the problem that the loss $L(h(\phi,\theta))$ is now dependent on
$\theta$ and hence Theorem 1 cannot be applied out-of-the-box. There are two
general approaches to solve this challenge: (i) similar to the vanilla
equilibrium RNNs, we can extend the state $\phi$ with an extra set of output
neurons to absorb the $\theta$-dependent decoder $h(\phi,\theta)$ inside of
the implicit layer; and (ii) we can use the extended version of the least-
control principle derived in Section S2.2 that can handle $\theta$-dependent
losses, to derive the parameter updates.

Augmenting the implicit layer. We can use an augmented state $\bar{\phi}$ that
concatenates the system states $\phi$ with an extra set of output neurons
$\tilde{\phi}$. Next we augment the implicit layer $f(\phi,\theta)$ towards

$\begin{split}\bar{f}(\bar{\phi},\theta)=\begin{bmatrix}f(\phi,\theta)\\\
h(\phi,\theta)-\tilde{\phi}\end{bmatrix}.\end{split}$ (144)

Now, we can use a fixed linear decoder $D=[0~{}\mathrm{Id}]$ to select the
output states $\tilde{\phi}$ for the loss $L(D\bar{\phi})$, which now is
independent of $\theta$. This augmented system is then controlled by
$\bar{\psi}$ which has the following dynamics when using the dynamic inversion
algorithm of Eq. 8:

$\displaystyle\dot{\bar{\psi}}=\frac{\partial\bar{f}}{\partial\bar{\phi}}(\bar{\phi},\theta)^{\top}\bar{\psi}+D^{\top}u$
(145)

Combining this with the system dynamics of Eq. 144 leads to the following set
of equilibrium conditions

$\begin{split}0&=f(\phi_{*},\theta)+\psi_{*}\\\
0&=-\tilde{\phi}_{*}+h(\phi_{*},\theta)+\tilde{\psi}_{*}\\\ 0&=\frac{\partial
f}{\partial\phi}(\phi_{*},\theta)^{\top}\psi_{*}+\frac{\partial
h}{\partial\phi}(\phi_{*},\theta)^{\top}\tilde{\psi}_{*}\\\
0&=-\tilde{\psi}_{*}+u_{*}\\\ 0&=\alpha u_{*}+\frac{\partial L}{\partial
y}(\phi_{*})\end{split}$ (146)

with $\tilde{\psi}$ the control applied to the augmented output neurons
$\tilde{\phi}$. Grouping the parameters $\theta$ into the encoder parameters
$\theta_{g}$, the implicit module parameters $\theta_{f}$, and the decoder
parameters $\theta_{h}$, we get the following parameter gradients:

$\frac{\mathrm{d}}{\mathrm{d}\theta_{f}}\mathcal{H}(\theta)=\psi_{*}^{\top}\frac{\partial
f}{\partial\theta_{f}}(\phi_{*},\theta),\quad\frac{\mathrm{d}}{\mathrm{d}\theta_{g}}\mathcal{H}(\theta)=\psi_{*}^{\top}\frac{\partial
f}{\partial\theta_{g}}(\phi_{*},\theta),\quad\frac{\mathrm{d}}{\mathrm{d}\theta_{h}}\mathcal{H}(\theta)=\tilde{\psi}_{*}^{\top}\frac{\partial
h}{\partial\theta_{h}}(\phi_{*},\theta)$ (147)

Note that due to the general treatment, if the implicit module, encoder or
decoder consist of multiple computational layers, the above partial
derivatives do not result in a local update of presynaptic activity times
postsynaptic control signal, but require instead to backpropagate the credit
assignment signal (the control) through the computational layers within the
module. The updates however remain first-order updates which is important for
computational efficiency, and they remain local in time in contrast to
recurrent backpropagation (c.f. Section S3). When a design requirement is to
have local updates in space, e.g. for biological plausibility, we can augment
the system state $\phi$ with states representing each computational layer in
the modules, such that each layer receives a dedicated control, making the
resulting updates local in time and space.

Using the least-control principle for $\theta$-dependent losses. In Section
S2.2, we extend the least-control principle to $\theta$-dependent losses, and
among others derive a first-order update when the decoder $h(\phi,\theta)$ is
$\theta$-dependent, which we restate below.

$\begin{split}\frac{\mathrm{d}}{\mathrm{d}\theta}\mathcal{H}(\theta)&=-\psi_{*}^{\top}\frac{\partial
f}{\partial\theta}(\phi_{*},\theta)-u_{*}^{\top}\frac{\partial
h}{\partial\theta}(\phi_{*},\theta)\end{split}$ (148)

This first-order learning rule equals the gradient when the inversion dynamics
(8) or energy-based dynamics (10) are used for computing the optimal control
$\psi_{\scalebox{0.42}{$\bigstar$}}$, in the limit of $\alpha\to 0$.

As the strategy of augmenting the implicit layer to have a loss independent of
$\theta$ only uses the theory explained in the main manuscript, we use this
strategy for our experiments.

#####  S6.4.4 Simulation details

Equilibrium RNNs. We use the following fixed point iterations for finding the
equilibrium of the inversion dynamics of Eq. 14:

$\begin{split}\phi_{\mathrm{next}}&=W\sigma(\phi)+Ux+\psi\\\
\psi_{\mathrm{next}}&=\sigma^{\prime}(\phi)S\phi+D^{\top}u\\\
u_{\mathrm{next}}&=(1-\alpha)u-\frac{\partial L}{\partial
y}(D\phi)^{\top}\end{split}$ (149)

The fixpoint iterations terminates for a given input when either the relative
change in the norm of the dynamic states $\Phi=\\{\phi\\}$ or
$\Phi=\\{\phi,\psi,u\\}$ is less than some threshold $\epsilon$, i.e.
$\frac{\|\Phi_{t}-\Phi_{t+1}\|^{2}}{\|\Phi_{t}\|\|\Phi_{t+1}\|}\leq\epsilon$,
or when some iteration budget max_steps is exhausted.

When we learn the feedback weights $S$ with Kolen-Pollack learning (LCP-DI
KP), we use the following weight updates for $S$ and $W$:

$\Delta W=\psi_{*}\sigma(\phi)^{\top}-\gamma W,~{}~{}~{}\Delta
S=\sigma(\phi)\psi_{*}^{\top}-\gamma S$ (150)

with $\gamma$ the weight decay parameter. When we do not train the feedback
weights $S$ (LCP-DI), we fix them towards the transpose of the recurrent
weights $S=W^{\top}$ and we use standard weight update for $W$ without weight
decay (12).

Deep equilibrium models. We generalize the above fixed point iterations to
more complex recurrent architectures, by using the abstract inversion dynamics
of Eq. 8 and the equilibrium equations of Eq. 146. Note that we directly
substitute $\tilde{\psi}_{*}=u_{*}$ in the fixed point equations, to eliminate
an equation.

$\begin{split}\phi_{\mathrm{next}}&=\phi+f(\phi,\theta)+\psi\\\
\tilde{\phi}_{\mathrm{next}}&=h(\phi,\theta)+u\\\
\psi_{\mathrm{next}}&=\psi+\frac{\partial
f}{\partial\phi}(\phi,\theta)^{\top}\psi+\frac{\partial
h}{\partial\phi}(\phi,\theta)^{\top}u\\\
u_{\mathrm{next}}&=(1-\alpha)u-\frac{\partial L}{\partial
y}(\tilde{\phi})^{\top}\end{split}$ (151)

We use the first-order updates discussed in Eq. 147 for training the DEQs,
with $\tilde{\psi}_{*}=u_{*}$.

####  S6.5 Additional experimental details

Here, we provide additional details on the experiments introduced in Section
3.1, including the architecture of the networks, the used hyperparameters,
possible data preprocessing and the used optimization algorithms.

#####  S6.5.1 Training feedforward networks and equilibrium RNNs on MNIST

We evaluate the least-control principle combined with the various feedback
circuits introduced in Section 3.1 on the MNIST digit classification task
[52], and compare their performance to recurrent backpropagation (RBP).

Data preprocessing and augmentation. No data preprocessing or augmentation was
employed in our experiment.

Feedforward network architecture. We use two hidden layers of 256 neurons to
have the same amount of parameters as the equilibrium RNN architecture used
for MNIST. In the hidden layers, we use $\tanh$ nonlinearities and on the
output layer a softmax nonlinearity.

RNN architecture. We use 256 recurrent units in the RNN used for the results
of table 1. We include a trained decoder $\tilde{D}$ inside the recurrent
layer, as detailed in Section S6.1 and Eq. 109. For the results shown in Fig.
2, we run the equilibrium RNN on MNIST with $[23,32,45,64,90,128,181,256]$
recurrent units. We use the $\tanh$ nonlinearity for all recurrent units, and
include bias terms according to Eq. 110.

Optimization. For the classification benchmark, we train the models for 30 of
the MNIST train set, with a minibatch size of 64 and with the Adam optimizer
kingma_adam_2015. When we learn the feedback weights $S$ with Kolen-Pollack
learning, we apply the weight decay after every gradient update, i.e. we apply
it outside of the optimizer.

For training the feedback weights $Q$ of the linear feedback controller, we
run the stochastic dynamics in Eq. 138 for max_steps timesteps, and we update
$Q$ following Eq. 139 by averaging $\Delta Q$ over the minibatch. For the
simulation of the stochastic dynamics, we use $\Delta t=\tau^{\prime 2}$,
$\tau=\tau^{\prime 2}$, $\tau_{\epsilon}=\tau^{\prime}$ and $\tau_{u}=1$ for
some hyperparameter $\tau^{\prime}$. Table S2 contains details on the used and
scanned hyperparameters.

Energy-based dynamics LCP-EBD models were trained by optimizing the energy for
max_steps steps using the Adam optimizer [88] with learning rate ebd_inner_lr.
Due to computational efficiency reasons, we initialized $\phi$ at the free-
dynamic equilibrium before the energy optimization for the RNN experiments.
For the feedforward networks, we initialized $\phi$ at zero. As initializing
the network to its free equilibrium results effectively in a two-phase
algorithm, we investigated whether the same fixed points are reached when
initializing at the free equilibrium or at a zero. We observed that both
initializations converge to the same fixed point, but that the zero
initialization needs many more iterations to reach the fixed point, compared
to the combined free-phase iterations and controlled dynamics iterations
starting from the free equilibrium.

For the overparametrization experiments, we use the same hyperparameters
except for a learning rate of $0.001$ for both RBP and LCP-DI, with a cosine
annealing scheduler loshchilov_sgdr_2017 annealing the learning rate down to
$0.0001$ over a total training time of 50 epochs.

For all RNN experiments, we use gradient clipping whereby the norm of the
gradient, computed over all gradients concatenated into a single vector, is
normalized to some value (chosen to be 10) whenever above it. The gradient
clipping was particularly crucial for a stable training of the RBP model.

Additional results. In table S3 we show additional experimental results for
the MNIST experiments. Note that the hyperparameter was chosen to maximise the
test set accuracy.

Robustness to $\alpha$. While our theory characterizes the strong control
limit $\alpha\rightarrow 0$, some interesting properties still hold outside
this regime. In Proposition S9, we showed that under a strict alignment
condition on the feedback mapping $Q$ (Eq. 76), the controlled steady state is
a local minimizer of the augmented energy
$F(\phi,\theta,\alpha)=||f(\phi,\theta)||^{2}+\alpha^{-1}L(\phi)$. This strict
alignment condition holds for energy-based (Eq. 10) and inversion (Eq. 8)
dynamics. This has two consequences. First, the update prescribed by Theorem 1
follows the gradient of some surrogate objective function (cf. the end of
Section S2.6) which is closely related to the least-control objective (cf.
Remark 1 in Section S2.5). Additionally, this surrogate objective (for
$\alpha\neq 0$) is still a meaningful one for the original learning problem:
similar results to Propositions 3 and 4 hold for this objective (c.f. Section
S2). However, the feedback mapping flexibility of Theorem 2 does not hold
anymore outside the strong control limit $\alpha\rightarrow 0$, as the column
space condition of Eq. 7 is now replaced by the strict alignment condition of
Eq. 76.

To corroborate these theoretical insights and to verify whether the least-
control principle exhibits robust performance for a wide range of values of
$\alpha$, we repeat the MNIST experiments for LCP-DI with a feedforward neural
network for varying values of the controller leakage $\alpha$. Note that for
$\alpha=0$, a softmax output layer combined with one-hot targets would lead to
an infinite control, as a softmax needs an infinite input to produce a one-hot
output. This can be remedied by either increasing $\alpha$ or by using soft
targets, i.e. taking $1-\epsilon$ as target for the correct class, and
$\epsilon/(N-1)$ for the other $N-1$ classes. As we want to include $\alpha=0$
in this line of experiments, we resort here to soft targets with
$\epsilon=0.01$, whereas for the experiments in the main text we increase the
controller leakage to $\alpha=0.1$ (c.f. Table S2). Table S1 shows that the
performance of LCP-DI is robust to a wide range of values for $\alpha$.

Table S1: Test set accuracy on MNIST for a feedforward network (2x256 neurons) trained by LCP-DI with different values for the controller leakage $\alpha$. We used the hyperparameters mentioned in Table S2, combined with forward Euler integration of the dynamical equations with timestep $\texttt{dt}=0.01$.  | Test accuracy
---|---
$\alpha=0$ | 98.09 %
$\alpha=0.01$ | 98.04 %
$\alpha=0.1$ | 98.10 %
$\alpha=1$ | 97.99 %
Table S2: Hyperparameter search space for the MNIST Feed forward and RNN experiment. A grid search with one seed for each configuration was done to find the best parameters, which are marked in bold. The same set of hyperparameters for both RNN and feedforward models yielded the best performance, except for the LCP-EBD model where the RNN model needed shorter inner optimization steps due the free-phase initialization. Hyperparameter | RBP | LCP-DI | LCP-DI (KP) | LCP-LF | LCP-EBD
---|---|---|---|---|---
num_epochs | 30 | 30 | 30 | 30 | 30
batch_size | 64 | 64 | 64 | 64 | 64
$\alpha$ | - | 0.1 | 0.1 | 0.1 | 0.1
lr | $\\{3.10^{-4},\mathbf{10^{-3}},3.10^{-3}\\}$ | $\\{3.10^{-4},\mathbf{10^{-3}},3.10^{-3}\\}$ | $\\{3.10^{-4},\mathbf{10^{-3}},3.10^{-3}\\}$ | $\\{3.10^{-4},\mathbf{10^{-3}},3.10^{-3}\\}$ | $\\{3.10^{-4},\mathbf{10^{-3}},3.10^{-3}\\}$
optimizer | ADAM | ADAM | ADAM | ADAM | ADAM
scheduler | $\\{none,\mathbf{cos}\\}$ | $\\{none,\mathbf{cos}\\}$ | $\\{none,\mathbf{cos}\\}$ | $\\{none,\mathbf{cos}\\}$ | $\\{none,\mathbf{cos}\\}$
max_steps | $\\{50,\mathbf{200},800\\}$ | $\\{50,200,\mathbf{800}\\}$ | $\\{200,\mathbf{800}\\}$ | $\\{200,\mathbf{800}\\}$ | 800 (FF) / 200 (RNN)
$\epsilon$ | $\\{10^{-5},\mathbf{10^{-4}},10^{-3}\\}$ | $\\{\mathbf{10^{-6}},10^{-5},10^{-4}\\}$ | $\\{\mathbf{10^{-6}},10^{-5},10^{-4}\\}$ | $\\{\mathbf{10^{-6}},10^{-5},10^{-4}\\}$ | -
$\gamma$ | - | - | $\\{\mathbf{10^{-6}},10^{-5},10^{-4},10^{-3}\\}$ | - | -
$s$ | - | - | - | $\\{\mathbf{0.01},0.1\\}$ | -
$\tau^{\prime}$ | - | - | - | $\\{0.1,\mathbf{0.2}\\}$ | -
$\tau_{Q}$ | - | - | - | $\\{100,1000,10000,\mathbf{100000}\\}$ | -
$\gamma_{Q}$ | - | - | - | $\\{10^{-4},10^{-3},\mathbf{10^{-2}},0.1\\}$ | -
ebd_inner_optimizer | - | - | - | - | ADAM
ebd_inner_lr | - | - | - | - | 0.01 (FF) / 0.001 (RNN)
| Train set | Test set
---|---|---
NLL | Accuracy | NLL | Accuracy
FF-LCP-LF | 0.00087 ± 0.00084 | 99.98 ± 0.03 | 0.1910 ± 0.0055 | 97.73 ± 0.07
FF-LCP | 0.00004 ± 0.00002 | 100.00 ± 0.00 | 0.1869 ± 0.0154 | 98.11 ± 0.07
FF-LCP-KP | 0.00003 ± 0.00002 | 100.00 ± 0.00 | 0.2041 ± 0.0173 | 98.14 ± 0.09
FF-LCP-EBD | 0.02645 ± 0.00183 | 99.90 ± 0.00 | 0.0848 ± 0.0020 | 98.00 ± 0.03
FF-BP | 0.00002 ± 0.00001 | 100.00 ± 0.00 | 0.0879 ± 0.0027 | 98.29 ± 0.14
RNN-LCP-LF | 0.0021 ± 0.0019 | 99.96 ± 0.04 | 0.1917 ± 0.0068 | 97.70 ± 0.11
RNN-LCP | 0.0121 ± 0.0013 | 99.80 ± 0.02 | 0.1630 ± 0.0007 | 97.58 ± 0.12
RNN-LCP-KP | 0.0105 ± 0.0011 | 99.84 ± 0.04 | 0.1739 ± 0.0012 | 97.75 ± 0.11
RNN-LCP-EBD | 0.0103 ± 0.0026 | 99.74 ± 0.05 | 0.1324 ± 0.0079 | 97.60 ± 0.15
RNN-RBP | 0.0036 ± 0.0047 | 99.94 ± 0.09 | 0.0961 ± 0.0022 | 97.87 ± 0.19
Table S3: Classification accuracy ($\%$) and negative log likelihood (NLL) on
the MNIST train and test set, of a Feed forward (FF) and RNN model trained by
backpropagation (BP/RBP) and least control (LCP). The reported mean±std is
computed over 3 seeds.

#####  S6.5.2 Training convolutional networks and deep equilibrium models on
CIFAR-10

Here, we train a convolutional network and Deep Equilibrium model and evaluate
their performances on the CIFAR-10 image classification benchmarks [60]. We
compare the performance of the models trained via either RBP or LCP.

Data preprocessing and augmentation. We normalize all images channel-wise by
subtracting the mean and dividing by the standard deviation, which are
computed on the training dataset. We employed no data augmentation in our
experiments.

Convolutional network architecture. We use a convolutional network with 3
convolutional layers of kernel 5-by-5 and channel sizes 96-128-256 followed by
two hidden fully connected layers of 2048 neurons. All layers are followed by
ReLU non linearity. We do not use batch normalization (for simplicity, we
dispense with normalization layers altogether). To keep the architecture as
simple as possible we do not use max-pooling units except after the first
convolutional layer, and resort for all other convolutional layers to simple
strided convolutions to implement downsampling.

Deep equilibrium model architecture. The model consists of 3 sub-modules that
are invoked sequentially: an encoder, an implicit module, and a decoder. The
encoder consists in a single convolutional layer with bias, followed by a
Batchnorm layer ioffe_batch_2015. The implicit module finds the fixed point of
the following equation:

$f(\phi,e,\theta)=-\phi+\mathrm{norm}(\mathrm{ReLU}(\phi+\mathrm{norm}(e+c_{2}(\mathrm{norm}(\mathrm{ReLU}(c_{1}(\phi)))))))$
(152)

where both $c_{1},c_{2}$ are convolutional layers, $\mathrm{norm}$ are group
normalization layers wu_group_2018, and $e$ is the output of the encoder
module $g(x,\theta)$. Finally, the decoder module contains a batchnorm layer,
followed by an average-pooling layer of kernel size 8-by-8, and a final linear
classification layer. All convolutional layers use 48 channels, a filter size
of 3-by-3, padding of 1, and are without biases unless specified otherwise.

Optimization. We train the model by finding the equilibrium of the controlled
dynamics using the fixed point iterations of Eq. 151 and we update the
parameters with the first order updates of Eq. 147. Note that the implicit
module contains multiple convolutional layers, hence the partial derivatives
in Eq. 147 backpropagate the credit assignment signal $\psi$ towards the
parameters of the different computational layers. We use automatic
differentiation for this. If local parameter updates in space are a design
requirement, one can augment the system state $\phi$ with additional states
for the output of each convolutional layer, resulting in dedicated control
signals for each computational layer, similarly to Section S6.4.3 where we
augment the state $\phi$ with output neurons to incorporate the decoder. Both
RBP and LCP-DI are trained on 75 epochs for the Deep Equilibrium model, while
the Convolutional network is trained on 30 epochs. For both models, we use the
stochastic gradient descent optimizer with a minibatch size of 64, combined
with an annealing learning rate using the cosine scheduler
loshchilov_sgdr_2017 which divides the learning rate by 10 over the training
epochs. We used a controller leakage of $\alpha=1$ to enable a faster
convergence to the controlled equilibrium and hence reduce the computational
requirements. Table S4 contains further details on the hyper parameters.

Additional results. In table S5 we show additional experimental results for
the CIFAR-10 experiments. Note that the hyperparameter was chosen to maximise
the test set accuracy. As can be seen, the training loss is not close to 0,
indicating that the model considered here is not overparametrized enough for
perfectly fitting the training set. Yet, LCP still manages to optimize for the
training loss, and crucially, in such a way that generalizes competitively.

Table S4: Hyperparameters used for the CIFAR-10 experiment, for the convolutional (CONV) and deep equilibrium model (DEQ) for backpropagation (-BP), recurrent backpropagation (-RBP) and LCP-DI training algorithms . Hyperparameter | CONV-BP | CONV-LCP-DI | DEQ-RBP | DEQ-LCP-DI
---|---|---|---|---
num_epochs | 30 | 30 | 75 | 75
batch_size | 64 | 64 | 64 | 64
$\alpha$ | - | 0.1 | - | 3
lr | 0.1 | 0.1 | 0.03 | 0.05
optimizer | SGD | SGD | SGD | SGD
scheduler | cos | cos | cos | cos
max_steps | - | 800 | 200 | 800
$\epsilon$ | - | $10^{-6}$ | $10^{-4}$ | $10^{-4}$
Gradient clipping | No | No | Yes | Yes
| Train set | Test set
---|---|---
NLL | Accuracy | NLL | Accuracy
DEQ-RBP | 0.2403±0.0027 | 92.50±0.14 | 0.6172±0.0097 | 80.14±0.20
DEQ-LCP | 0.3218±0.0062 | 89.55±0.17 | 0.5874±0.0077 | 80.26±0.17
Table S5: Classification accuracy ($\%$) and negative log likelihood (NLL) on
the CIFAR-10 train and test set, of a DEQ model trained by recurrent
backpropagation (RBP) and least control (LCP). The reported mean±std is
computed over 3 seeds.

#####  S6.5.3 Neural implicit representations

Inspired by the recent success of using implicit layers for implicit neural
representations [10], we investigate using the least-control principle in this
setup for learning. Implicit neural representations (INRs) are a quickly
emerging field in computer vision and beyond with the aim to represent any
signal through a continuous and differentiable neural network. Here we
strictly follow the experimental setup of Huang et al. [10] and leverage the
continuity of implicit neural representations to generalize on trained images.
In particular, given an Image $I$, we train the INR to represent $25\%$ of the
pixels (we simply jump over every second pixel and every second row) while
testing the networks prediction on $25\%$ of unknown pixels and measuring the
peak signal-to-noise ratio (PSNR). This can be regarded as simple image
inpainting.

Architecture. Compared to Huang et al. [10], we slightly modify the implicit
layer resulting in the following free dynamics

$f_{\text{INR}}(\phi;x)=-\phi+\mathrm{ReLU}(W\phi)+Ux.$ (153)

We include a trained decoder matrix $\tilde{D}$ inside the implicit layer, as
discussed in Section S6.1. We use fixed point iterations on the inversion
dynamics (151) for finding an optimal control, and use the weight updates for
the equilibrium RNNs (117). Note that we did not include spectral
normalization or a learning rate scheduler to train our models as in Huang et
al. [10]. Despite these simplifications we obtain comparable results across
methods and models.

Hyperparameters. Note the small discrepancy between the INR models trained
with RBP and LCP-DI in Table S7. Interestingly, we had an easier time scanning
hyperparameters to obtain competitive performance and stable forward dynamics
when using LCP-DI compared to RBP, for which we had to include gradient
clipping and weight decay. See Table S7 for the hyperparameters found when
repeating the the grid search identical to the one seen in Table S2. Following
Huang et al. [10], we include an extra hyperparameter (Omega), which is a
scalar multiplier applied to all neural activity, and we scale the
initialization of the weights accordingly with $1/$Omega.

Additional results. Table S6 shows the full experimental results for learning
neural implicit representations on the natural images and text dataset
introduced by Huang et al. [10]. As we changed the network architecture
slightly, we include the SIREN and iSIREN results of Huang et al. [10] for
comparison. We see that the least-control principle achieves competitive
performance on this task, compared to RBP and the results of Huang et al.
[10].

Table S6: Peak signal-to-noise ratio (PSNR; in dB) for all models on the natural image and text generalization task [10]. The reported mean ± std is taken over the individual PSNR of the 16 images. SIREN and iSIREN results are taken from Huang et al. [10]. Note that while our INR model is trained without spectral normalization, we reach competitive results. Interestingly, LCP does perform on par with RBP while not relying on weight-decay, learning-rate scheduling or gradient clipping. | Natural | Text
---|---|---
1L-256D | 1L-512D | 1L-256D | 1L-512D
SIREN (input inj.) | 22.88 ± 3.0 | 24.52 ± 3.28 | 24.54 ± 2.19 | 25.69 ± 2.18
iSIREN | 24.28 ± 3.37 | 24.92 ± 3.58 | 26.06 ± 2.18 | 26.81 ± 2.09
INR-RBP | 24.72 ± 3.90 | 25.47 ± 4.16 | 25.34 ± 1.67 | 26.53 ± 2.40
INR-LCP | 24.25 ± 2.72 | 25.11 ± 2.90 | 26.71 ± 2.40 | 27.53 ± 2.00
Table S7: Hyperparameter used for the INR experiments. Hyperparameter | RBF-Natural | LCP-Natural | RBF-Text | LCP-Text
---|---|---|---|---
iteration per image | 500 | 500 | 500 | 500
$\alpha$ | - | 0 | - | 0
lr | 0.0001 | 0.0001 | 0.0001 | 0.0001
optimizer | Adam | Adam | Adam | Adam
max_steps | 50 | 50 | 50 | 50
$\epsilon$ | $10^{-5}$ | $10^{-5}$ | $10^{-5}$ | $10^{-5}$
weight decay | 0.1 | 0.1 | 0.1 | $0.1$
Omega | 15 | 35 | 35 | 50
Gradient clipping | Yes | No | Yes | No

####  S6.6 Computational cost

Here, we include a computational comparison between LCP-DI and RBP on the
MNIST experiments for equilibrium RNNs. The table below shows the amount of
iterations (phase length), the computation time for one single iteration
(silicon time; obtained by using the NVIDIA/PyProf profiler offline), and the
total training time for both algorithms. Note that RBP has two phases, whereas
LCP only one. We see that LCP requires more iterations compared to the
combined two phases of RBP. Furthermore, RBP requires significantly fewer
computations per iteration, as LCP needs to compute both the system states and
control states every iteration, resulting in more matrix-vector products.
Consequently, LCP requires significantly more compute compared to RBP on this
line of experiments, when simulating on standard digital hardware. Hence,
despite its desirable single-phase property, the added computational cost per
iteration results in a significantly higher computational cost of LCP. We
observed that increasing the controller leakage $\alpha$ decreases the number
of iterations required for reaching the controlled equilibrium significantly,
and hence also decreases the computational cost.

Table S8: Computational cost | LCP-DI | RBP
---|---|---
Average 1st phase length | 483 | 77
Average silicon time per iteration (1nd phase) (ns) | 221634 | 28735
Average 2nd phase length | - | 72
Average silicon time per iteration (2nd phase) (ns) | - | 32479
Total training time (s) | 35280 | 1683

####  S6.7 Related work

Meulemans et al. [30] introduced an approach for training feedforward neural
networks with local update rules by minimizing control. The mechanistic
implementation of this method can be considered as a specific instance of the
least-control principle applied to feedforward neural networks, where we use
direct linear feedback connectivity from an output controller to the network
to approximate an optimal control (c.f. Section S6.3). Importantly however,
the least-control principle is the first to put the idea of learning by
minimizing control on a strong and general theoretical foundation, as the
theory of Meulemans et al. [30] contains a problematic assumption.

When computing the gradient of the surrogate loss $\mathcal{H}$ w.r.t.
$\theta$ with the implicit function theorem, Meulemans et al. 2022 consider
the direct feedback weights $Q$ to be independent of $\theta$. However, after
taking an update, the feedback weights need to change to satisfy the column
space condition (Meulemans et al. 2022, Condition 1) again for the new
parameter setting $\theta$. Hence, the loss $\mathcal{H}=\|Qu\|^{2}$ is
affected by this change in Q, which was not incorporated in the gradient
calculation, making Strong-DFC only a heuristically-derived method. We were
not able to incorporate a $\theta$ -and $\phi$ dependent feedback mapping $Q$
into the theoretical framework following the methods of Meulemans et al. 2022,
using the implicit function theorem alone. We could only resolve these
fundamental issues by formulating least-control as an optimal control problem.

Furthermore, the least-control principle provides a way to learn general
dynamical systems that reach an equilibrium. These systems perform
computations both in space and time and cannot in general be reduced to a
feedforward computational graph. Consequently, the least-control principle
provides credit assignment both in space and time, significantly widening the
range of systems that can be learned compared to Meulemans et al. [30].

Finally, we introduce various different controller designs compatible with the
least-control principle, that move beyond the direct linear feedback
controller introduced by Meulemans et al. [30]. Importantly, the inversion
dynamics of Eq. 8 and the energy-based dynamics of Eq. 10 satisfy the column
space condition of Theorem 2 exactly and thereby lead to gradient-following
updates, in contrast to the direct linear feedback that cannot satisfy the
column space condition perfectly and hence results in only approximate
gradient updates.

###  Appendix S7 The least-control principle for meta-learning

In Section 3.2, we consider a meta-learning problem in which the objective is
to improve the performance of a learning algorithm as new tasks are
encountered. To do so, we equip a neural network with fast weights $\phi$ that
are learned on each task, and slow parameters $\theta=\\{\omega,\lambda\\}$
that parametrize the learning algorithm. The fast weights $\phi$ quickly learn
how to solve a task $\tau$ and $\lambda$ determines how strongly $\phi$ is
pulled towards the slow weights $\omega$, which consolidate the knowledge that
has been previously acquired. The capabilities of the learning algorithm are
then measured by evaluating $L^{\mathrm{eval}}_{\tau}(\phi^{*}_{\tau})$, the
performance of the learned neural network on data from task $\tau$ that was
held out during learning. Meta-learning finally corresponds to improving the
learning performance on many tasks by adjusting the meta-parameters
$\theta=\\{\omega,\lambda\\}$. This can be formalized through the following
bilevel optimization problem [35, 68]:

$\min_{\theta}\,\mathbb{E}_{\tau}\\!\left[L^{\mathrm{eval}}_{\tau}(\phi^{*}_{\tau})\right]\quad\mathrm{s.t.}\enspace\phi_{\tau}^{*}=\operatorname*{arg\,min}_{\phi}L^{\mathrm{learn}}_{\tau}(\phi)+\sum_{i}\frac{\lambda_{i}}{2}(\phi_{i}-\omega_{i})^{2}.$
(154)

Note that compared to Eq. 16 in the main text, we now consider the more
general setting in which each synapse $i$ has its own attraction strength
$\lambda_{i}$ and meta-learn them.

We now derive the update rules prescribed by the least-control principle,
review the meta-learning algorithms we are comparing against and detail our
experimental setup.

####  S7.1 Derivation of the meta-parameter updates

The first step of the derivation is to replace the $\operatorname*{arg\,min}$
in Eq. 154 by the associated stationarity condition:

$\nabla_{\phi}L^{\mathrm{learn}}_{\tau}(\phi^{*}_{\tau})+\lambda(\phi^{*}_{\tau}-\omega)=0.$
(155)

This equation can be obtained from a dynamical perspective: let us assume that
our learning algorithm is gradient descent, so that

$\dot{\phi}=f(\phi,\theta)=-\nabla_{\phi}L^{\mathrm{learn}}_{\tau}(\phi)-\lambda(\phi-\omega).$
(156)

The stationarity conditions therefore correspond to
$f(\phi^{*}_{\tau},\theta)=0$ and we can apply the least control principle.

To do so, we run the controlled dynamics

$\dot{\phi}=-\nabla_{\phi}L^{\mathrm{learn}}_{\tau}(\phi)-\lambda(\phi-\omega)+u,\quad\mathrm{and}\quad\dot{u}=-\alpha
u-\nabla_{\phi}L_{\tau}^{\mathrm{eval}}(\phi),$ (157)

and note $\phi_{*}^{\tau}$ and $u_{*}^{\tau}$ the equilibrium. Note that this
corresponds to taking $Q=\mathrm{Id}$ and $h(\phi)=\phi$ in the general
feedback dynamics of Eq. 6. Theorem 2 therefore guarantees that
$\phi_{*}^{\tau}$ is an optimally controlled state and $u_{*}^{\tau}$ an
optimal control in the limit $\alpha\rightarrow 0$, as long as
$\partial_{\phi}f(\phi_{*}^{\tau},\theta)$ and
$\partial_{y}^{2}L(h(\phi_{*}^{\tau}))$ are invertible, and the column space
condition

$\mathrm{col}\left[Q\right]=\mathrm{row}\left[\frac{\partial
h}{\partial\phi}(\phi_{*}^{\tau})\frac{\partial
f}{\partial\phi}(\phi_{*}^{\tau},\theta)^{-1}\right]$ (158)

is satisfied. The column space conditions here rewrites

$\mathrm{col}[\mathrm{Id}]=\mathrm{row}\left[\left(\frac{\partial^{2}L^{\mathrm{learn}}_{\tau}}{\partial\phi^{2}}(\phi_{*}^{\tau})+\lambda\mathrm{Id}\right)^{-1}\right]\\!,$
(159)

and is automatically satisfied through the assumption that the Jacobian of $f$
is invertible. To summarize, the controlled dynamics are optimal in the limit
$\alpha\rightarrow 0$, if
$\partial_{y}^{2}L^{\mathrm{learn}}_{\tau}(\phi^{\tau}_{*})+\lambda\mathrm{Id}$
and $\partial_{y}^{2}L^{\mathrm{eval}}_{\tau}(\phi^{\tau}_{*})$ are
invertible.

We now put ourselves in a regime in which the conditions detailed above are
satisfied, so that $\phi_{*}^{\tau}=\phi_{\scalebox{0.42}{$\bigstar$}}^{\tau}$
and $u_{*}^{\tau}=u_{\scalebox{0.42}{$\bigstar$}}^{\tau}$. Theorem 1 then
gives

$\Delta\theta=-\left(\frac{\mathrm{d}}{\mathrm{d}\theta}\mathcal{H}(\theta)\right)^{\top}=\frac{\partial
f}{\partial\theta}(\phi_{*}^{\tau},\theta)^{\top}u_{*}^{\tau}=-\frac{\partial
f}{\partial\theta}(\phi_{*}^{\tau},\theta)^{\top}f(\phi_{*}^{\tau},\theta),$
(160)

as $\psi=Qu$ is here equal to $u$. A simple calculation gives

$\begin{split}\frac{\partial
f}{\partial\omega}(\phi_{*}^{\tau},\theta)&=\mathrm{diag}(\lambda)\\\
\frac{\partial
f}{\partial\lambda}(\phi_{*}^{\tau},\theta)&=-\mathrm{diag}(\phi_{*}^{\tau}-\omega),\end{split}$
(161)

where $\mathrm{diag}(x)$ is a diagonal matrix whose diagonal values are the
elements of the vector $x$. It further implies that the updates of the
parameters are local, as

$\begin{split}\Delta\omega&=\lambda
u_{*}^{\tau}=\lambda\left(\nabla_{\phi}L_{\tau}^{\mathrm{learn}}(\phi_{*}^{\tau})+\lambda(\phi_{*}^{\tau}-\omega)\right)\\\
\Delta\lambda&=-(\phi_{*}^{\tau}-\omega)u_{*}^{\tau}=-(\phi_{*}^{\tau}-\omega)\left(\nabla_{\phi}L_{\tau}^{\mathrm{learn}}(\phi_{*}^{\tau})+\lambda(\phi_{*}^{\tau}-\omega)\right)\\!,\end{split}$
(162)

the multiplications being performed element wise.

####  S7.2 Comparison to existing meta-learning algorithms

We here review existing meta-learning rules that we compare against in Section
3.2. We separate those methods in two categories, the one relying on implicit
differentiation and the Reptile algorithm [76].

Implicit differentiation methods. The first set of methods we are reviewing
aims at directly minimizing the bilevel optimization problem (154) using the
implicit gradient

$\begin{split}\frac{\mathrm{d}}{\mathrm{d}\theta}L^{\mathrm{eval}}_{\tau}(\phi^{*}_{\tau})&=\frac{\partial
L^{\mathrm{eval}}_{\tau}}{\partial\phi}(\phi_{\tau}^{*})\left(\frac{\partial^{2}L_{\tau}^{\mathrm{train}}}{\partial\phi^{2}}(\phi_{\tau}^{*})+\mathrm{diag}(\lambda)\right)^{-1}\frac{\partial^{2}}{\partial\theta\partial\phi}\left[\sum_{i}\frac{\lambda_{i}}{2}(\phi^{*}_{\tau,i}-\omega_{i})^{2}\right]\\!.\end{split}$
(163)

In practice, this quantity is intractable because of the inverse of the
Hessian and has to be estimated. This lead to different algorithms. T1-T2 [75]
uses a first-order approximation of the implicit gradient by simply ignoring
the inverse Hessian term:

$\Delta\omega^{\mathrm{T1-T2}}_{i}=\lambda_{i}\frac{\partial
L^{\mathrm{eval}}_{\tau}}{\partial\phi_{i}}(\phi_{\tau}^{*}),\quad\mathrm{and}\quad\Delta\lambda_{i}^{\mathrm{T1-T2}}=-(\phi_{\tau,i}^{*}-\omega_{i})\frac{\partial
L^{\mathrm{eval}}_{\tau}}{\partial\phi_{i}}(\phi_{\tau}^{*}).$ (164)

Alternatively, a better estimate of the gradient can be obtained by
iteratively minimizing the quadratic form

$\delta\mapsto\delta\frac{1}{2}\left(\frac{\partial^{2}L_{\tau}^{\mathrm{train}}}{\partial\phi^{2}}(\phi_{\tau}^{*})+\mathrm{diag}(\lambda)\right)\delta^{\top}+\delta\frac{\partial
L^{\mathrm{eval}}_{\tau}}{\partial\phi}(\phi^{*}_{\tau})$ (165)

and updating the meta-parameters with

$\Delta\omega^{\mathrm{Impl}}_{i}=-\lambda_{i}\delta^{*}_{\tau,i},\quad\mathrm{and}\quad\Delta\lambda^{\mathrm{Impl}}_{i}=(\phi_{\tau,i}^{*}-\omega_{i})\delta^{*}_{\tau,i},$
(166)

where $\delta^{*}_{\tau}$ is the result of the minimization procedure. Using
gradient descent to minimize the quadratic form (165) yields the recurrent
backpropagation algorithm [33, 34, RBP]liao_reviving_2018, also known as
Neumann series approximation in this context [77], we discussed in more
details in Section S3. Gradient descent is a generic algorithm and more
specialized optimization algorithms exist for quadratic forms, such as the
conjugate gradients one used by Rajeswaran et al. [68, iMAML]. Note that the
updates (162), (164) and (166) share some similar structure, but do not
incorporate the teaching signal in the same one. The update of T1-T2 does not
correspond to any gradient, but the other two do. The estimation of the
implicit gradient through Eq. 166 requires computing second derivatives,
whereas the other two updates don’t need to.

The implicit gradient can be approximated in a totally different way, using
the equilibrium propagation theorem [24], as done by the contrastive meta-
learning rule [35]. The idea is to first minimize the so-called augmented loss

$\mathcal{L}_{\tau}(\phi,\theta,\beta):=L^{\mathrm{learn}}_{\tau}(\phi)+\sum_{i}\frac{\lambda_{i}}{2}(\phi_{i}-\omega_{i})^{2}+\beta
L^{\mathrm{eval}}_{\tau}(\phi)$ (167)

for two different values of $\beta$ ($\beta=0$ and some small positive value,
note $\phi_{\tau,0}^{*}$ and $\phi_{\tau,\beta}^{*}$ the solutions) and then
update the meta-parameters with

$\Delta\omega^{\mathrm{CML}}_{i}=\frac{\phi_{\tau,\beta,i}^{*}-\phi_{\tau,0,i}^{*}}{\beta},\quad\mathrm{and}\quad\Delta\lambda^{\mathrm{CML}}_{i}=-\frac{(\phi_{\tau,\beta,i}^{*}-\omega_{i})^{2}-(\phi_{\tau,0,i}^{*}-\omega_{i})^{2}}{2\beta}.$
(168)

The equilibrium propagation theorem guarantees that those updates will
converge to the true gradient when $\beta\rightarrow 0$. Both the contrastive
meta-learning rule and the rule prescribed by our principle are first-order,
as that they don’t require computing second derivative. However, the former
approximates the implicit gradient when the latter approximates the gradient
of the least-control objective $\mathcal{H}$. Interestingly, as we mentioned
in Section 2.4, the augmented solution $\phi^{*}_{\beta}$ and the controlled
equilibrium $\phi_{*}^{\tau}$ are both stationary points of the augmented loss
$\mathcal{L}_{\tau}$ (taking $\beta=\alpha^{-1}$ for $\phi_{*}^{\tau}$).
However, the contrastive meta-learning update (168) is justified in the weak
nudging limit ($\beta\rightarrow 0$), whereas the update (162) is justified in
the perfect control limit ($\alpha=\beta^{-1}\rightarrow 0$).

We refer to zucchet_beyond_2022 for a more detailed review of those methods.

Reptile. The Reptile algorithm [76] meta-learns the initialization of a neural
network such that it can adapt to a new task $\tau$ in very few gradient
descent steps. Let us denote $\omega$ the parameters of the network at
initialization. The Reptile algorithm first minimizes the learning loss
$L^{\mathrm{learn}}_{\tau}(\phi)$ and then update $\omega$ through

$\Delta\omega^{\mathrm{Reptile}}=\phi_{\tau}^{*}-\omega,$ (169)

with $\phi_{\tau}^{*}$ the result of the learning loss minimization. Contrary
to the methods presented above, this algorithm is only heuristically motivated
and does not rely on any theoretical foundations. Still, it performs
surprisingly well on many few-shot learning tasks.

####  S7.3 Experimental details

We compare the LCP to existing meta-learning and few-shot learning algorithms
on the Omniglot few-shot image classification benchmark [78]. We follow the
standard training and evaluation used in prior works [74, 68]. We focus on the
20-way 1-shot setting, and only meta-learn $\theta=\\{\omega\\}$ while
choosing a global $\lambda$ as a hyperparameter. The iMAML result in Table 3
is taken from Rajeswaran et al. [68] while those of Reptile and FOMAML are
taken from Nichol et al. [76].

Data preprocessing and augmentation. We resize the Omniglot images to
$28\times 28$ and augment the images by random rotation by multiples of $90$
degrees during meta-training.

Architecture. For our classifier, we use the same architecture as in previous
works [74, 68], which consists in 4 convolutional modules with a $3\times 3$
convolutions and 64 filters, each followed by a batch normalization layer, a
ReLU nonlinearity, and a $2\times 2$ max-pooling layer. The output of the
module is then flattened and fed into a final classification layer followed by
a softmax activation.

Optimization. We train all models on 100000 tasks, split into meta batches of
size 16. Each task consists in a 20-way classification problem with a single
example per class provided for the adaptation phase, and 15 examples per class
for evaluating the adaptation. For T1-T2, the adaptation consists in 100
gradient steps with learning rate $\eta=0.03$. For LCP, we use 100 steps of
the discretized version of the dynamics (157):

$\begin{split}\phi_{\mathrm{next}}&=\phi+\eta\left(-\nabla_{\phi}L^{\mathrm{learn}}_{\tau}(\phi)-\lambda(\phi-\omega)+u\right)\\\
u_{\mathrm{next}}&=u+\eta\left(-\alpha
u-\nabla_{\phi}L_{\tau}^{\mathrm{eval}}(\phi)\right)\\!,\end{split}$ (170)

where $\eta=0.03$. The Adam optimizer takes the $\Delta\omega$ as input and
updates the meta-parameters, with learning rate 0.001. Finally, the models are
evaluated on 100 test tasks. See Table S9 for an overview of the
hyperparameters.

Table S9: Hyperparameter used for the 20-way 1-shot Omniglot experiment. Hyperparameter | T1-T2 | LCP
---|---|---
num_tasks | 100000 | 100000
meta_batch_size | 16 | 16
optimizer_outer | ADAM | ADAM
lr_outer | 0.001 | 0.001
steps_inner | 100 | 100
$\eta$ | 0.03 | 0.03
$\alpha$ | - | 0.1
$\lambda$ | 0.001 | 0.001

unsrtnat references