SMSupplementary References

# Constrained Predictive Coding as a Biologically Plausible Model of the
Cortical Hierarchy

Siavash Golkar$\thanks{Equal contribution}\hskip 3.0pt^{\,1}$
&Tiberiu Tesileanu$\hskip 3.0pt^{\,1}$11footnotemark: 1
&Yanis Bahroun 1
Anirvan M. Sengupta 2,3,4 Dmitri B. Chklovskii 1,5
${}^{1\,}$Center for Computational Neuroscience, Flatiron Institute
${}^{2\,}$Center for Computational Mathematics, Flatiron Institute
${}^{3\,}$Center for Computational Quantum Physics, Flatiron Institute
${}^{4\,}$Department of Physics and Astronomy, Rutgers University
${}^{5\,}$Neuroscience Institute, NYU Medical Center
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>Equal contribution

###### Abstract

Predictive coding (PC) has emerged as an influential normative model of neural
computation with numerous extensions and applications. As such, much effort
has been put into mapping PC faithfully onto cortex, but there are issues that
remain unresolved or controversial. In particular, current implementations
often involve separate value and error neurons and require symmetric forward
and backward weights across different brain regions. These features have not
been experimentally confirmed. In this work, we show that the PC framework in
the linear regime can be modified to map faithfully onto the cortical
hierarchy in a manner compatible with empirical observations. Employing a
disentangling-inspired constraint on hidden-layer neural activities, we derive
an upper bound for the PC objective. Optimization of this upper bound leads to
an algorithm that shows the same performance as the original objective and
maps onto a biologically plausible network. The units of this network can be
interpreted as multi-compartmental neurons with non-Hebbian learning rules,
with a remarkable resemblance to recent experimental findings. There exist
prior models which also capture these features, but they are phenomenological,
while our work is a normative derivation. Notably, the network we derive does
not involve one-to-one connectivity or signal multiplexing, which the
phenomenological models required, indicating that these features are not
necessary for learning in cortex. The normative nature of our algorithm in the
simplified linear case allows us to prove interesting properties of the
framework and analytically understand the computational role of our network’s
components. The parameters of our network have natural interpretations as
physiological quantities in a multi-compartmental model of pyramidal neurons,
providing a concrete link between PC and experimental measurements carried out
in cortex.

##  1 Introduction

Over the past decades, predictive coding (PC), a normative framework for
learning representations that maximize predictive power, has played an
important role in computational neuroscience [1, 2]. Initially proposed as an
unsupervised learning paradigm in the retina [3], it has since been expanded
to the supervised regime [4] with arbitrary graph topologies [5, 6]. The PC
framework has been analyzed in many contexts [7, 8] and has found many
applications, from clinical neuroscience [9, 10] to memory storage and
retrieval [11]. We refer the reader to [12, 13] for recent reviews.

Figure 1: Schematic architecture of the predictive coding network (PCN) and
our covariance-constrained network (BioCCPC). (a) PCN from [4], figure adapted
from [13]. The intra-layer connectivity is to be one-to-one, while the inter-
layer connectivity is symmetric. (b) Our BioCCPC network. There is no
requirement for symmetric weights across layers, and the connectivity within
layers can be diffuse.

Predictive coding is viewed as a possible theory of cortical computation, and
many parallels have been drawn with the known neurophysiology of cortex [12].
While the initial works proposed a biologically plausible network [1, 2], the
connection with cortex was more closely examined in [14], where the PC module
was mapped onto a cortical-column microcircuit. However, there are aspects of
this mapping that have proved controversial. Among these are the requirement
of multiple redundant cortical operations, the symmetric connectivity pattern,
the one-to-one connectivity of value and error neurons, and also the
requirement that feedback connections be inhibitory [12, 15, 16], as sketched
in Fig. 1a . The presence of separate error and value neurons has itself been
called into question [15].

The PC-based neural circuits also do not account for more recent experimental
findings highlighting the details of computation in the cortex [17, 18, 19,
20, 21, 22, 23, 24, 25, 26]. For example, the learning dynamics of multi-
compartmental pyramidal neurons has been closely investigated. In particular,
it was observed that the plasticity of the synapses of the basal compartment
is driven by the activity in the apical tuft of the neuron by so-called
calcium plateau potentials [20, 19, 24], leading to non-Hebbian learning rules
[25]. These experiments have motivated the development of several models of
microcircuits with multi-compartmental neurons [27, 28, 29, 30, 31, 32, 33].
In a number of cases, it has been shown that these models can replicate
learning similar to the backpropagation algorithm under specific assumptions
[34, 35, 30]. However, because of their rather phenomenological nature,
detailed analysis of these models is in many cases challenging and one must
resort to numerical simulations, rendering the task of understanding the role
of various neurophysiological quantities difficult. To this date, a normative
framework (PC or otherwise) that can explain these experimental findings is
still lacking.

In this work, we show how to modify the PC framework to make it compatible
with the aforementioned experimental observations while maintaining much of
the analytical control afforded to use by a normative approach. Inspired by
prior work which explored the effects of finding decorrelated representations
[36, 37, 38, 39, 40, 41], we add a decorrelating inequality constraint to the
covariance of the latent representations of PC. Using this constraint, we
derive an upper bound for the PC objective. By working in the linear regime,
we can prove several useful properties of our algorithm. We show that the
learning algorithm derived from this upper bound does not suffer from the
issues of prior implementations and naturally maps onto the hierarchical
structure of the cortical pyramidal neurons.

#### Contributions

  * •

We introduce BioCCPC, a novel biologically plausible algorithm derived as a
modification of predictive coding (PC), and show that it closely matches the
known physiology of the cortical hierarchy.

  * •

We find an analytic upper bound to the PC objective by imposing a
decorrelation-inspired inequality constraint on the latent space.

  * •

We interpret the different parameters of the algorithm in terms of the
conductances and leaks of the separate compartments of the pyramidal neuron.
We find that the neural compartmental conductances encode the variances of the
PC framework, and the somatic leak maps onto a thresholding mechanism of the
associated eigenvalue problem.

##  2 Related work and review of predictive coding

The backpropagation algorithm [42, 34] is the predominant tool in the training
of deep neural networks but is generally not considered biologically plausible
[43]. Over the years, many authors have explored biologically plausible
approximations or alternatives to backpropagation [44, 45, 46, 47, 48, 49, 30,
50, 51, 52, 31, 53, 54, 55, 56, 57, 58] (for a more complete review see [59]).
These approaches generally fall into two categories. First are normative
approaches, such as Predictive Coding [1, 2, 3], Target Propagation and
variations [44, 45, 47], Equilibrium Propagation [49, 53] and others [54, 51,
55, 56], where one starts from a mathematically motivated optimization
objective. These methods, by virtue of their normative derivation, have a firm
grounding in theory; however, they do not fully conform to the experimental
observations in the brain (see below and Sec. 4). The second approach is
driven by biology, with network structures and learning rules inspired by
experimental findings [28, 48, 60, 30, 31]. While these works mostly conform
to the experimentally observed findings, they are more challenging to analyze
because of their conjectured phenomenological nature.

The goal of the present paper is _not_ to propose yet another biologically
plausible alternative to backpropagation. It is rather to demonstrate that the
normative framework of predictive coding, when combined with a constraint, can
indeed closely match experimental observations. For this reason, in this work,
we focus on comparing our method with previous implementations of predictive
coding and do not concern ourselves with other biologically plausible
alternatives to backpropagation. The relationship between the PC framework and
backpropagation was explored in [4, 61, 62, 63, 64]. The advantages of PC over
backpropagation were highlighted in [13].

We now turn to a brief review of the predictive coding framework from [4]
before delving into the details of our method.

#### Notation.

Bold upper case ${\bf M}$ and lower case variables ${\bf v}$ denote matrices
and vectors, respectively. By upper case letter ${\bf X},{\bf Y},{\bf Z}$, we
denote data matrices in $\mathbb{R}^{d\times T}$, where $d$ and $T$ are the
dimensions of the relevant variable and the number of samples. Lower case
${\bf x},{\bf y},{\bf z}$, denote the relevant quantities of a single sample,
and ${\bf x}_{t},{\bf y}_{t},{\bf z}_{t}$, denote the $t^{\text{th}}$ sample.
$\|{\bf M}\|^{2}_{F}$ denotes the Frobenius norm of M.

###  2.1 Review of predictive coding

#### Probabilistic model.

In this section, we review the supervised predictive coding algorithm from
[4]. The derivation starts from a probabilistic model for supervised learning,
which parallels the architecture of an artificial neural network (ANN) with
$n+1$ layers. In this model, the neurons of each layer are random variables
${\bf z}^{(l)}$ (denoting the vector of activations in the $l^{\text{th}}$
layer) with layers $0$ and $n$, respectively, denoting the input and output
layers of the network. We assume that the joint probability of the latent
variables factorizes in a Markovian structure

$p({\bf z}^{(0)},{\bf z}^{(1)},\cdots,{\bf z}^{(n)})=p({\bf z}^{(n)}|{\bf
z}^{(n-1)})\times p({\bf z}^{(n-1)}|{\bf z}^{(n-2)})\times\cdots\times p({\bf
z}^{(0)})$

with the relationship between the random variables of adjacent layers given
by:

$p\bigl{(}{\bf z}^{(l)}|{\bf z}^{(l-1)}\bigr{)}=\mathcal{N}\bigl{(}{\bf
z}^{(l)};\bm{\mu}^{(l)},\bm{\Sigma}^{(l)}\bigr{)}\,,\qquad\text{with}\quad\bm{\mu}^{(l)}={\bf
W}^{(l-1)}f({\bf z}^{(l-1)}),\quad\bm{\Sigma}^{(l)}=\sigma^{(l)2}{\bf I}\,.$
(1)

The mean of the probability density on layer $l$ mirrors the activity of the
analogous ANN given by $\bm{\mu}^{(l)}={\bf W}^{(l-1)}f({\bf z}^{(l-1)}),$
where ${\bf W}^{(l-1)}$ are the weights connecting layers $l-1$ and $l$. The
objective function is then given by the negative log-likelihood of the joint
distribution function:

$\displaystyle L=-\sum_{t}\log p({\bf z}_{t}^{(0)},\dotsc,{\bf z}_{t}^{(n)})$
$\displaystyle=\frac{1}{2}\sum_{l}\frac{\bigl{\|}{\bf Z}^{(l)}-{\bf
W}^{(l-1)}f({\bf
Z}^{(l-1)})\bigr{\|}_{F}^{2}}{\sigma^{(l)2}}+\text{const}~{},$ (2)

where we have switched to the data matrix notation for brevity and assumed the
variances $\sigma^{(l)2}$ are fixed hyperparameters. In the following, we
refer to $L$, Eq. (2), without the constant term.

#### Learning.

Learning in the Whittington and Bogacz framework takes place in two steps [4].
First, the values of the random variables are determined by finding the most
probable configuration of the joint distribution function when both the input
and output layers are conditioned on the given input and output (${\bf
z}^{(0)}={\bf x}$ and ${\bf z}^{(n)}={\bf y}$):

${\bf z}^{*(1)},\dotsc,{\bf z}^{*(n-1)}=\underset{{\bf z}^{(1)},\dotsc,{\bf
z}^{(n-1)}}{\operatorname{arg}\,\operatorname{min}}\;L({\bf z}^{(0)}={\bf
x},{\bf z}^{(n)}={\bf y}).$ (3)

The solution to this minimization problem can be found via gradient descent,
which we evaluate component-wise for clarity as:

$\dot{z}^{(l)}_{j}=-\eta\partial_{z_{j}}L=\eta(-\varepsilon^{(l)}_{j}+\sum_{i}\varepsilon^{(l+1)}_{i}W_{ij}^{(l)}f^{\prime}(z_{j}^{(l)})),\quad\text{with}\quad\varepsilon^{(l)}_{i}=\frac{z_{i}^{(l)}-\mu_{i}^{(l)}}{\sigma^{(l)2}}\,,$
(4)

where $\eta$ is the gradient descent step size. The second step is then to
minimize the objective with respect to the weights while keeping the
previously obtained neuron values fixed. This corresponds to optimizing the
value of the loss at the MAP estimate and can also be carried out by gradient
descent, leading to:

$\delta W_{ij}^{(l)}\propto-\frac{\partial L}{\partial
W_{ij}^{(l)}}=\varepsilon^{*(l+1)}_{i}f\bigl{(}z^{*(l)}_{j}\bigr{)},$ (5)

where the star superscript denotes equilibrium values of the latent variable
dynamics obtained in the first step (Eq. (3)). This algorithm can be
implemented by a biologically plausible network as described in [4]; see
Figure 1a. However, as discussed in Section 1, its mapping onto the cortex has
proved controversial. For further details regarding these steps see [4].

##  3 A constrained predictive coding framework

In this section, we introduce and discuss our novel covariance-constrained
predictive coding (CCPC) model that builds upon the supervised PC paradigm of
[4]. Our model also straightforwardly extends to the unsupervised PC framework
of [1, 2]. For simplicity, we work in the linear regime ($f(x)=x$), which
allows us to prove different properties of our framework.

###  3.1 Derivation of an upper bound objective

#### Reduction to a sum of objectives.

We start by reducing the optimization problem, Eq. (2), into a set of
overlapping sub-problems, which will allow us to break the symmetry between
feedforward and feedback weights. To do so, we first introduce a copy of the
terms containing the weights ${\bf W}^{(1)}$ to ${\bf W}^{(n-2)}$, denoted by
${\bf W}_{a}^{(l)}$ and ${\bf W}_{b}^{(l)}$ respectively, as

$\begin{split}\min_{{\bf Z},{\bf W}}L&=\min_{{\bf Z},{\bf
W}}\frac{1}{2}\sum_{l=1}^{n}\frac{\bigl{\|}{\bf Z}^{(l)}-{\bf W}^{(l-1)}{\bf
Z}^{(l-1)}\bigr{\|}_{F}^{2}}{\sigma^{(l)2}}=\min_{{\bf Z},{\bf W}_{a},{\bf
W}_{b}}\hat{L}\\\ \hat{L}&=\min_{{\bf Z},{\bf W}_{a},{\bf
W}_{b}}\frac{1}{2}\frac{\bigl{\|}{\bf Z}^{(1)}-{\bf W}_{b}^{(0)}{\bf
Z}^{(0)}\bigr{\|}_{F}^{2}}{\sigma^{(1)2}}+\frac{1}{2}\frac{\bigl{\|}{\bf
Z}^{(n)}-{\bf W}_{a}^{(n-1)}{\bf
Z}^{(n-1)}\bigr{\|}_{F}^{2}}{\sigma^{(n)2}}\\\
&\qquad+\frac{1}{4}\sum_{l=2}^{n-1}\left[\frac{\bigl{\|}{\bf Z}^{(l)}-{\bf
W}_{b}^{(l-1)}{\bf
Z}^{(l-1)})\bigr{\|}_{F}^{2}}{\sigma^{(l)2}}+\frac{\bigl{\|}{\bf Z}^{(l)}-{\bf
W}_{a}^{(l-1)}{\bf
Z}^{(l-1)})\bigr{\|}_{F}^{2}}{\sigma^{(l)2}}\right].\end{split}$ (6)

For consistency, we rename ${\bf W}^{(0)}$, ${\bf W}^{(n-1)}$ to ${\bf
W}_{b}^{(0)}$, ${\bf W}_{a}^{(n-1)}$, respectively. Introducing these copies
does not change the optimization111This can be directly verified by finding
the optima for ${\bf W}$’s before and after the change. We have ${\bf
W}^{(l)}={\bf W}_{a}^{(l)}={\bf W}_{b}^{(l)}={\bf Z}^{(l+1)}{\bf
Z}^{(l)\top}({\bf Z}^{(l+1)}{\bf Z}^{(l)\top})^{-1}$. Plugging these back into
Eq. (6) we see that the equality holds. However, in the next step, since we
treat ${\bf W}_{a}$’s and ${\bf W}_{b}$’s differently, ${\bf W}_{a}={\bf
W}_{B}$ will no longer hold. but will help us avoid weight sharing in the
steps below. We now pair the terms two by two as

$\displaystyle\min_{{\bf Z},{\bf W}_{a},{\bf
W}_{b}}\hat{L}=\frac{1}{2}\sum_{l=1}^{n-1}\left[g_{b}^{(l)}\bigl{\|}{\bf
Z}^{(l)}-{\bf W}_{b}^{(l-1)}{\bf
Z}^{(l-1)}\bigr{\|}_{F}^{2}+g_{a}^{(l)}\bigl{\|}{\bf Z}^{(l+1)}-{\bf
W}_{a}^{(l)}{\bf Z}^{(l)}\bigr{\|}_{F}^{2}\right]~{},$ (7)

with $g_{a}^{(n-1)}=1/\sigma^{(n)2}$, $g_{b}^{(1)}=1/\sigma^{(1)2}$, and
$g_{a}^{(l-1)}=g_{b}^{(l)}=1/(2\sigma^{(l)2})$ for $l=2,\dotsc,n-1$.

Weight sharing occurs here from terms like ${\bf z}^{(l+1)\top}{\bf W}{\bf
z}^{(l)}$, obtained from expanding the squared norms in Eq. (7). Indeed, the
gradient descent dynamics with respect to ${\bf z}^{(l+1)}$ (resp. ${\bf
z}^{(l)}$) leads to terms of the form ${\bf W}{\bf z}^{(l)}$ in $\dot{\bf
z}^{(l+1)}$ (resp. ${\bf W}^{\top}{\bf z}^{(l+1)}$ in $\dot{\bf z}^{(l)}$),
which use the same weights ${\bf W}$. Thanks to the doubling of the weights,
we can avoid this problem by optimizing each term in the sum in Eq. (7)
separately. In other words,

$\displaystyle\min_{{\bf Z},{\bf W}_{a},{\bf
W}_{b}}\hat{L}\leq\sum_{l=1}^{n-1}\min_{{\bf Z}^{(l)},{\bf W}_{a}^{(l)},{\bf
W}_{b}^{(l-1)}}\hat{L}^{(l)}~{},$ (8)
$\displaystyle\text{where}~{}~{}\hat{L}^{(l)}\equiv\frac{1}{2}\sum_{l=0}^{n}\left[g_{b}^{(l)}\bigl{\|}{\bf
Z}^{(l)}-{\bf W}_{b}^{(l-1)}{\bf
Z}^{(l-1)}\bigr{\|}_{F}^{2}+g_{a}^{(l)}\bigl{\|}{\bf Z}^{(l+1)}-{\bf
W}_{a}^{(l)}{\bf Z}^{(l)}\bigr{\|}_{F}^{2}\right].$

This inequality holds simply because we are no longer finding the minimum of
the full objective $\hat{L}$. We are instead finding the minimum of each
component separately and then evaluating $\sum_{l}\hat{L}^{(l)}=\hat{L}$. This
splits the $(n+1)$-layer optimization problem into a set of $3$-layer
optimizations, in each of which only the middle layer is being optimized.
Note, however, that these are overlapping, so the different optimization
problems need to be solved self-consistently. We make this precise in the
supplementary materials section and show that it provides an upper bound for
our objective $L$ (SM Sec. A). Separating the objective function in this
manner eliminates the weight sharing problem for ${\bf W}_{b}$, but the
problem remains for ${\bf W}_{a}$. This can be seen by noting that the
gradient descent equations only depend on ${\bf W}_{b}$ but depend on both
${\bf W}_{a}$ and ${\bf W}_{a}^{\top}$. We address this problem next.

#### Whitening constraint.

The idea of decorrelating internal representation has been widely used for
unsupervised tasks, often motivated by neuroscience [36, 37, 38]. In the case
of deep learning, the main motivations were improved convergence speed and
generalization [39, 40, 41]. Decorrelation has also been used to circumvent
the weight transport problem [55]. Inspired by these observations, we
introduce the inequality constraint $\frac{1}{T}{\bf Z}^{(l)}{\bf
Z}^{(l)\top}\preceq{\bf I}$, imposing an upper bound on the eigenvalues of the
covariance matrix. In Sec. 5 and SM Sec. B we show that under some assumptions
this inequality constraint leads to whitening the latent variables. Adding
this inequality constraint to the objective we have:

$\displaystyle\min_{{\bf Z}^{(l)},{\bf W}_{a}^{(l)},{\bf
W}_{b}^{(l)}}\hat{L}^{(l)}\leq$ $\displaystyle\min_{\begin{subarray}{c}{\bf
Z}^{(l)},{\bf W}_{a}^{(l)},{\bf W}_{b}^{(l)}\\\ {\bf Z}^{(l)}{\bf
Z}^{(l)\top}\preceq T\times{\bf I}\end{subarray}}\hat{L}^{(l)}$
$\displaystyle\leq$ $\displaystyle\min_{\begin{subarray}{c}{\bf Z}^{(l)},{\bf
W}_{a}^{(l)},{\bf W}_{b}^{(l)}\\\ {\bf Z}^{(l)}{\bf Z}^{(l)\top\preceq
T\times{\bf I}}\end{subarray}}\frac{1}{2}\biggl{[}g_{b}^{(l)}\Bigl{\|}{\bf
Z}^{(l)}-{\bf W}_{b}^{(l-1)}{\bf
Z}^{(l-1)}\Bigr{\|}_{F}^{2}+c^{(l)}\bigl{\|}{\bf Z}^{(l)}\bigr{\|}_{F}^{2}$
$\displaystyle\qquad+g_{a}^{(l)}\operatorname{Tr}\Bigl{(}-2{\bf
Z}^{(l+1)\top}{\bf W}_{a}^{(l)}{\bf Z}^{(l)}\Bigr{)}+g_{a}^{(l)}T\,{\bf
W}_{a}^{(l)\top}{\bf W}_{a}^{(l)}$

The first inequality results from the fact that the parameters of the new
optimization problem are a subset of the parameters of the original
optimization problem. In the second inequality we have expanded the second
term of $\hat{L}^{(l)}$ and used the constraint $\frac{1}{T}{\bf Z}^{(l)}{\bf
Z}^{(l)\top}\preceq{\bf I}$ to replace ${\bf Z}^{(l)}{\bf Z}^{(l)\top}$ with
$T$, making the term larger, resulting again in an upper bound. This
replacement eliminates the weight sharing problem in ${\bf W}_{a}$, as the
gradient descent equations now only depend on ${\bf W}_{a}^{\top}$ and not on
${\bf W}_{a}$ (cf. Eqs. (10) and (11)).

The inequality can be implemented by using a positive-definite Lagrange
multiplier ${\bf Q}^{(l)\top}{\bf Q}^{(l)}$:

$\min_{{\bf Z}^{(l)},{\bf W}_{a}^{(l)},{\bf W}_{b}^{(l-1)}}\max_{{\bf
Q}^{(l)}}\frac{1}{2}\biggl{[}g_{b}^{(l)}\Bigl{\|}{\bf Z}^{(l)}-{\bf
W}_{b}^{(l-1)}{\bf Z}^{(l-1)}\Bigr{\|}_{F}^{2}+c^{(l)}\bigl{\|}{\bf
Z}^{(l)}\bigr{\|}_{F}^{2}\\\ +g_{a}^{(l)}\operatorname{Tr}\Bigl{(}-2{\bf
Z}^{(l+1)\top}{\bf W}_{a}^{(l)}{\bf Z}^{(l)}\Bigr{)}+g_{a}^{(l)}T\,{\bf
W}_{a}^{(l)\top}{\bf W}_{a}^{(l)}\\\ +\operatorname{Tr}{\bf Q}^{(l)\top}{{\bf
Q}^{(l)}}\bigl{(}{{\bf Z}^{(l)}}{\bf Z}^{(l)\top}-T{\bf
I}\bigr{)}+\text{const}\biggr{]}\,,$ (9)

where we have added an additional quadratic term in ${\bf Z}$ as a
regularizer. This can also be interpreted as an extra term in the prior of the
latent space. Equation (9) is the final form of the objective, which we then
optimize using gradient descent/ascent.

###  3.2 Neural dynamics and learning rules.

Similar to the PCN of [4], the dynamics of our network during learning
proceeds in two steps. First the neural dynamics is derived by taking gradient
steps of $\hat{L}^{(l)}$ from Eq. (9) with respect to ${\bf z}^{(l)}$:

$\displaystyle\dot{\bf z}^{(l)}=g_{b}^{(l)}{\bf W}_{b}^{(l-1)}{\bf
z}^{(l-1)}+g_{a}^{(l)}{\bf W}_{a}^{(l)\top}{\bf
z}^{(l+1)}-\bigl{(}g_{b}^{(l)}+c^{(l)}\bigr{)}{\bf z}^{(l)}-g_{a}^{(l)}{\bf
Q}^{(l)\top}{\bf n}^{(l)},$ (10)

where we have defined the variables ${\bf n}^{(l)}=(1/g_{a}^{(l)}){\bf
Q}^{(l)}{\bf z}^{(l)}$ for each layer of the network.

The weight updates are derived via stochastic gradient descent of the loss
given in Eq. (9) after the neural dynamics have reached equilibrium. These are
given by

$\displaystyle\delta{\bf W}_{b}^{(l)}$
$\displaystyle\propto\Big{[}g_{a}^{(l+1)}\bigl{(}{\bf W}_{a}^{(l+1)\top}{\bf
z}^{(l+2)}-{\bf Q}^{(l+1)\top}{\bf n}^{(l+1)}\bigr{)}-c^{(l+1)}{{\bf
z}^{(l+1)}}\Big{]}{{\bf z}^{(l)\top}}\,,$ (11a) $\displaystyle\delta{\bf
W}_{a}^{(l)\top}$ $\displaystyle\propto{\bf z}^{(l)}{{\bf z}^{(l+1)\top}}-{\bf
W}_{a}^{(l)\top}\,,$ (11b) $\displaystyle\delta{\bf Q}^{(l)}$
$\displaystyle\propto{\bf n}^{(l)}{\bf z}^{(l)\top}-{\bf Q}^{(l)}\,.$ (11c)

We used the neural dynamics equilibrium equation for ${\bf z}^{(l)}$ to
simplify the weight update for ${\bf W}_{b}^{(l)}$. This yields our online
algorithm (Alg. 1 in SM Sec. E), with the architecture shown in Fig. 1b. The
algorithm can be implemented in a biologically plausible neural network as in
Fig. 2; see Sec. 4.

##  4 Biological implementation and comparison with experimental observations

In this section, we introduce a biologically plausible neural circuit that
implements the CCPC algorithm, denoted by BioCCPC. We also demonstrate that
the details of this circuit resemble the neurophysiological properties of
pyramidal cells in the neocortex and the hippocampus.

#### Neural architecture.

The algorithm for BioCCPC (Alg. 1) summarized by the neural dynamics from Eq.
(10) and weight update rules from Eqs. (11) can be implemented by a neural
circuit with schematic shown in Fig. 2.

The activity of the $n-1$ hidden layers, $\\{{\bf z}^{(l)}\\}_{l=1}^{n-1}$, is
encoded as the outputs of $n-1$ sets of neurons, representing pyramidal
neurons of different cortical regions. The matrices ${\bf W}_{a}^{(l)}$ and
${\bf W}_{b}^{(l)}$ are encoded as the weights of synapses between the
pyramidal neurons of adjacent layers. Explicitly, the matrix ${\bf
W}_{a}^{(l)}$ (resp. ${\bf W}_{b}^{(l-1)}$) is the efficacy of the synapse
connecting ${\bf z}^{(l+1)}$ (resp. ${\bf z}^{(l-1)}$) to the pyramidal
neurons ${\bf z}^{(l)}$. Because of the disjoint nature of these two inputs,
we model these as synapsing respectively onto the distal (apical tuft) and
proximal (mostly basal) dendrites of the pyramidal neurons, respectively; see
Fig. 2. This is reminiscent of cortical pyramidal neurons, which also have two
integration sites, the proximal compartment comprised of the basal and
proximal apical dendrites providing inputs to the soma, and the distal
compartment comprised of the apical dendritic tuft [65, 22]. These two
compartments receive excitatory inputs from two separate sources [66, 18].

Similarly, the auxiliary variables ${\bf n}^{(l)}$ are represented by the
activity of interneurons in each cortical region.222There are multiple types
of interneurons targeting pyramidal cells [67, 68]. The interneurons of
BioCCPC most closely resemble the somatostatin-expressing interneurons, which
preferentially inhibit the apical dendrites. The ${\bf Q}^{(l)}$ synaptic
weights are encoded in the weights of synapses connecting ${\bf n}^{(l)}$ to
${\bf z}^{(l)}$, while ${\bf Q}^{(l)\top}$ models the weights of synapses from
${\bf z}^{(l)}$ to ${\bf n}^{(l)}$. In a biological setting, the implied
equality of weights of synapses from ${\bf z}^{(l)}$ to ${\bf n}^{(l)}$ and
the transpose of those from ${\bf n}^{(l)}$ to ${\bf z}^{(l)}$ can be
guaranteed by applying the same Hebbian learning rule [55]. However, unlike
previous work, we do not require interneurons of one layer to be connected to
the pyramidal neurons of another layer.

Figure 2: Schematic of the biological implementation of the BioCCPC algorithm.
Left: overall connectivity pattern, showing all-to-all connections between
pyramidal neurons in consecutive layers, and between pyramidal and
interneurons in the same layer. Right: zoom in on one neuron and its immediate
neighbors, showing the membrane potentials and synaptic weights that are
relevant for our algorithm.

#### Neural dynamics and canonical components.

The neural dynamics of the pyramidal neurons of our circuit given by Eq. (10)
can be recast as the dynamics of a three compartment neuron:

$\displaystyle\tau\dot{\bf z}^{(l)}=-g^{(l)}_{\text{lk}}{\bf
z}^{(l)}+g_{a}^{(l)}({\bf v}_{a}^{(l)}-{\bf z}^{(l)})+g_{b}^{(l)}({\bf
v}_{b}^{(l)}-{\bf z}^{(l)}),$ (12)

where we have defined the apical and basal compartmental membrane potentials

${\bf v}_{a}^{(l)}={\bf W}_{a}^{(l)\top}{\bf z}^{(l+1)}-{\bf Q}^{(l)\top}{\bf
n}^{(l)}\;,\qquad{\bf v}_{b}^{(l)}={\bf W}_{b}^{(l-1)}{\bf z}^{(l-1)}.$ (13)

Here we defined the leak conductance
$g_{\text{lk}}^{(l)}=c^{(l)}-g_{a}^{(l)}$. Interestingly, this is the same
neural dynamics of a three-compartment neuron posited previously in [28, 30],
which we derived in a normative way. Here, $g_{a}$ (resp. $g_{b}$) are
conductances between the apical (resp. basal) compartments and the soma, and
$g_{\text{lk}}$ is the somatic leak. Our derivation clarifies the relationship
between the physiological quantities and the objective function. Indeed, the
basal and apical conductances play the role of the inverse variances in the
predictive coding objective. These inverse variances, generally referred to as
‘precisions’ in the PC literature, have previously been argued to be encoded
in lateral inhibitions [2] or top-down attention [69, 70]. In our model, they
are simply encoded in the compartmental conductances. In Sec. 5 we show that
the leak term also has a functional role and determines a threshold for
dynamically learning more diverse variables [71, 72].

#### Synaptic weight updates.

Matching observations in the cortex, in our circuit [18, 73, 74, 65, 22], the
basal and apical weights (${\bf W}_{b}$ and ${\bf W}_{a}$) of our algorithm
are updated differently. The basal synaptic weights of the pyramidal neurons,
given by the elements of ${\bf W}_{b}^{(l)}$, are updated by the product of
two factors represented in the corresponding post- and pre-synaptic neurons
(Eq. 11a):

$\delta{\bf W}_{b}^{(l-1)}\propto\Big{[}g_{a}^{(l)}{\bf
v}^{(l)}_{a}-c^{(l)}{\bf z}^{(l)}\Big{]}{\bf z}^{(l-1)\top}\,.$ (14)

The first term in parentheses, $g_{a}^{(l)}{\bf v}^{(l)}_{a}$, is proportional
to the total apical current, given by the difference between the excitatory
synaptic current in the apical tuft, ${\bf W}_{a}^{(l)\top}{\bf z}^{(l+1)}$,
and the inhibitory current induced by interneurons synapsing onto the distal
compartment, ${\bf Q}^{(l)\top}{\bf n}^{(l)}$. Biologically, this factor can
be approximated by the calcium plateau potential traveling down the apical
shaft. The calcium plateau potential has been experimentally seen to drive
plasticity of the basal synapses matching the derived update rule of our
circuit [18, 20, 23, 24, 25]. The second factor is simply the neural output.
Because the total update is not purely dependent on the action potentials of
the pre- and post-synaptic neurons, such plasticity is called non-Hebbian
[25]. The synaptic learning rule for the apical weights ${\bf W}_{a}^{(l)}$
(Eq. 11b) and the synapses connecting the interneurons and the pyramidal
neurons ${\bf Q}^{(l)}$ (Eq. 11c) are simply Hebbian, also matching
experimental observations in the cortex [19, 75].

#### Features and limitations.

Here we summarize the main differences between the biological implementations
of CCPC and PC. In PC, the forward and backward weights (in our notation ${\bf
W}_{b}$ and ${\bf W}_{a}^{\top}$) are symmetric, whereas, in CCPC, they are
not constrained to be symmetric and are generically not so (see Sec. 6).
Another difference is that in PC, the connection between value neurons and
error neurons is one-to-one and fixed; this constraint does not exist between
the pyramidal and interneurons of our model, which are no longer one-to-one,
and have plastic connectivity.

While we have shown the improvement in biological realism of our model over
the traditional PC network, this was done in the linear case. However, recent
work on deep linear networks [76] has provided many insights into the learning
dynamics of deep networks. Some of the properties discovered in the deep
linear network, like “balancedness" of weights [77, 78], generalize for
certain nonlinear networks [78]. We believe several of our observations will
generalize in a similar fashion.

There are other aspects that can also be improved. For example, the
connectivity between pyramidal neurons and interneurons of CCPC is required to
be symmetric. This symmetry can be achieved via Hebbian learning [55];
however, it would be interesting to explore if such symmetry is indeed
required. In CCPC, symmetric connectivity is only required between neurons of
the _same_ cortical region, which is a much less stringent biological
requirement than that of symmetric weights between _different_ cortical
regions, as in PC. Another shortcoming of our model is that while the teaching
signal for the basal synapses in CCPC is signed and graded, in the cortex
these signals are generally believed to be stereotypical [18].

##  5 Theoretical arguments and interpretation

In this section, we summarize some of the salient theoretical features of our
framework and describe their biological significance. Further details and
proofs are given in SM Sec. B.

###  5.1 Implicit error computation

In backpropagation, as well as in PC, the learning algorithm computes a loss
(or a local error) which is then used to compute updates to the weights. In
CCPC, no error or loss is explicitly computed. So how does CCPC learn? In the
following proposition, we show that in the $g_{a}\to 0$ limit, the algorithm
implicitly computes an error in its weight updates. This quasi feed-forward or
weak nudging limit has been used to explore the learning dynamics of
biologically plausible networks [30] and is related to the ‘fixed prediction’
assumption [4, 61, 62].

######  Proposition 1.

Assume that the learning rules (Eqs. 11) are at equilibrium and we receive a
new datapoint given by ${\bf x}_{T+1},{\bf y}_{T+1}$. For $c^{(l)}=0$ and in
the limit of ${\epsilon}\equiv g_{a}/g_{b}\to 0$, the leading term in
${\epsilon}$ for the forward weight updates $\delta{\bf W}_{b}^{(l-1)}$ for
this new sample is given by

$\displaystyle\delta{\bf W}_{b}^{(l-1)}\propto{\bf v}_{a,T+1}^{(l)}{{\bf
z}_{T+1}^{(l-1)\top}}=$ $\displaystyle\;{\epsilon}^{n-l}\left[{{\bf
W}_{a}^{(l)\top}\cdots{{\bf W}_{a}^{(n-1)\top}}}({\bf y}_{T+1}-\tilde{\bf
y}_{T+1})\right]{{\bf
z}_{T+1}^{(l-1)\top}}+\mathcal{O}({{\epsilon}^{n-l+1}}),$

where $\tilde{\bf y}_{T+1}={\bf Y}{\bf Z}^{(l-1)\top}\left({\bf Z}^{(l-1)}{\bf
Z}^{(l-1)\top}\right)^{-1}{\bf z}_{T+1}^{(l-1)}$ is the optimal linear
inferred value for ${\bf y}_{T+1}$.

#### Biological interpretation

Recall that the learning rule for the basal synapses is proportional to ${\bf
v}_{a}^{(l)}$ (Eq. (14)). We argued in the previous section that this quantity
can be implemented as the calcium plateau potential. In this light, Prop. 1
shows that the calcium plateau potential is implicitly encoding the prediction
error. This is remarkable as in our algorithm, neither the prediction nor the
prediction error is explicitly computed. However, via the combined learning
dynamics of the synapses and the neural dynamics of the interneurons, the
final apical signal implicitly encodes the prediction error. This is in spirit
similar to difference target propagation (DTP) [47], however, in DTP, the
difference between the forward pass and the target is _explicitly_ computed
and backpropagated.

###  5.2 Adaptive latent dimension

For optimization within a single hidden layer, we show that the system
performs adaptive latent dimension discovery by thresholding eigenvalues of a
Gramian including contributions from both adjacent layers. If we take the loss
function in (9), setting $g_{a}=1,g_{b}=\epsilon$, and optimizing over ${\bf
W}_{a},{\bf W}_{b}$, we get the objective $\operatorname{Tr}[(c+\epsilon){\bf
Z}^{\top}{\bf Z}-{\bf Z}^{\top}{\bf Z}\big{[}\epsilon{\bf X}^{\top}({\bf
X}{\bf X}^{\top})^{-1}{\bf X}+{\bf Y}^{\top}{\bf Y}\big{]}$ for ${\bf Z}$. The
next theorem indicates the solution:

######  Theorem 1.

Let the concatenated matrix $\bm{\Xi}=[\epsilon^{1/2}{\bf C}_{x}^{-1/2}{\bf
X},{\bf Y}]\in\mathbb{R}^{(d_{X}+d_{Y})\times T}$ have SVD
$\bm{\Xi}=\sum_{\alpha=1}^{(d_{X}+d_{Y})}{\bf u}_{\alpha}\lambda_{\alpha}{\bf
v}_{\alpha}^{\top}$ with $\\{{\bf u}_{\alpha}\in\mathbb{R}^{(d_{X}+d_{Y})}\\}$
and $\\{{\bf v}_{\alpha}\in\mathbb{R}^{T}\\}$ both being sets of orthonormal
vectors, with the convention that $(\lambda_{\alpha})$ are sorted in
decreasing order. We consider the minimization

$\min_{{\bf Z}\in\mathbb{R}^{d\times T},\,\,{\bf Z}{\bf Z}^{\top}\preccurlyeq
I_{d}}\operatorname{Tr}\left[(c+\epsilon){\bf Z}^{\top}{\bf Z}-{\bf
Z}^{\top}{\bf Z}\big{[}\epsilon{\bf X}^{\top}({\bf X}{\bf X}^{\top})^{-1}{\bf
X}+{\bf Y}^{\top}{\bf Y}\big{]}\right].$ (15)

Then one of the optimal solutions is given by $\hat{\bf
Z}=\sum_{\alpha=1}^{D}{\bf
w}_{\alpha}\mathbbm{1}(\lambda_{\alpha}>c+\epsilon){\bf v}_{\alpha}^{\top}$
where $\\{{\bf w}_{\alpha}\in\mathbb{R}^{d}\\}$ is an arbitrary set of
orthonormal vectors, with $D=\min(d,d_{X}+d_{Y})$.

Thus, the latent modes are found by thresholding the eignevalues [72] of the
sum of Gramians [71] in Eq. (15). From this theorem, we derive the
computational limits of ${\epsilon}\to\infty$ and ${\epsilon}\to 0$ in the
supplementary materials. We also show that the neurons of the narrowest layer
get whitened and the other layers become low-rank. For further details of
these statements and all proofs see SM Sec. B.

#### Biological interpretation

In our model, the number of neurons in the intermediate layers is fixed.
However, the dimensionality of the representation encoded in these
intermediate layers is determined dynamically according to the spectrum of the
eigenvalue problem specified in Thm. 1. In particular, it is the quantity
$c^{(l)}$, which encodes the somatic compartment leak conductance, that
biologically sets an adaptive threshold for the dimensionality of the latent
representation. In other words, lower (resp. higher) somatic leak corresponds
to higher (resp. lower) sensitivity to the small eigenvalues of the above
eigenvalue problem. For details of the interpretation of the eigenvalue
problem in different limits (${\epsilon}\to\infty$ and ${\epsilon}\to 0$), see
SM Sec. B.

##  6 Numerical experiments

Figure 3: Our algorithm (BioCCPC) compares well to less biologically plausible
predictive-coding (PC) implementations while enforcing an inequality
constraint on the covariance of hidden-layer activations. Simulations run on
the MNIST dataset using a network with two hidden layers, sizes 50 and 5. The
shaded area represents the 95% range out of 40 runs. (a) Evolution of the
predictive-coding loss, Eq. (6), on a validation set during learning. (b), (c)
Evolution of the eigenvalues of the covariance matrix of hidden-layer
activations in the two hidden layers. The dashed lines show the scale of the
constraint; see SM Sec. A. Inset: covariance matrix at the end of training.
Note how the activity is whitened in the second hidden layer (c), and low-rank
in the first hidden layer (b).

We tested our BioCCPC algorithm on the MNIST [79] as well as other
datasets.333The code is available at https://github.com/ttesileanu/bio-pcn. We
compared our results with a standard predictive coding network [4] with linear
activation (PC) and with a version that includes a covariance constraint
(PC+cons; Sec. C). We used unit variances, $\sigma^{(l)2}=1$, in the PC
networks, and matching conductances in our algorithm,
$g_{b}^{(1)}=g_{a}^{(n)}=1$, $g_{a}^{(l)}=g^{(l)}=\tfrac{1}{2}$ for other $l$.

Our method finds solutions that are almost optimal in terms of predictive-
coding loss on a held-out validation set (see Fig. 3a) and are generally
better than the solutions obtained from simply adding the covariance
constraint to a predictive-coding network. Learning is slower in our network
than in unconstrained PCN, as it takes time for the covariance constraint to
be enforced. It is possible that different initialization schemes could
improve the learning speed, but we leave this for future work. Note that since
we focus on linear networks, metrics based on prediction accuracy are not very
informative for the datasets we consider, which are not linearly separable.

The covariance constraint is saturated in the narrowest layer—in our example,
the second hidden layer, of size 5—see Fig. 3c. The activation in the other
hidden layer is rank-restricted by the narrowest layer, as described above,
and so it is only whitened within a 5D subspace, Fig. 3b. In SM Sec. C, we
provide further numerical experiments, exploring the effects of thresholding,
the size of the network, and performance on other datasets. In all tested
scenarios, BioCCPC achieves the same validation loss as PC (Figs. 5 and 6 of
SM Sec. C). We also show that the forward and backward weights, ${\bf W}_{b}$
and ${\bf W}_{a}^{\top}$, are not symmetric, in accordance with biological
observations, and explore their relationship.

In SM Sec. D we compare PC and BioCCPC with regards to generalization. We find
that there is very little generalization gap for either case (Fig. 9). This is
due to the linearity of our model, which acts as a strong regularizer. Because
of this, the potential generalization benefits of whitening the latent
variables are not observed. Therefore, the primary benefit of the covariance
constraint in this work is the elimination of weight sharing in ${\bf W}_{a}$
(cf. Sec. 3.1).

##  7 Conclusion

We have derived a biologically plausible algorithm for covariance-constrained
predictive coding and have shown that it avoids the criticisms of PC and has
many features in common with recent experimental observations. Our algorithm
does not have symmetric forward and backward connectivity matrices, and does
not need one-to-one connectivity between any neurons. Furthermore, we showed
that the learning rules of our algorithm closely match experimental
observations, and we connected our parameters to neurophysiological
quantities. Using the simplicity afforded to us because of the linearity, we
showed a number of interesting properties of our algorithm, including the
effect of the somatic leak term, which acts as an adaptive thresholding
mechanism. We also showed how our algorithm implicitly computes a difference
term that it uses for learning. We argued that this term could be encoded in
the calcium plateau potential, which is known to affect the learning of basal
synapses. We hope that this concrete connection between our normative model
and observation can further our understanding of cortical computation by
guiding future experiments.

## Acknowledgments and Disclosure of Funding

We would like to thank Colin Bredenberg for guidance and helpful comments
throughout the project. SG, TT, YB, and DC were internally funded by the
Simons Foundation. AS was partly funded by Simons Foundation Neuroscience
grant SF 626323 during this work.

## Broader impact

We do not foresee any potential societal harm caused by our work. We hope our
contribution moves forward the discussion on the relationship between multi-
layer networks and cortical structure.

## References

  * [1] Rajesh P.N. Rao and Dana H. Ballard.  Predictive coding in the visual cortex: A functional interpretation of some extra-classical receptive-field effects.  Nature Neuroscience, 1999.
  * [2] Karl Friston.  A theory of cortical responses.  Philosophical Transactions of the Royal Society B: Biological Sciences, 2005.
  * [3] Mandyam Veerambudi Srinivasan, Simon Barry Laughlin, A. Dubs, and George Adrian Horridge.  Predictive coding: a fresh view of inhibition in the retina.  Proceedings of the Royal Society of London. Series B. Biological Sciences, 216(1205):427–459, 1982.
  * [4] James C.R. Whittington and Rafal Bogacz.  An approximation of the error backpropagation algorithm in a predictive coding network with local hebbian synaptic plasticity.  Neural Computation, 29(5):1229–1262, may 2017.
  * [5] Tommaso Salvatori, Luca Pinchetti, Beren Millidge, Yuhang Song, Tianyi Bao, Rafal Bogacz, and Thomas Lukasiewicz.  Learning on arbitrary graph topologies via predictive coding, 2022.
  * [6] André Ofner and Sebastian Stober.  Differentiable generalised predictive coding, 2021.
  * [7] Martin Boerlin, Christian K Machens, and Sophie Denève.  Predictive coding of dynamical variables in balanced spiking networks.  PLoS computational biology, 9(11):e1003258, 2013.
  * [8] Jonathan Kadmon, Jonathan Timcheck, and Surya Ganguli.  Predictive coding in balanced neural networks with noise, chaos and delays.  Advances in neural information processing systems, 33:16677–16688, 2020.
  * [9] Ryan Smith, Paul Badcock, and Karl J Friston.  Recent advances in the application of predictive coding and active inference models within clinical neuroscience.  Psychiatry and Clinical Neurosciences, 75(1):3–13, 2021.
  * [10] Paul F. Kinghorn, Beren Millidge, and Christopher L. Buckley.  Habitual and reflective control in hierarchical predictive coding, 2021\.
  * [11] Tommaso Salvatori, Yuhang Song, Yujian Hong, Lei Sha, Simon Frieder, Zhenghua Xu, Rafal Bogacz, and Thomas Lukasiewicz.  Associative Memories via Predictive Coding.  In M Ranzato, A Beygelzimer, Y Dauphin, P S Liang, and J Wortman Vaughan, editors, Advances in Neural Information Processing Systems, volume 34, pages 3874–3886. Curran Associates, Inc., 2021.
  * [12] Beren Millidge, Anil Seth, and Christopher L Buckley.  Predictive coding: a theoretical and experimental review, 2021.
  * [13] Beren Millidge, Tommaso Salvatori, Yuhang Song, Rafal Bogacz, and Thomas Lukasiewicz.  Predictive coding: Towards a future of deep learning beyond backpropagation?, 2022.
  * [14] André M. Bastos, W. Martin Usrey, Rick A Adams, George R. Mangun, Pascal Fries, and Karl J. Friston.  Canonical microcircuits for predictive coding.  Neuron, 76:695–711, 2012.
  * [15] Micha Heilbron and Maria Chait.  Great Expectations: Is there Evidence for Predictive Coding in Auditory Cortex?  Neuroscience, 389:54–73, 2018.
  * [16] Naoki Kogo and Chris Trengove.  Is predictive coding theory articulated enough to be testable?  Frontiers in Computational Neuroscience, 9, 2015.
  * [17] Nelson Spruston.  Pyramidal neurons: dendritic structure and synaptic integration.  Nature Reviews Neuroscience, 9(3):206–221, 2008.
  * [18] Matthew Larkum.  A cellular mechanism for cortical associations: an organizing principle for the cerebral cortex.  Trends in Neurosciences, 36(3):141–151, 2013.
  * [19] Per Jesper Sjöström and Michael Häusser.  A cooperative switch determines the sign of synaptic plasticity in distal dendrites of neocortical pyramidal neurons.  Neuron, 51(2):227–238, 2006.
  * [20] Nace L Golding, Nathan P Staff, and Nelson Spruston.  Dendritic spikes as a mechanism for cooperative long-term potentiation.  Nature, 418(6895):326–331, 2002.
  * [21] Frédéric Gambino, Stéphane Pagès, Vassilis Kehayas, Daniela Baptista, Roberta Tatti, Alan Carleton, and Anthony Holtmaat.  Sensory-evoked LTP driven by dendritic plateau potentials in vivo.  Nature, 515(7525):116–119, 2014.
  * [22] Guy Major, Matthew E. Larkum, and Jackie Schiller.  Active properties of neocortical pyramidal neuron dendrites.  Annual Review of Neuroscience, 36(1):1–24, 2013.  PMID: 23841837.
  * [23] Katie C Bittner, Christine Grienberger, Sachin P Vaidya, Aaron D Milstein, John J Macklin, Junghyup Suh, Susumu Tonegawa, and Jeffrey C Magee.  Conjunctive input processing drives feature selectivity in hippocampal ca1 neurons.  Nature neuroscience, 18(8):1133, 2015.
  * [24] Katie C Bittner, Aaron D Milstein, Christine Grienberger, Sandro Romani, and Jeffrey C Magee.  Behavioral time scale synaptic plasticity underlies ca1 place fields.  Science, 357(6355):1033–1036, 2017.
  * [25] Jeffrey C. Magee and Christine Grienberger.  Synaptic plasticity forms and functions.  Annual Review of Neuroscience, 43, 2020.
  * [26] Jason Hardie and Nelson Spruston.  Synaptic depolarization is more effective than back-propagating action potentials during induction of associative long-term potentiation in hippocampal pyramidal neurons.  Journal of Neuroscience, 29(10):3233–3241, 2009.
  * [27] Konrad P Körding and Peter König.  Supervised and unsupervised learning with two sites of synaptic integration.  Journal of Computational Neuroscience, 11(3):207–215, 2001.
  * [28] Robert Urbanczik and Walter Senn.  Learning by the Dendritic Prediction of Somatic Spiking.  Neuron, 81(3):521–528, 2014.
  * [29] Jordan Guergiuev, Timothy P Lillicrap, and Blake A Richards.  Deep learning with segregated dendrites.  arXiv preprint arXiv:1610.00161, 2016.
  * [30] João Sacramento, Rui Ponte Costa, Yoshua Bengio, and Walter Senn.  Dendritic cortical microcircuits approximate the backpropagation algorithm.  Advances in Neural Information Processing Systems, 2018-Decem(Nips):8721–8732, oct 2018.
  * [31] Alexandre Payeur, Jordan Guerguiev, Friedemann Zenke, Blake A Richards, and Richard Naud.  Burst-dependent synaptic plasticity can coordinate learning in hierarchical circuits.  bioRxiv, 2020.
  * [32] Albert Gidon, Timothy Adam Zolnik, Pawel Fidzinski, Felix Bolduan, Athanasia Papoutsi, Panayiota Poirazi, Martin Holtkamp, Imre Vida, and Matthew Evan Larkum.  Dendritic action potentials and computation in human layer 2/3 cortical neurons.  Science, 367(6473):83–87, jan 2020.
  * [33] Aaron D. Milstein, Yiding Li, Katie C. Bittner, Christine Grienberger, Ivan Soltesz, Jeffrey C. Magee, and Sandro Romani.  Bidirectional synaptic plasticity rapidly modifies hippocampal representations independent of correlated activity.  bioRxiv, 2020.
  * [34] David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams.  Learning internal representations by error propagation.  In David E. Rumelhart and James L. Mcclelland, editors, Parallel Distributed Processing: Explorations in the Microstructure of Cognition, Volume 1: Foundations, pages 318–362. MIT Press, Cambridge, MA, 1986.
  * [35] Ian Goodfellow, Yoshua Bengio, and Aaron Courville.  Deep learning.  MIT press, 2016.
  * [36] Seungjin Choi, Andrzej Cichocki, Hyung-Min Park, and Soo-Young Lee.  Blind source separation and independent component analysis: A review.  Neural Information Processing-Letters and Reviews, 6(1):1–57, 2005\.
  * [37] Cengiz Pehlevan and Dmitri B Chklovskii.  Neuroscience-inspired online unsupervised learning algorithms: Artificial neural networks.  IEEE Signal Processing Magazine, 36(6):88–96, 2019.
  * [38] Adrian A Wanner and Rainer W Friedrich.  Whitening of odor representations by the wiring diagram of the olfactory bulb.  Nature neuroscience, 23(3):433–442, 2020.
  * [39] Guillaume Desjardins, Karen Simonyan, Razvan Pascanu, and Koray Kavukcuoglu.  Natural neural networks, 2015.
  * [40] Yuki Fujimoto and Toru Ohira.  A neural network model with bidirectional whitening.  In Artificial Intelligence and Soft Computing, pages 47–57. Springer International Publishing, 2018.
  * [41] Ping Luo.  Learning deep architectures via generalized whitened neural networks.  In International Conference on Machine Learning, pages 2238–2246. PMLR, 2017.
  * [42] P. J. Werbos.  Beyond Regression: New Tools for Prediction and Analysis in the Behavioral Sciences.  PhD thesis, Harvard University, 1974.
  * [43] Francis Crick.  The recent excitement about neural networks.  Nature, 337(6203):129–132, 1989.
  * [44] Yann Le Cun.  Learning Process in an Asymmetric Threshold Network.  In Disordered Systems and Biological Organization, pages 233–240. Springer, 1986.
  * [45] Y. Le Cun and Franccoise Fogelman-Soulié.  Modèles connexionnistes de l’apprentissage.  Intellectica. Revue de l’Association pour la Recherche Cognitive, 1987.
  * [46] John F. Kolen and Jordan B. Pollack.  Backpropagation without weight transport.  Proceedings of 1994 IEEE International Conference on Neural Networks (ICNN’94), 3:1375–1380 vol.3, 1994.
  * [47] Dong-Hyun Lee, Saizheng Zhang, Antoine Biard, and Yoshua Bengio.  Target Propagation.  Iclr 15, 1(3):1–16, 2014.
  * [48] Timothy P. Lillicrap, Daniel Cownden, Douglas B. Tweed, and Colin J. Akerman.  Random synaptic feedback weights support error backpropagation for deep learning.  Nature Communications, 7:1–10, 2016.
  * [49] Benjamin Scellier and Yoshua Bengio.  Equilibrium propagation: Bridging the gap between energy-based models and backpropagation.  Frontiers in Computational Neuroscience, 11(May), 2017.
  * [50] Dmitry Krotov and John J. Hopfield.  Unsupervised learning by competing hidden units.  Proceedings of the National Academy of Sciences, 116(16):7723–7731, 2019.
  * [51] Mohamed Akrout, Collin Wilson, Peter C. Humphreys, Timothy Lillicrap, and Douglas Tweed.  Deep Learning without Weight Transport.  arXiv, apr 2019.
  * [52] Jordan Guerguiev, Konrad Kording, and Blake Richards.  Spike-based causal inference for weight alignment.  arXiv preprint arXiv:1910.01689, 10 2019.
  * [53] Maxence Ernoult, Julie Grollier, Damien Querlioz, Yoshua Bengio, and Benjamin Scellier.  Equilibrium Propagation with Continual Weight Updates.  arXiv preprint arXiv:2005.04168, pages 1–37, 2020.
  * [54] Will Xiao, Honglin Chen, Qianli Liao, and Tomaso A. Poggio.  Biologically-plausible learning algorithms can scale to large datasets.  CoRR, abs/1811.03567, 2018.
  * [55] Siavash Golkar, David Lipshutz, Yanis Bahroun, Anirvan Sengupta, and Dmitri Chklovskii.  A simple normative network approximates local non-Hebbian learning in the cortex.  In H Larochelle, M Ranzato, R Hadsell, M F Balcan, and H Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 7283–7295. Curran Associates, Inc., 2020.
  * [56] Siavash Golkar, David Lipshutz, Yanis Bahroun, Anirvan M. Sengupta, and Dmitri B. Chklovskii.  A biologically plausible neural network for local supervision in cortical microcircuits, 2020.
  * [57] Yoshua Bengio.  Deriving differential target propagation from iterating approximate inverses, 2020.
  * [58] Jyotibdha Acharya, Arindam Basu, Robert Legenstein, Thomas Limbacher, Panayiota Poirazi, and Xundong Wu.  Dendritic Computing: Branching Deeper into Machine Learning.  Neuroscience, 489:275–289, 2022.
  * [59] Timothy P Lillicrap, Adam Santoro, Luke Marris, Colin J Akerman, and Geoffrey Hinton.  Backpropagation and the brain.  Nature Reviews Neuroscience, 21(6):335–346, 2020.
  * [60] João Sacramento, Rui Ponte Costa, Yoshua Bengio, and Walter Senn.  Dendritic error backpropagation in deep cortical microcircuits.  arXiv preprint arXiv:1801.00062, pages 1–37, 2017.
  * [61] Beren Millidge, Alexander Tschantz, and Christopher L. Buckley.  Predictive coding approximates backprop along arbitrary computation graphs, 2020.
  * [62] Yuhang Song, Thomas Lukasiewicz, Zhenghua Xu, and Rafal Bogacz.  Can the Brain Do Backpropagation? — Exact Implementation of Backpropagation in Predictive Coding Networks.  In H Larochelle, M Ranzato, R Hadsell, M F Balcan, and H Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 22566–22579. Curran Associates, Inc., 2020.
  * [63] Tommaso Salvatori, Yuhang Song, Thomas Lukasiewicz, Rafal Bogacz, and Zhenghua Xu.  Predictive coding can do exact backpropagation on convolutional and recurrent neural networks.  CoRR, abs/2103.03725, 2021.
  * [64] Robert Rosenbaum.  On the relationship between predictive coding and backpropagation, 2021\.
  * [65] Matthew E. Larkum, Thomas Nevian, Maya Sandler, Alon Polsky, and Jackie Schiller.  Synaptic integration in tuft dendrites of layer 5 pyramidal neurons: A new unifying principle.  Science, 325(5941):756–760, 2009.
  * [66] Hiroto Takahashi and Jeffrey C Magee.  Pathway interactions and synaptic plasticity in the dendritic tuft regions of CA1 pyramidal neurons.  Neuron, 62(1):102–111, 2009.
  * [67] Thomas Klausberger, Peter J Magill, László F Márton, J David B Roberts, Philip M Cobden, György Buzsáki, and Peter Somogyi.  Brain-state-and cell-type-specific firing of hippocampal interneurons in vivo.  Nature, 421(6925):844–848, 2003.
  * [68] Therese Riedemann.  Diversity and function of somatostatin-expressing interneurons in the cerebral cortex.  International Journal of Molecular Sciences, 20(12):2952, 2019.
  * [69] Karl Friston.  The free-energy principle: A unified brain theory?, 2010.
  * [70] Ryota Kanai, Yutaka Komura, Stewart Shipp, and Karl Friston.  Cerebral hierarchies: predictive processing, precision and the pulvinar.  Philosophical Transactions of the Royal Society B: Biological Sciences, 370(1668):20140169, 2015.
  * [71] David Lipshutz, Yanis Bahroun, Siavash Golkar, Anirvan M. Sengupta, and Dmitri B. Chklovskii.  A Biologically Plausible Neural Network for Multichannel Canonical Correlation Analysis.  Neural Computation, 33(9):2309–2352, 08 2021.
  * [72] Cengiz Pehlevan and Dmitri Chklovskii.  A Normative Theory of Adaptive Dimensionality Reduction in Neural Networks.  In C Cortes, N Lawrence, D Lee, M Sugiyama, and R Garnett, editors, Advances in Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015.
  * [73] Charles D Gilbert and Wu Li.  Top-down influences on visual processing.  Nature Reviews Neuroscience, 14(5):350–363, 2013.
  * [74] Georg B Keller and Thomas D Mrsic-Flogel.  Predictive processing: a canonical cortical computation.  Neuron, 100(2):424–435, 2018.
  * [75] Björn M. Kampa, Johannes J. Letzkus, and Greg J. Stuart.  Dendritic mechanisms controlling spike-timing-dependent synaptic plasticity.  Trends in Neurosciences, 30(9):456–463, 2007.
  * [76] Andrew M Saxe, James L McClelland, and Surya Ganguli.  Exact solutions to the nonlinear dynamics of learning in deep linear neural networks.  arXiv preprint arXiv:1312.6120, 2013.
  * [77] Sanjeev Arora, Nadav Cohen, and Elad Hazan.  On the optimization of deep networks: Implicit acceleration by overparameterization.  In International Conference on Machine Learning, pages 244–253. PMLR, 2018.
  * [78] Simon S Du, Wei Hu, and Jason D Lee.  Algorithmic regularization in learning deep homogeneous models: Layers are automatically balanced.  Advances in Neural Information Processing Systems, 31, 2018.
  * [79] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.

Supplementary Materials

This is the supplementary materials section for the NeurIPS 2022 submission
titled “Constrained Predictive Coding as a Biologically Plausible Model of the
Cortical Hierarchy”.

##  Appendix A Derivation of upper bound

In this section we give a more detailed derivation of our framework, using a
slightly modified order of steps compared to the main text. This should help
clarify some aspects of the final algorithm.

Similar to the main text, we start with the predictive coding objective and
introduce the weight doubles. We arrive at

$\min_{{\bf Z},{\bf W}}L=\min_{{\bf Z},{\bf W}_{a},{\bf
W}_{b}}\frac{1}{2}\sum_{l=1}^{n-1}\left[g_{b}^{(l)}\bigl{\|}{\bf Z}^{(l)}-{\bf
W}_{b}^{(l-1)}{\bf Z}^{(l-1)}\bigr{\|}_{F}^{2}\right.\\\
\left.+\,g_{a}^{(l+1)}\bigl{\|}{\bf Z}^{(l+1)}-{\bf W}_{a}^{(l)}{\bf
Z}^{(l)}\bigr{\|}_{F}^{2}+c^{(l)}\bigl{\|}{\bf Z}^{(l)}{\bf
Z}^{(l)\top}\bigr{\|}^{2}_{F}\right],$ (16)

where we have also introduced a quadratic cost (a prior on ${\bf Z}$) with
coefficient $c^{(l)}$. We now impose the constraint on the empirical
covariance, $\frac{1}{T}{\bf Z}^{(l)}{\bf Z}^{(l)\top}\preceq\rho^{(l)}{\bf
I}$, of the internal variables ${\bf Z}^{(l)}$444The factors $\rho^{(l)}$
which set the scale of the constraint were set to 1 in the derivation in the
main text for brevity. However, it is sensible for this scale factor to vary
from layer to layer and possibly depend on the layer width. In Sec. C, we
explore the effects of this parameter on the performance of the network.:

$\min_{{\bf Z},{\bf W}}L\leq\min_{\begin{subarray}{c}{\bf Z},{\bf W}_{a},{\bf
W}_{b}\\\ \tfrac{1}{T}{\bf Z}^{(l)}{\bf Z}^{(l)\top}\preceq\rho^{(l)}{\bf
I}\end{subarray}}\frac{1}{2}\sum_{l=1}^{n-1}\left[g_{b}^{(l)}\bigl{\|}{\bf
Z}^{(l)}-{\bf W}_{b}^{(l-1)}{\bf Z}^{(l-1)}\bigr{\|}_{F}^{2}\right.\\\
\left.+\,g_{a}^{(l+1)}\bigl{\|}{\bf Z}^{(l+1)}-{\bf W}_{a}^{(l)}{\bf
Z}^{(l)}\bigr{\|}_{F}^{2}+c^{(l)}\bigl{\|}{\bf Z}^{(l)}{\bf
Z}^{(l)\top}\bigr{\|}^{2}_{F}\right].$ (17)

Expanding this and using the constraint we again upper bound the objective:

$\min_{{\bf Z},{\bf W}}L\leq\min_{\begin{subarray}{c}{\bf Z},{\bf W}_{a},{\bf
W}_{b}\\\ \tfrac{1}{T}{\bf Z}^{(l)}{\bf Z}^{(l)\top}\preceq\rho^{(l)}{\bf
I}\end{subarray}}\frac{1}{2}\sum_{l=1}^{n-1}\left[g_{b}^{(l)}\bigl{\|}{\bf
Z}^{(l)}-{\bf W}_{b}^{(l-1)}{\bf
Z}^{(l-1)}\bigr{\|}_{F}^{2}+g_{a}^{(l+1)}\bigl{\|}{\bf Z}^{(l+1)}{\bf
Z}^{(l+1)\top}\bigr{\|}^{2}_{F}\right.\\\
\left.+\,g_{a}^{(l+1)}\operatorname{Tr}\bigl{(}-2{\bf Z}^{(l+1)\top}{\bf
W}_{a}^{(l)}{\bf Z}^{(l)}+T\rho^{(l)}{\bf W}_{a}^{(l)\top}{\bf
W}_{a}^{(l)}\bigr{)}+c^{(l)}\bigl{\|}{\bf Z}^{(l)}{\bf
Z}^{(l)\top}\bigr{\|}^{2}_{F}\right],$ (18)

We next implement the inequality constraint using positive-definite Lagrange
multipliers:

$\displaystyle\min_{{\bf Z},{\bf W}}L\leq$ $\displaystyle\min_{{\bf Z},{\bf
W}_{a},{\bf W}_{b}}\max_{\bf Q}\sum_{l=1}^{n-1}\tilde{L}^{(l)},$ (19)

where

$\begin{split}\tilde{L}^{(l)}&=\frac{1}{2}\left[g_{b}^{(l)}\bigl{\|}{\bf
Z}^{(l)}-{\bf W}_{b}^{(l-1)}{\bf
Z}^{(l-1)}\bigr{\|}_{F}^{2}+c_{q}^{(l)}\operatorname{Tr}{\bf Q}^{(l)\top}{{\bf
Q}^{(l)}}\bigl{(}{{\bf Z}^{(l)}}{\bf Z}^{(l)\top}-T\rho^{(l)}{\bf
I}\bigr{)}\right.\\\
&\qquad\qquad\left.+\,g_{a}^{(l+1)}\operatorname{Tr}\bigl{(}-2{\bf
Z}^{(l+1)\top}{\bf W}_{a}^{(l)}{\bf Z}^{(l)}+T\rho^{(l)}{\bf
W}_{a}^{(l)\top}{\bf W}_{a}^{(l)}\bigr{)}\right.\\\
&\qquad\qquad\left.+\,c^{(l)}\bigl{\|}{\bf Z}^{(l)}{\bf
Z}^{(l)\top}\bigr{\|}^{2}_{F}+g_{a}^{(l+1)}\bigl{\|}{\bf Z}^{(l+1)}{\bf
Z}^{(l+1)\top}\bigr{\|}^{2}_{F}\right],\end{split}$ (20)

where we have added a positive constant $c_{q}^{(l)}$ defining the coefficient
of the inequality constraint. This constant does not affect the solution of
the optimization problem but does affect the details of the neural dynamics.
Motivated by biological arguments555Since the feedback of the interneurons is
received and integrated in the apical compartment, we expect the apical-to-
soma conductance to also be present here., we set this constant equal to
$g_{a}^{(l)}$. In doing so, the interneuron variables no longer need a
$(1/g_{a}^{(l)})$ in their definition, as in the main text, and are instead
given by ${\bf n}^{(l)}={\bf Q}^{(l)}{\bf z}^{(l)}$.

A final upper bound to (19) can be found if each $\tilde{L}^{(l)}$ is
optimized individually. Explicitly, we perform gradient descent of each
variable, ${\bf Z}^{(l)},{\bf W}_{a}^{(l)},{\bf W}_{b}^{(l)}$, only with
respect to $\tilde{L}^{(l)}$ and not with respect to the entire objective $L$:

$\begin{split}\delta{\bf Z}^{(l)}=\nabla_{{\bf Z}^{(l)}}L&\;\to\;\nabla_{{\bf
Z}^{(l)}}\tilde{L}^{(l)},\\\ \delta{\bf W}_{a}^{(l)}=\nabla_{{\bf
W}_{a}^{(l)}}L&\;\to\;\nabla_{{\bf W}_{a}^{(l)}}\tilde{L}^{(l)},\\\ \delta{\bf
W}_{b}^{(l)}=\nabla_{{\bf W}_{b}^{(l)}}L&\;\to\;\nabla_{{\bf
W}_{b}^{(l)}}\tilde{L}^{(l)}.\end{split}$ (21)

This affords us an upper bound because we are no longer finding the actual
minimum of $L$. Symbolically this can be written as in (9):

$\min_{{\bf Z},{\bf W}}L\leq\sum_{l=1}^{n-1}\min_{{\bf Z},{\bf W}_{a},{\bf
W}_{b}}\max_{\bf Q}\tilde{L}^{(l)},$ (22)

where the precise definition of this equation is given by the update equations
(21). Note that these updates are no longer curl-free and hence cannot be
found as the gradient of a single objective function.

##  Appendix B Theoretical section proofs

Note that the definition of ${\epsilon}$ in Proposition 1 is different from
the definition in other sections.

###  B.1 Proof of Proposition 1

To prove Proposition 1, we first show that the latent variables ${\bf
z}^{(l)}$ have a specific form. As in the main text, we define
${\epsilon}^{(l)}=g_{a}^{(l)}/g_{b}^{(l)}$ and, for simplicity, we assume
${\epsilon}^{(l)}={\epsilon}$ is the same for all layers. For simplicity also,
we assume $c^{(l)}$ (equivalently $d^{(l)}$ in previous section) is equal to
zero, and use $c_{q}^{(l)}=g_{a}^{(l)}$ as in the previous section. This
choice makes the notation simpler without changing our results. Equation (10)
becomes

$\displaystyle\dot{\bf z}^{(l)}=g_{b}^{(l)}{\bf W}_{b}^{(l-1)}{\bf
z}^{(l-1)}+g_{a}^{(l)}{\bf W}_{a}^{(l)\top}{\bf z}^{(l+1)}-g_{b}^{(l)}{\bf
z}^{(l)}-g_{a}^{(l)}{\bf Q}^{(l)\top}{\bf Q}^{(l)}{\bf z}^{(l)}.$ (23)

######  Proposition 2.

Assume the neural dynamics given in Eq. (23) are at equilibrium. Then the
latent variables take the form

${\bf z}^{(n-k)}={\bf B}_{k}\left({\bf W}_{b}^{(n-k-1)}{\bf
z}^{(n-k-1)}+{\epsilon}^{k}{\bf A}_{k}{{\bf W}_{a}^{(n-1)\top}}{\bf
y}\right),$ (24)

where ${\bf B}_{k}$ and ${\bf A}_{k}$ are independent of $y$ and defined
recursively via

$\displaystyle{\bf A}_{k+1}$ $\displaystyle={{\bf W}_{a}^{(n-k-1)\top}}{\bf
B}_{k}{\bf A}_{k},$ $\displaystyle{\bf B}_{k+1}$
$\displaystyle=\left(1-{\epsilon}{\bf F}^{(n-k-1)}{{\bf
W}_{a}^{(n-k-1)\top}}{\bf B}_{k}{\bf W}_{b}^{(n-k-1)}\right)^{-1}{\bf
F}^{(n-k-1)},$

with ${\bf F}^{(l)}=\left({\bf I}+{\epsilon}{{\bf Q}^{(l)\top}}{\bf
Q}^{(l)}\right)^{-1}$ and

${\bf B}_{1}={\bf F}^{(n-1)}\;,\quad{\bf A}_{1}=1.$

###### Proof.

We prove this via induction. Starting from the neural dynamics Eq. (23), at
equilibrium ($\dot{\bf z}^{(l)}=0$) we have:

$\displaystyle{\bf z}^{(l)}={\bf F}^{(l)}\left({\bf W}_{b}^{(l-1)}{\bf
z}^{(l-1)}+{\epsilon}{{\bf W}_{a}^{(l)\top}}{\bf z}^{(l+1)}\right).$

Setting $l=n-1$ we get:

$\displaystyle{\bf z}^{(n-1)}={\bf F}^{(n-1)}\left({\bf W}_{b}^{(n-2)}{\bf
z}^{(n-2)}+{\epsilon}{{\bf W}_{a}^{(n-1)\top}}{\bf y}\right),$

which satisfies the proposed form in Eq. (24) for $k=1$. Now, we assume that
the form in Eq. (24) holds for $z^{(n-k)}$, and show that it also holds for
$z^{(n-k-1)}$. We have

$\displaystyle{\bf z}^{(n-k-1)}=$ $\displaystyle{\bf F}^{(n-k-1)}\left({\bf
W}_{b}^{(n-k-2)}{\bf z}^{(n-k-2)}+{\epsilon}{{\bf W}_{a}^{(n-k-1)\top}}{\bf
z}^{(n-k)}\right)$ $\displaystyle=$ $\displaystyle{\bf F}^{(n-k-1)}\left({\bf
W}_{b}^{(n-k-2)}{\bf z}^{(n-k-2)}+{\epsilon}{{\bf W}_{a}^{(n-k-1)\top}}{\bf
B}_{k}\left({\bf W}_{b}^{(n-k-1)}{\bf z}^{(n-k-1)}+{\epsilon}^{k}{\bf
A}_{k}{{\bf W}_{a}^{(n-1)\top}}{\bf y}\right)\right).$

Collecting ${\bf z}^{(n-k-1)}$ and multiplying by the inverse on the left we
arrive at:

$\displaystyle{\bf z}^{(n-k-1)}=$ $\displaystyle\left(1-{\epsilon}{\bf
F}^{(n-k-1)}{{\bf W}_{a}^{(n-k-1)\top}}{\bf B}_{k}{\bf
W}_{b}^{(n-k-1)}\right)^{-1}{\bf F}^{(n-k-1)}\times$ $\displaystyle\hskip
115.0pt\times\left({\bf W}_{b}^{(n-k-2)}{\bf
z}^{(n-k-2)}+{\epsilon}^{k+1}{{\bf W}_{a}^{(n-k-1)\top}}{\bf B}_{k}{\bf
A}_{k}{{\bf W}_{a}^{(n-1)\top}}{\bf y}\right),$

as proposed. ∎

Note that, whereas Proposition 1 holds only in the ${\epsilon}\to 0$ limit,
Proposition 2 holds for general ${\epsilon}$. Using this result we show that
the relationship between the activity, feedback and feedforward inputs of the
$l^{\text{th}}$ module has a simple form given in the following corollary:

######  Corollary 1.

The leading order in ${\epsilon}$ of the dependence of ${\bf z}^{(l)}$ and
${\bf z}^{(l+1)}$ on ${\bf z}^{(l-1)}$ and ${\bf y}$, which has non-zero
dependence on ${\bf y}$, is given by

$\displaystyle{\bf z}^{(l)}$ $\displaystyle={\bf C}_{1}^{(l)}{\bf
z}^{(l-1)}+\mathcal{O}({\epsilon}^{n-l}),$ $\displaystyle{\bf z}^{(l+1)}$
$\displaystyle={\bf C}_{2}^{(l)}{\bf z}^{(l-1)}+{\epsilon}^{n-l-1}{{\bf
W}_{a}^{(l+1)\top}}\cdots{{\bf W}_{a}^{(n-1)\top}}{\bf
y}+\mathcal{O}({\epsilon}^{n-l}),$

where ${\bf C}_{1}^{(l)}$ and ${\bf C}_{2}^{(l)}$ are respectively the
expansions of ${\bf B}_{n-l}{\bf W}_{b}^{(l-1)}$ and ${\bf B}_{n-l-1}{\bf
W}_{b}^{(l)}{\bf B}_{n-l}{\bf W}_{b}^{(l-1)}$ up to order
${\epsilon}^{n-l-1}$.

###### Proof.

The result follows directly from applying Eq. (24) twice and then dropping all
but the leading dependence on ${\epsilon}$ in the coefficient of the ${\bf y}$
term. ∎

Now, we restate and prove Proposition 1, which can be viewed as a
generalization of the results in [55]. Note that in this proposition, the
exact form of ${\bf C}_{1}^{(l)}$ and ${\bf C}_{2}^{(l)}$ is not important.
All that matters is that both ${\bf z}^{(l)}$ and ${\bf z}^{(l+1)}$ have a
linear dependence on ${\bf z}^{(l-1)}$. Since the dependence of ${\bf
z}^{(l)}$ on ${\bf y}$ is subleading compared to the dependence of ${\bf
z}^{(l+1)}$ on ${\bf y}$, this can be dropped in the limit of ${\epsilon}\to
0$.

Proposition 1. Assume that the learning rules (Eqs. 11) are at equilibrium and
we receive a new datapoint given by ${\bf x}_{T+1},{\bf y}_{T+1}$. For
$c^{(l)}=0$ and in the limit of ${\epsilon}\equiv g_{a}/g_{b}\to 0$, the
leading term in ${\epsilon}$ for the forward weight updates $\delta{\bf
W}_{b}^{(l-1)}$ for this new sample is given by

$\displaystyle\delta{\bf W}_{b}^{(l-1)}\propto{\bf v}_{a,T+1}^{(l)}{{\bf
z}_{T+1}^{(l-1)\top}}=$ $\displaystyle\;{\epsilon}^{n-l}\left[{{\bf
W}_{a}^{(l)\top}\cdots{{\bf W}_{a}^{(n-1)\top}}}({\bf y}_{T+1}-\tilde{\bf
y}_{T+1})\right]{{\bf
z}_{T+1}^{(l-1)\top}}+\mathcal{O}({{\epsilon}^{n-l+1}}),$

where $\tilde{\bf y}_{T+1}={\bf Y}{\bf Z}^{(l-1)\top}\left({\bf Z}^{(l-1)}{\bf
Z}^{(l-1)\top}\right)^{-1}{\bf z}^{(l-1)}_{T+1}$ is the optimal linear
inferred value for ${\bf y}_{T}$.

###### Proof.

The statement follows from two applications of the update rule for ${\bf
W}_{b}$ given in Eq. (11a) for ${\bf W}_{b}^{(l-1)}$. Putting $c^{(l)}=0$ we
have:

$\delta{\bf W}_{b}^{(l-1)}\propto{\epsilon}\bigl{(}{\bf W}_{a}^{(l)\top}{\bf
z}^{(l+1)}-{\bf Q}^{(l)\top}{\bf Q}^{(l)}{\bf z}^{(l)}\bigr{)}{{\bf
z}^{(l-1)\top}}.$ (25)

The equilibrium condition given by $\delta{\bf W}_{b}=0$ is given by averaging
over the entire dataset:

${\bf Q}^{(l)\top}{\bf Q}^{(l)}{\bf Z}^{(l)}{\bf Z}^{(l-1)\top}={\bf
W}_{a}^{(l)\top}{\bf Z}^{(l+1)}{\bf Z}^{(l-1)\top}.$ (26)

Now we plug the derived values in Corollary 1. After right-multiplying by
$({\bf Z}^{(l-1)}{\bf Z}^{(l-1)\top})^{-1}$ we get:

${\bf Q}^{(l)\top}{\bf Q}^{(l)}{\bf C}_{1}={\bf W}_{a}^{(l)\top}{\bf
C}_{2}+{\epsilon}^{n-l-1}{\bf W}_{a}^{(l)\top}\cdots{{\bf
W}_{a}^{(n-1)\top}}{\bf Y}{\bf Z}^{(l-1)\top}\big{(}{\bf Z}^{(l-1)}{\bf
Z}^{(l-1)\top}\big{)}^{-1}.$ (27)

Now, we look at the update for a new sample $T+1$. Plugging the above as well
as the equations from Corollary 1 into the update rule we see that the ${\bf
C}_{2}^{(l)}$ dependent term cancels and we are only left with the term
depending on ${\bf y}$:

$\displaystyle\delta{\bf W}_{b}^{(l-1)}\propto$
$\displaystyle{\epsilon}^{n-l}\bigl{(}{\bf W}_{a}^{(l)\top}{\bf
z}_{T+1}^{(l+1)}-{\bf Q}^{(l)\top}{\bf Q}^{(l)}{\bf
z}_{T+1}^{(l)}\bigr{)}{{\bf z}_{T+1}^{(l-1)\top}}$ $\displaystyle=$
$\displaystyle\;{\epsilon}^{n-l}\left[{{\bf W}_{a}^{(l)\top}\cdots{{\bf
W}_{a}^{(n-1)\top}}}({\bf y}_{T+1}-\tilde{\bf y}_{T+1})\right]{{\bf
z}^{(l-1)\top}},$

with $\tilde{\bf y}_{T+1}$ as defined in the proposition. ∎

###  B.2 Proof of Theorem 1 and consequences

In this section we look at some properties of the within-layer computations of
CCPC. We will focus on the computations within the $l^{\text{th}}$ layer and
for clarity of notation, for this subsection only, we rename ${\bf
Z}^{(l)}\to{\bf Z}$, ${\bf Z}^{(l+1)}\to{\bf Y}$, and ${\bf Z}^{(l-1)}\to{\bf
X}$. We will assume $\rho^{(l)}=1$ for brevity. The results have straight-
forward generalizations for the case where $\rho^{(l)}\neq 1$.

Theorem 1. Let the concatenated matrix $\bm{\Xi}={\bf U}\bm{\Lambda}{\bf
V}^{\top}=[{\epsilon}^{1/2}{\bf C}_{x}^{-1/2}{\bf X},{\bf
Y}]\in\mathbb{R}^{(d_{X}+d_{Y})\times T}$ have singular value decomposition
$\bm{\Xi}=\sum_{\alpha=1}^{(d_{X}+d_{Y})}{\bf u}_{\alpha}\lambda_{\alpha}{\bf
v}_{\alpha}^{\top}$, with $\\{{\bf
u}_{\alpha}\in\mathbb{R}^{(d_{X}+d_{Y})}\\}$ and $\\{{\bf
v}_{\alpha}\in\mathbb{R}^{T}\\}$ both being sets of orthonormal vectors, and
singular values $(\lambda_{\alpha})$ sorted in decreasing order. We consider
the minimization

$\min_{\begin{subarray}{c}{\bf Z}\in\mathbb{R}^{d\times N}\\\ {\bf Z}{\bf
Z}^{\top}\preccurlyeq
I_{d}\end{subarray}}\operatorname{Tr}\left[(c+{\epsilon}){\bf Z}^{\top}{\bf
Z}-{\bf Z}^{\top}{\bf Z}\big{[}{\epsilon}{\bf X}^{\top}({\bf X}{\bf
X}^{\top})^{-1}{\bf X}+{\bf Y}^{\top}{\bf Y}\big{]}\right].$ (28)

Then one of the optimal solutions is given by $\hat{\bf
Z}=\sum_{\alpha=1}^{D}{\bf
w}_{\alpha}\mathbbm{1}(\lambda_{\alpha}>c+{\epsilon}){\bf v}_{\alpha}^{\top}$
where $\\{{\bf w}_{\alpha}\in\mathbb{R}^{d}\\}$ is an arbitrary set of
orthonormal vectors, with $D=\min(d,d_{X}+d_{Y})$.

###### Proof.

Let the SVD of ${\bf Z}$ be ${\bf W}\bm{\Lambda}_{\bf Z}{\bf V}_{\bf
Z}^{\top}$. Minimizing over ${\bf Z}$ amounts to finding the matrices ${\bf
W}$, $\bm{\Lambda}_{\bf Z}$, and ${\bf V}_{\bf Z}$ that give the smallest
value of the objective. Because of the constraint ${\bf Z}{\bf
Z}^{\top}\preccurlyeq I_{d}$, $\bm{\Lambda}_{\bf Z}$ satisfies
$0\leq\lambda_{{\bf Z}\alpha}\leq 1$. Plugging this form into the optimization
objective we have

$\operatorname{Tr}\left[(c+{\epsilon}){\bf Z}^{\top}{\bf Z}-{\bf Z}^{\top}{\bf
Z}\big{[}{\epsilon}{\bf X}^{\top}({\bf X}{\bf X}^{\top})^{-1}{\bf X}+{\bf
Y}^{\top}{\bf
Y}\big{]}\right]=\operatorname{Tr}\left[(c+{\epsilon})\bm{\Lambda}_{\bf
Z}^{\top}\bm{\Lambda}_{\bf Z}-\bm{\Lambda}_{\bf Z}^{\top}\bm{\Lambda}_{\bf
Z}\tilde{\bf V}^{\top}\bm{\Lambda}^{\top}\bm{\Lambda}\tilde{\bf V}\right]$

where $\tilde{\bf V}={\bf V}^{\top}{\bf V}_{\bf Z}$. Note that ${\bf W}$ does
not appear in this expression and therefore remains undetermined in this
optimization objective.

${\bf V}$ only appears in the second term and hence can be found by optimizing
the second term alone. This term is of the form
$-\sum_{\alpha,\beta}\lambda_{{\bf
Z}\alpha}^{2}|\tilde{V}_{\alpha\beta}|^{2}\lambda_{\beta}^{2}$. The
$\tilde{\bf V}$ that minimizes the term is the one that pairs
$\lambda_{\beta}$ and $\lambda_{{\bf Z}\alpha}$, in the sorted order. Assuming
$(\lambda_{{\bf Z}\alpha})$ are also sorted in the decreasing order, this is
achieved by setting ${\bf V}_{\bf Z}={\bf V}$.

Plugging this in, the entire expression simplifies to

$-\sum_{\alpha}(\lambda_{\alpha}-c-{\epsilon})\lambda_{{\bf Z}\alpha}.$

Since $0\leq\lambda_{{\bf Z}\alpha}\leq 1$, it is clear that, for a minimum,
if $(\lambda_{\alpha}-c-{\epsilon})<0$, $\lambda_{{\bf Z}\alpha}=0$, while
$(\lambda_{\alpha}-c-{\epsilon})>0$ implies $\lambda_{{\bf Z}\alpha}=1$. In
case of equality, the choice of $\lambda_{{\bf Z}\alpha}$ does not matter.
Hence $\lambda_{{\bf Z}\alpha}=\mathbbm{1}(\lambda_{\alpha}>c+{\epsilon})$ is
a valid minimum. The vectors $\\{{\bf w}_{\alpha}\\}$ are the columns of ${\bf
W}$, while $\\{{\bf v}_{\alpha}\\}$ are the columns of ${\bf V}$. Hence, the
result. ∎

From this theorem, we derive two corollaries, describing the different limits
of ${\epsilon}\to 0$ and ${\epsilon}\to\infty$. The proofs of these
corollaries follow directly from Theorem 1 which maps the within-layer
computation of CCPC onto a PCA problem, and the known properties of PCA.

######  Corollary 2.

If the input sample covariance matrix, ${\bf C}_{x}$, is full rank, the
internal representation ${\bf Z}$ in the limit ${\epsilon}\rightarrow 0$ is
driven by the principal components (scores) of ${\bf Y}$, i.e.,

$\displaystyle{\bf Z}^{*}={\bf W}_{a}^{\top}{\bf Y}+{\epsilon}{\bf W}_{b}{\bf
X},$ (29)

with ${\bf W}_{a}$ the top principal directions of ${\bf Y}$ and ${\bf W}_{b}$
the maximally correlated directions between ${\bf X}$ and the columns of ${\bf
W}_{a}$.

######  Corollary 3.

If the input sample covariance matrix, ${\bf C}_{x}$, is full rank, the
internal representation ${\bf Z}$ in the limit ${\epsilon}\rightarrow\infty$
is driven by the maximally correlated directions between ${\bf
C}_{x}^{-1/2}{\bf X}$ and ${\bf Y}$, i.e.,

$\displaystyle{\bf Z}^{*}={\bf W}_{a}^{\top}{\bf Y}+{\bf W}_{b}{\bf X},$ (30)

with ${\bf W}_{a},{\bf W}_{b}$ projecting onto a common lower-dimensional
subspace so that the projections are maximally correlated.

###  B.3 Rank structure of CCPC

In this section again we set $c^{(l)}=0$ and $\rho^{(l)}=1$, and assume
$g_{a}=1,g_{b}=\epsilon$ is the same for all layers. The generalization to
other cases is straightforward. We demonstrate a series of rank constraints on
the intermediate layers of CCPC. The first of these is very general and
requires no assumptions. The following statements become gradually more tight
but also require more assumptions.

In this section, we define $r_{x}$ and $r_{y}$ to be the ranks of the the
covariances of ${\bf X}$ and ${\bf Y}$, and $r^{(l)}$ to be the rank of the
covariance of ${\bf Z}^{(l)}$ at the equilibrium of the CCPC update equations.
We take $r^{(0)}=r_{x}$ and $r^{(n)}=r_{y}$.

######  Corollary 4.

The rank of the $l^{\text{th}}$ layer is bounded by $r^{(l+1)}+r^{(l-1)}$:

$\forall l\;:\;r^{(l)}\leq r^{(l+1)}+r^{(l-1)}.$ (31)

###### Proof.

This follows directly from Theorem 1 which describes $Z^{(l)}$ as the whitened
PCA of the concatenation of ${\bf Z}^{(l-1)}$ and ${\bf Z}^{(l+1)}$. ∎

When the adjacent layers are whitened, i.e., the non-zero eigenvalues of
$Z^{(l+1)}Z^{(l+1)\top}$ and $Z^{(l-1)}Z^{(l-1)\top}$ are all one, we can
prove a tighter form of the rank constraint. This is relevant for the
intermediate layers of CCPC other than the final hidden layer. This
proposition is relevant for the final hidden layer only if the supervisory
input ${\bf Y}$ is white.

######  Proposition 3.

Let $r_{c}^{(l)}$ be the rank of the covariance matrix ${\bf Z}^{(l-1)}{\bf
Z}^{(l+1)\top}$. The rank of intermediate layers $l\leq n-2$ is given by:

$\forall l\leq
n-2\;:\;r^{(l)}\leq\begin{cases}r_{c}^{(l)}\leq\min\\{r^{(l+1)},r^{(l-1)}\\}&{\epsilon}\geq
1,\\\ r^{(l+1)}&{\epsilon}<1.\end{cases}$ (32)

###### Proof.

This is a slight generalization of [71]. It follows from the eigenvalue
structure of the concatenated data matrix given in Theorem 1, where only a
subset of the eigenvalues exceed the threshold given by ${\epsilon}$. The
number of eigenvalues that exceed this threshold is given by Eq. (32). ∎

Applying this proposition to all the layers of a CCPC network, we arrive at
the following corollary:

######  Corollary 5.

Assume that ${\bf Y}$ is white and ${\epsilon}\geq 1$. In a CCPC network with
$c^{(l)}\geq 0$, $\rho^{(l)}=1$, and layer widths $w^{(l)}$, the rank of the
latent variables ${\bf z}^{(l)}$ is less than or equal to
$\min\\{r_{x},r_{y},w^{(1)},\cdots,,w^{(n-1)}\\}$.

In other words, the rank of the latent variables is bounded by the minimum of
the input/output ranks and the bottleneck of the network.

###  B.4 Alignment and relation to CCA

In the limit where the internal loss terms in Eq. (6) are zero, we effectively
have a deep linear network [76]. A lot is known about deep linear networks,
trained by gradient descent of weights. One of the remarkable properties of
such networks is that the weights remain balanced, namely ${\bf
W}^{(l)\top}{\bf W}^{(l)}-{\bf W}^{(l-1)}{\bf W}^{(l-1)\top}$ remains a
constant, under gradient descent. In the case this difference is zero, the
left singular vectors of ${\bf W}^{(l-1)}$ match the right singular vectors of
${\bf W}^{(l)}$ ji2019gradient. In addition, for two class linear
discriminants, ${\bf W}^{(l-1)}$ ultimately become rank 1, under gradient
descent ji2019gradient.

##  Appendix C Numerical simulations

Figure 4: Dependence of BioCCPC performance on constraint scale. Unless
otherwise noted, this is the predictive-coding loss (Eq. 6) after 500,000
iterations. (a) A large network with two hidden layers of sizes 50 and 5
trained on MNIST. (b) The same large network as in (a) run on the Mediamill
dataset. See text. (c) A small net (single hidden layer of size 5) on
Mediamill. (d) Mean-squared prediction error for the large net on Mediamill,
showing that the best prediction performance may occur at a different choice
of $\rho$ compared to the lowest predictive-coding loss.

###  C.1 Effect of constraint scale

The BioCCPC network uses an inequality constraint on the covariance of the
hidden variables, $\frac{1}{T}{\bf Z}^{(l)}Z^{(l)\top}\preceq\rho^{(l)}{\bf
I}$, and this constraint is saturated at convergence.666More precisely, the
constraint is saturated only within the subspace that ${\bf Z}^{(l)}$
occupies. In general, this is only the case in the layer with the smallest
dimension, as the activity in the other layers is restricted to a subspace of
dimension given by the smallest layer; see Corollary 5. The scale,
$\rho^{(l)}$, is therefore important: a scale that is too small or too big
forces the hidden-layer activities into a range that can hinder learning, as
seen in Fig. 4a. Interestingly, the BioCCPC algorithm with a well-chosen
constraint scale converges faster than a plain PC network, that is, it can
reach a lower PC loss given the same number of training samples (Fig. 4a).

The predictive-coding loss, Eq. (6), penalizes prediction errors in all
layers. If we think about the task as supervised learning, however, it might
make more sense to focus only on the accuracy of the prediction in the last
layer. The best output prediction error may require a different constraint
scale compared to the lowest PC loss, cf. Fig. 4b and Fig. 4d.

###  C.2 Effect of dataset choice

Figure 5: Comparison of our algorithm with PC on the Mediamill dataset.
Simulations run using a network with two hidden layers, sizes 50 and 5. Shaded
area represents the variation over 10 runs. (a) Evolution of the predictive-
coding loss, Eq. (6), on a validation set during learning. (b) and (c)
Evolution of the eigenvalues of the covariance matrix of hidden-layer
activations in the two hidden layers. The dashed lines show the scale of the
constraint; see SM Sec. A. Inset: covariance matrix at end of training. Note
how activity is whitened in the second hidden layer (c), and low-rank in the
first hidden layer (b). Note also that the network performance is good even
though the constraint is not yet fully converged.

While we derived our algorithm by starting from the supervised PC framework of
[4], the same derivation would hold for the unsupervised or self-supervised PC
networks. In these settings, PC can be thought of as a framework which
combines information from the inputs at the first and last layer, rather than
a method for predicting an output from an input. In the linear case, this has
much in common with biologically plausible circuits that perform canonical
correlation analysis (CCA) [55, 71]. We therefore tested the algorithm on
Mediamill snoek2006challenge, a dataset consisting of two views of the same
scene, one given by the video and the other by text annotation. The results
(Figs. 4b-d and Fig. 5) show similar trends as those obtained on MNIST. In
particular, an appropriate choice of constraint scale leads to results that
are competitive with less biologically plausible implementations of PC.

In Fig. 6 we show the evolution of the validation-set loss for four other
datasets: FashionMNIST, CIFAR10, CIFAR100, and the LFW (Labeled Faces in the
Wild). The first three are supervised-learning scenarios and we use a network
similar to the one we employed for MNIST. The LFW dataset contains pairs of
views of the same subject, and so we employ a symmetric autoencoder-like
architecture for this self-supervised task. As before, we see that the BioCCPC
algorithm converges to a solution that has predictive-coding loss comparable
to the minimum obtained from traditional PC implementations (see Fig. 6).

Figure 6: Comparison of our algorithm with PC on a variety of datasets. The
name of the dataset and the architecture we employed are shown below each
figure. Shaded area represents the variation over 10 runs.

###  C.3 Network architecture

We tested BioCCPC across different architectures, such as a shallow network
with a small hidden-layer dimension (Fig. 4) and a deeper network with two
hidden layers, one of which is large-dimensional (Fig. 4a,b,d). In both cases
BioCCPC compares well to PC.

###  C.4 Lack of symmetry in feedforward and feedback weights

One of the difficulties with a biological implementation of predictive coding
is that it requires symmetric weights across layers. This is not so much an
issue of plausibility—indeed, symmetric weights can be obtained simply by
using the same Hebbian rules to train both feedforward and feedback
connections—but rather a problem with realism: weights in many areas of the
brain are not in fact seen to be symmetric.

Figure 7: Comparison of feedforward and feedback weights after learning on a
network with hidden layers of sizes 50 and 5 trained on the Mediamill dataset.
Figure 8: Comparison of feedforward and feedback weights after learning for
different network sizes. This is trained on MNIST. (a) There is less symmetry
in a smaller network with hidden layers of sizes 10 and 5. (b) There is more
symmetry in a larger network, with hidden layers 25 and 5.

Our algorithm explicitly breaks the symmetry between the weights, leading to
solutions that generically alleviate this difficulty. Indeed, in some cases
there is very little correlation between feedforward and feedback weights; see
Fig. 7. However, in some cases we see that the forward and backward weights
converge to similar values even though the learning rules are not symmetric.
We suspect that this can occur in cases where there is a large imbalance
between the width of adjacent layers, the larger the imbalance (i.e. the lower
rank the wider layer is) the more correlation we see in our experiments; see
Fig. 8. In biologically relevant neural networks, we do not generally see this
large imbalance in the width of adjacent layers and therefore we do not expect
to see high correlation between the forward and backward weights.

###  C.5 Predictive coding with constraint

In the main text we considered a predictive-coding network that imposes a
covariance constrained but we did not make explicit what this entails.
Explicitly, we augmented the Whittington-Bogacz objective, Eq. (2), with a
Lagrange multiplier enforcing an inequality constraint on the covariance of
hidden variables, similar to what we do in BioCCPC:

$\displaystyle L_{\text{cons}}$
$\displaystyle=\frac{1}{2}\sum_{l}\frac{\bigl{\|}{\bf Z}^{(l)}-{\bf
W}^{(l-1)}{\bf
Z}^{(l-1)}\bigr{\|}_{F}^{2}}{\sigma^{(l)2}}+\operatorname{Tr}{\bf
Q}^{(l)\top}{\bf Q}^{(l)}({\bf Z}^{(l)}{\bf Z}^{(l)\top}-T\rho^{(l)}{\bf I}),$
(33)

and optimized using stochastic gradient descent. This lets us distinguish the
effect of the constraint from the effect of the approximations we used in the
derivation.

###  C.6 Details of simulation runs

We initialized our the weights in our networks using uniform random numbers in
the range $[-a,a]$, where $a=1/\sqrt{\text{number of columns}}$. This is
similar to Xavier initialization glorot2010understanding except the scale of
initial weights depends only on the number of columns instead of both rows and
columns. We have found this initialization to provide slightly better results.

We optimized hyperparameters using Optuna optuna_2019 with the tree-structured
Parzen Estimator (TPE) algorithm bergstra2011algorithms. We ran the
optimization for 100 trials and repeated it 10 times with different random
seeds for every combination of dataset, network architecture, and constraint
scale $\rho$ (with the exception of the plain PC algorithm which does not have
a constraint). In each case we optimized the learning rate $\eta$ used in the
weight updates and the learning rate $\tau^{-1}$ of the fast dynamics (see
Algorithm 1). For networks that have a constraint, we allowed a different
learning rate $f\eta$ for the ${\bf Q}$ dynamics, and optimized for the factor
$f$. Finally, for BioCCPC network, we allowed for different learning rates for
the forward and backward weights, ${\bf W}_{b}$ and ${\bf W}_{a}$, and
optimized for the ratio between these learning rates. Hyperparameter
optimization runs were trained for 500 batches with batch size 100 and were
evaluated on a held-out validation set comprising 500 samples. Each run was
repeated 4 times with different random initialization to account for
stochasticity in the learning dynamics.

The simulation runs that were used to make the figures in the paper were based
on the hyperparameter optimization results that yielded the lowest predictive-
coding loss. Since hyperparameter optimization tends to push the learning rate
up to the brink of instability, the learning rate used in the simulations was
lowered by multiplying by a certain “safety” factor $\alpha$. We generally
used $\alpha=0.8$.

We ran our simulations on the CPU. Each run took from 1 to 4 minutes when run
single-threaded on a 24-core Intel Xeon CPU E5-2643 at 3.40GHz running Linux.
We also used a cluster to run the hyperparameter optimization for different
architectures in parallel.

##  Appendix D Generalization error

Figure 9: There is little overfitting in the linear setting that we are
investigating. The plots compare the training- and validation-set predictive-
coding losses for both the less biologically plausible predictive-coding
implementation from [4] (PC) and for our algorithm (BioCCPC).

The linear constraint that we impose on our network is a strong regularizer.
Figure 9 compares the training and validation losses during training on six
different datasets (as described above). We see that there is very little
overfitting in all cases.

##  Appendix E BioCCPC Pseudocode

Algorithm 1 :  Covariance Constrained predictive coding algorithm (BioCCPC)

input: ${\bf x}_{t}$, ${\bf y}_{t}$ $\triangleright$ new sample and previous
weight matrices

input: ${\bf W}_{a}^{(l)},{\bf W}_{b}^{(l-1)},{\bf Q}^{(l)}\qquad\forall
l\in\\{1,\ldots,n-1\\}$

${\bf z}_{t}^{(0)}\leftarrow{\bf x}_{t},\;{\bf z}_{t}^{(l)}\leftarrow{\bf
W}_{b}^{(l-1)}{\bf z}_{t}^{(l-1)}$ $\triangleright$ initialize latents via
forward pass

run until convergence for all $l\in\\{1,n-1\\}$ $\triangleright$ neural
dynamics

${\bf n}^{(l)}\leftarrow(1/g_{a}^{(l)}){\bf Q}^{(l)}{\bf z}_{t}^{(l)}$

${\bf v}_{a}^{(l)}\leftarrow{\bf W}_{a}^{(l)\top}{\bf z}_{t}^{(l+1)}-{\bf
Q}^{(l)\top}{\bf n}^{(l)}$ , ${\bf v}_{b}^{(l)}\leftarrow{\bf
W}_{b}^{(l-1)}{\bf z}_{t}^{(l-1)}$

${\bf z}_{t}^{(l)}\leftarrow{\bf
z}_{t}^{(l)}+\tau^{-1}\left[-g_{\text{lk}}^{(l)}{\bf
z}_{t}^{(l)}+g_{a}^{(l)}({\bf v}_{a}^{(l)}-{\bf z}_{t}^{(l)})+g_{b}^{(l)}({\bf
v}_{b}^{(l)}-{\bf z}_{t}^{(l)})\right]$

update for all $l\in\\{1,n-1\\}$ $\triangleright$ synaptic weight updates

${\bf W}_{a}^{(l)}\leftarrow{\bf W}_{a}^{(l)}+\eta\bigl{(}{\bf
z}_{t}^{(l+1)}{\bf z}_{t}^{(l)\top}-{\bf W}_{a}^{(l)}\bigr{)}$

${\bf W}_{b}^{(l-1)}\leftarrow{\bf
W}_{b}^{(l-1)}+\eta\Big{[}g_{a}^{(l)}\bigl{(}{\bf W}_{a}^{(l)\top}{\bf
z}_{t}^{(l+1)}-{\bf Q}^{(l)\top}{\bf n}^{(l)}\bigr{)}-c^{(l)}{{\bf
z}_{t}^{(l)}}\Big{]}{{\bf z}_{t}^{(l-1)\top}}$

${\bf Q}^{(l)}\leftarrow{\bf Q}^{(l)}+\eta\bigl{(}{\bf n}^{(l)}{\bf
z}_{t}^{(l)\top}-{\bf Q}^{(l)}\bigr{)}$

output: ${\bf z}_{t}^{(l)}\qquad\forall l\in\\{0,\ldots,n\\}$ $\triangleright$
internal representations and updated weights

output: ${\bf W}_{a}^{(l)},{\bf W}_{b}^{(l-1)},{\bf Q}^{(l)}\qquad\forall
l\in\\{1,\ldots,n-1\\}$

unsrt biblio