# Tightening the Biological Constraints
on Gradient-Based Predictive Coding

Nick Alonso<EMAIL_ADDRESS>1234-5678-9012 and  Emre Neftci<EMAIL_ADDRESS>University of California, IrvineP.O. Box 1212IrvineCAUSA92617

(2021)

###### Abstract.

Predictive coding (PC) is a general theory of cortical function. The local,
gradient-based learning rules found in one kind of PC model have recently been
shown to closely approximate backpropagation (Whittington and Bogacz, 2017;
Millidge et al., 2020; Song et al., 2020). This finding suggests that this
gradient-based PC model may be useful for understanding how the brain solves
the credit assignment problem. The model may also be useful for developing
local learning algorithms that are compatible with neuromorphic hardware. In
this paper, we modify this PC model so that it better fits biological
constraints, including the constraints that neurons can only have positive
firing rates and the constraint that synapses only flow in one direction. We
also compute the gradient-based weight and activity updates given the modified
activity values. We show that, under certain conditions, these modified PC
networks perform as well or nearly as well on MNIST data as the unmodified PC
model and networks trained with backpropagation.

Predictive Coding, Gradient-based Learning, Local Learning, Synaptic
Plasticity

††copyright: acmcopyright††journalyear: 2021††doi:
10.1145/1122445.1122456††conference: Oakridge ’21: International Conference on
Neuromorphic Systems; July 27–29, 2021; Oakridge, IL††booktitle: Oakridge, IL
’21: International Conference on Neuromorphic Systems, July 27–29, 2021,
Oakridge, IL††price: 15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Computing
methodologies Neural networks

##  1\. Introduction

Predictive coding (PC) is a general theory of the function of top-down and
bottom-up processing in the neocortex (Rao and Ballard, 1999a; Spratling,
2008; Friston and Kiebel, 2009). According to the theory, a central function
of top-down projections in the cortex, which are connections that lead from
higher to lower level areas of the cortical hierarchy, is to predict neural
activity. Differences between the predictions and the actual activity is
computed and encoded in error neurons, which feed the errors back up the
cortical hierarchy from lower to higher level areas. Learning dynamics and
neural activity updates both seek to minimize the prediction errors. While
empirical testing of PC theory is needed to validate it, the theory has seen
significant empirical support and is consistent with much of what we know
about neuroanatomy (Huang and Rao, 2011; Kok and de Lange, 2015; Walsh et al.,
2020), making it a potential basis for a general theory of how the neocortex
learns and performs inference.

Recent computational work on PC theory has shown that, under certain
conditions, a seminal model of PC, originally developed by Rao and Ballard
(Rao and Ballard, 1999b), learns in a way that is very similar to, or in some
cases exactly similar to, the gradient back-propagation of errors algorithm
(Whittington and Bogacz, 2017; Millidge et al., 2020; Song et al., 2020). This
finding is significant from the point of view of neuroscience because PC
networks learn using local, Hebbian-like learning rules, and therefore PC
theory could provide a basis for an explanation of how the cortex solves the
credit assignment problem. The theory may also advance neuromorphic computing
by providing a path toward local learning algorithms that are compatible with
neuromorphic hardware (Davies et al., 2018; Qiao et al., 2015; Friedmann et
al., 2016).

In this work, we modify Rao and Ballard’s (Rao and Ballard, 1999b) PC model
with the goal of making the model more tightly constrained by biology. Doing
so will yield both a model that better fits with neurobiology and will be a
step toward implementing gradient-based PC algorithms in spiking neural
networks. We make and test several modifications to the model. First, the
weights that propagate errors in the original PC model are transposes of the
forward weights. To avoid this biologically implausible transport of weights,
we replace the weight transpose with a separate weight matrix. We test both
using the weight matrix untrained (random feedback) and trained using a rule
from (Kolen and Pollack, 1994). Second, the neurons in Rao and Ballard model
encode continuous valued firing rates and the values of activity neurons
(green nodes in figure 1) are allowed to be negative. Firing rates can only be
positive in artificial spiking neurons and biological neurons, so we prevent
values from going negative using a ReLU activation function and develop a
method for preventing the loss of gradient information due to the ReLU
function. Third, it is necessary to report negative errors, but neural
activities cannot be negative. We discuss several ways errors can be encoded
in spiking neurons with only (positive) firing rate values. We also compute
the gradients for these error encoding schemes for firing rate neurons and use
them to update the model.

Figure 1. Standard Predictive Coding Model

We test these modifications on a supervised learning task using the MNIST
dataset. We find that certain versions of the modified PC model works as well
as the original PC model and backpropagation. This suggests the PC model may
be a favorable basis on which to develop more neuromorphic and biologically-
constrained spiking neural network models that use gradient-based learning.

##  2\. Gradient-Based Predictive Coding

Term | Value
---|---
$a_{l}$ | activity at level $l$
$e_{l}$ | prediction error at level $l$
$p_{l}$ | prediction at level $l$
$W_{l}$ | synaptic weights from $l$ to $l+1$

In this section, we describe the standard gradient-based PC model based on Rao
and Ballard’s model (Rao and Ballard, 1999b). We choose to focus on this model
because, under certain conditions, the model was demonstrated to closely
approximate gradient backpropagation (Whittington and Bogacz, 2017; Millidge
et al., 2020; Song et al., 2020). Song et al. (Song et al., 2020), for
example, showed that the first non-zero weight update that occurs to each
weight matrix during supervised training is exactly identical to those made by
backpropagation. Whittington et al. (Whittington and Bogacz, 2017) showed that
after the network converges, the equations that define how errors are
propagated through the network match those of the backpropagation algorithm.
This relation is important because it provides a formal link between PC and
the minimization of a global objective function through gradient ascent, which
is a method that works well in artificial neural networks (ANN) at scale.
Other predictive coding that lack this connection have no guarantees on their
ability to minimize global objectives.

There are two kinds of neurons in the PC model: activity neurons and
prediction error neurons. Activity neurons represent features of the input,
while error neurons encode differences between predictions of activity neurons
and the actual values of activity neurons. Predictions of activities at level
$l$ are passed from level $l-1$ according to the following equation:

(1) $p_{l}=W_{l-1}a_{l-1}.$

A local prediction error is then computed and encoded in error neurons. The
error is the difference between the activity and the prediction, which is
first passed through an element-wise non-linearity $f$.

(2) $e_{l}=a_{l}-f(p_{l}).$

The weight matrices $W_{l}$ are trained using this local error according to:

(3) $\begin{split}W_{l}\leftarrow\alpha\Delta W_{l}+W_{l},\text{ where }\Delta
W_{l}&=e_{l+1}f^{\prime}(p_{l+1})a_{l}^{T}.\end{split}$

$\Delta W_{l}$ is multiplied by a small learning rate $\alpha$ before being
added to the weights, where $0<\alpha<1$. This update is equivalent to taking
the gradient of the squared prediction error $\frac{1}{2}e_{l}^{2}$ with
respect to the weights that produced the prediction (Rao and Ballard, 1999b).

In addition to updating their weights, PC networks also optimize their
activities. Activities optimize the same cost, i.e. the squared prediction
errors at each level. At each time-step, they update their activities
according to:

(4) $\Delta a_{l}=W_{l}^{T}e_{l+1}f^{\prime}(p_{l+1})-e_{l}.$

Here, $f^{\prime}$ is the gradient of the non-linearity. Updating activities
at level $l$ is equivalent to performing gradient ascent over the activities
using the gradient of the squared prediction error at level $l$ and the next
level $l+1$. Because one is performing gradient ascent, the computed gradient
must be multiplied by a small inference rate $\beta$, where $0<\beta<1$,
before adding to the activities. One typically wants the activities to update
quicker than the weights, so $\beta$ should be larger than $\alpha$. This is
because one wants activities to converge quickly to an optimal representation
of the input, while learning (i.e. weight change) should occur slowly to
ensure generalization ability.

##  3\. Predictive Coding Constrained

We now describe our modified PC network. Each modification to the network is
made to make the model fit better with neurobiology. However, we wish to make
these changes without losing the ability of the network to propagate gradients
and perform similarly to backpropagation.

###  3.1. Predictive Coding without Weight Transport

Equation 4 shows that errors propagated back through the network require using
the transpose of the forward weights. Using the transpose of forward weights
to propagate errors is generally considered biologically implausible as it
implies a bidirectional synapse. Transposing a weight matrix is also
problematic from a neuromorphic hardware implementation perspective because it
requires the duplication and synchronization of the weights on the pre-
synaptic and the post-synaptic side.

The weight transport problem has been tackled in the context of gradient-
backpropagation for conventional deep networks using approximations such as
random feedback weights (Lillicrap et al., 2016) and local loss functions
(Mostafa et al., 2018; Akrout et al., 2019; Kunin et al., 2020). These
approximations underperform slightly compared to exact gradient back-
propagation, but do not require a symmetric transpose of the network weights.

We propose two modifications to PC networks that avoid this problem, which
build off of previous work done on the weight transport problem for
backpropagation.

In the first network model, which we call Rand-PC, we replace the $W^{T}$ in
equation 4 with a random matrix $B$ of the same size. Replacing weight
transposes with random feedback matrices has been shown to work reasonably
well for approximating backpropagation in conventional deep networks
(Lillicrap et al., 2016; Xiao et al., 2018). In the results section, we show
it works nearly as well as backpropagation in our predictive coding networks.

In the second network model, which we call the Kollen-Pollack PC (KP-PC), we
use a method developed by Kolen and Pollack (Kolen and Pollack, 1994), and
extended by Akrout et al. (Akrout et al., 2019), to train backward matrices so
that they converge to the transposes of the forward matrices. Like Rand-PC,
with KP-PC we replace the weight transpose in equation 4 with a separate
weight matrix $B_{l}$. This weight matrix is then trained with the transpose
of the update to the forward $W_{l}$. Kollen and Pollack build on the simple
idea that, if one makes the same weight update to two matrices, $W$ and $B$,
the two matrices will grow more similar. However, due to numerical precision
errors, the two matrices in practice eventually diverge. To solve this problem
in conventional ANN, Kolen and Pollack add a small decay term to both updates.
They show that $W$ and $B$ will eventually converge to (nearly) the same
matrix.

(5) $\begin{split}\Delta W_{l}&=A_{l}-\gamma W_{l},\,\Delta
B_{l}=A_{l}^{T}-\gamma B_{l},\\\ \end{split}$

where $0<\gamma<1$ is a small decay rate and $A_{l}$ is the weight adjustment.
Acrout et al. (Akrout et al., 2019) point out that these weight updates are
local for the forward and backward matrices, if each weight matrix is
connected to the same two populations of neurons. Interestingly, as can be
seen from the equations in section 2 and figure 1, PC networks require just
this structure, where forward weights take input from activities at level $l$
and output to error neurons at level $l+1$, while backward weights to the
opposite.

###  3.2. Biologically Constrained Activity Neurons

Equation 4 shows that PC networks allow for negative neural activity values
because updates are linear. In biological neurons and artificial spiking
neurons, an activation $a$ typically represents the firing rate of a neuron,
which cannot be negative. Thus, in a neural modelling and neuromorphic
computing setting it is desirable to have neurons with non-negative
activities. This raises the question of how well PC models work under the
constraint that activities are non-negative. To test this, we enforce positive
activation values by passing the activities through a rectified linear unit
(ReLU) function after each activity update to ensure there are no negative
values.

Additionally, to help prevent a loss of gradient information due to the ReLU
function, we add a bias term to the model. In equation 2, for example, the
error $e_{l}$ is now equal to $a_{l}-f(p_{l})+b$. For the subtraction and
division threshold schemes, one simply replaces $f(p_{l})$ with $f(p_{l})+b$
in the equations where $f(p_{l})$ is present. No other changes are needed.
This simple addition to the model will help prevent activities from going
negative and will therefore help prevent gradient information from being lost
when negative activities are zeroed out. The bias will do this by forcing
predictions to be greater than zero. This will make sure the network
activities are initialized to non-negative values. It will also help
activities maintain non-negative during activity updates, through the top-down
error (second term in equation 2). The top-down error pulls the activities
toward the predictions at the same level. If predictions are always non-zero,
they will always pull activities toward non-zero values. If activities a non-
negative, the gradients they accumulate will not be lost when passed through a
ReLU function.

###  3.3. Biologically Constrained Error Neurons

Like the activity neurons, error neurons, according to equation 2, can take on
negative values when activity values are over-predicted. For the same reasons
listed in the previous section, it is desirable to have error neurons with
only positive firing rate values. This raises the question of how prediction
errors could be encoded in spiking patterns with non-negative firing rates.

Rao and Ballard and several other scientists (e.g.(Keller and Mrsic-Flogel,
2018)) hypothesize there could be two kinds of prediction error neurons in the
brain. One kind spikes in response to over-predicted values (i.e. negative
errors according to equation 2). The other kind spikes in response to under-
predicted values (i.e. positive errors according to equation 2). This encoding
scheme implies that generally, when prediction errors are large the error
neuron activities will be large, and while errors are small error neuron
activity is small. These are the sort of error neurons that are assumed to
exist in the standard PC model. Implementing this encoding scheme could, for
example, use one set of error neurons that encode only the positive values of
$a_{l}-f(p_{l})$ (underpredictions) and another set of error neurons encoding
only the positive values of $f(p_{l})-a_{l}$ (overpredictions). We will call
this sort of error encoding scheme subtractive separated encoding since the
subtractive error is encoded in two separate kinds of error neurons.

An alternative way to encode errors in spike trains involves setting error
neurons to have a baseline activity rate, which is dampened when activities
are over-predicted and excited when they are under-predicted. When subtraction
is used to dampen and excite error neurons, we call this encoding scheme the
subtractive threshold encoding scheme. There may be several ways to implement
this in spiking neurons. For example, one could use a constant input of
current that causes the error neurons to fire at a certain rate even when
errors are absent.

Here we show how a subtractive threshold scheme can be developed in firing
rate neurons so weight and activity updates are equivalent to the original Rao
and Ballard equations. Consider the following error neuron encoding scheme

(6)
$e_{l}^{*}=\frac{2}{e_{\max}}((a_{l}-f(p_{l})-e_{\min})=\frac{2(e_{l}-e_{\min})}{e_{\max}}.$

Here, $e_{\min}$ is the minimum possible value of $e_{l}$, and $e_{\max}$ is
the maximum possible value of $e_{l}-e_{\min}$. There will be a minimum value
if activities are clamped to have a minimum value and $f$ is a squashing non-
linearity (e.g. sigmoid). If there is no maximum one can replace $e_{\max}$
with an approximate $e_{\max}$ value to obtain a similar result. For example,
we apply the sigmoid function to $p_{l}$, and the ReLU function to the
activities. If activities never go above 1.1, then $e_{\min}=-1$ and
$e_{\max}=2.1$. This equation forces all error values to be between 0 and 2,
with a baseline firing rate of 1. One can replace the 2 with another value
which determines the maximum firing rate.

In spiking neurons, the $e_{\min}$ can be thought of as a constant excitatory
current source that, absent any other inputs, causes the neurons to
spontaneously fire at a constant rate. The $\frac{2}{e_{\max}}$ term can also
be seen as a constant input current that performs some form of normalization,
which is well known to be a pervasive computation throughout the brain
(Carandini and Heeger, 2012).

Weight and activity updates are computed by replacing the $e_{l}$ term in
equation 2 and 4 with $(\frac{e_{\max}}{2}e_{l}^{*}+e_{\min})$, since
$e_{l}=(\frac{e_{\max}}{2}e_{l}^{*}+e_{\min})$. With these replacements,
updates remain local (since they only depend on the local errors, local
activities, and constants) and are equivalent to the original equations.

We note that these two error encoding schemes with biological or spiking
neurons can be obtained as special cases of population decoding (Dayan and
Abbott, 2001), which can be conveniently realized with e.g. the neural
engineering framework (Eliasmith and Anderson, 2003). Our focus on these
special cases is due to their efficiency, as they require one or two neurons
per encoded error value compared to population of neurons per error value in
the general case.

A third way to encode mismatches between predictions and activities is an
encoding scheme that involves dividing the activities by the predictions (or
vice versa). Spratling (Spratling, 2008), for example, developed a firing rate
model of predictive coding which uses a division error term that divides
activities (element-wise) by the predictions. His model was able to replicate
some neurophysiological data, including fine-grained calcium imaging data that
seemed to show the existence of neurons in the mouse primary visual cortex
that were sensitive to mismatches between actual and predicted visual flow
(Spratling, 2019). We call encoding schemes of this form division mismatch
encoding.

Spratling’s particular model is not formulated to minimize a global loss using
its gradients w.r.t. weights and activities. We are interested in working
within the framework of gradient ascent as it has proven effective in large
neural networks. We thus present an alternative firing rate model with
division mismatch encoding that maintains gradient updates on activities and
errors. We compute values of the mismatch neurons as follows:

(7) $e_{l}^{**}=\sqrt{\frac{a_{l}+\epsilon}{f(p_{l})+\epsilon}}.$

The $\epsilon$ is a small constant that prevents division by zero and helps
prevent exploding gradients. The term
$\frac{a_{l}+\epsilon}{f(p_{l})+\epsilon}$ ranges between 0 and $\infty$
(assuming activities and predictions are positive). Over-predictions range
between $0$ and $1$, while under-predictions range between $1$ and $\infty$.
We find that, although not necessary, adding the square root function improves
learning slightly and makes activity updates more stable. When activities
equal the predictions, $e_{l}^{**}$ will equal one, so we develop a cost
function that measures the difference between $e_{l}^{**}$ and $1$ (see
appendix). Weights and activities are then updated in proportion to the
gradient of this new cost function (see appendix).

Data: ($X_{B},Y_{B}$) = dataset X and targets Y formed into minibatches B

begin

// For each minibatch perform n activity updates then a weight update

for  _$(x_{b},y_{b})\in(X_{B},Y_{B})$_ do

// Initialize

$a_{0}\leftarrow x_{b}$

for  _$l=1$ to $L$_ do

$p_{l}\leftarrow W_{l-1}a_{l-1}$

$a_{l}\leftarrow f(p_{l})$

end for

$a_{L}\leftarrow y_{b}$

// Perform $n$ Activity Updates

for  _$i=0$ to $n$_ do

for  _$l=1$ to $L$_ do

// Compute Errors

equation 2, 6, or 7

// Activity Update

equation 4

// Prediction Update

equation 2

end for

end for

// Train Weight Matrices

for  _$l=0$ to $L-1$_ do

// Compute Errors then update weights

equation 2, 6, or 7

equation 3

end for

end for

end

Algorithm 1 Supervised Predictive Coding Training

##  4\. Results

Predictive coding networks can be trained for both supervised and self-
supervised learning tasks. In self-supervised learning tasks, the activities
of the bottom level activities (level 0 in figure 1) are set to the data (e.g.
image), while the activities of the top level (level 4 in figure 1) are set to
some constant (e.g. 1). Alternatively, the top level can be removed entirely.
The network weights can then be trained online after each update to the
activities, or a single weight update can occur after the network activities
converge (Bogacz, 2017).

In what follows, we show the results of training our models using supervised
methods, which were developed in (Whittington and Bogacz, 2017). In supervised
learning, the bottom level activities of the PC network are clamped to the
target, while the top-level activities are clamped to the input data.
Activities are initialized by propagating predictions down level by level to
level 0. Activities of hidden layers are set to the predictions at each level:
$a_{l}^{initial}=f(p_{l}^{initial})$. After initialization, activities are
optimized. Weight updates can occur online or after the activities reach
convergence (See Algorithm 1). We found activity convergence is slow for
MNIST, so instead of waiting for convergence we update weights after $20$
activity updates. For fashion-MNIST, the activities converge more quickly, so
we update weights after $7$ activity updates. When tested, the the activities
at level $4$ are clamped to the test image, then predictions propagate down to
the output layer (level $0$), where the output is compared to the target.

Inference rates are set to .1 for MNIST and .025 for Fashion-MNIST. Learning
rates are set to .001. Adam optimizers are used to train all weight matrices.
All models use a network with fully connected layers of size of
784-300-300-10. We train the networks to classify images in the MNIST and
Fashion-MNIST data set. Each dataset consists of gray-scale images of size
28x28. There are 60,000 training images and 10,000 test images in each
dataset. We show the classification accuracies on the test sets below. 111Code
for models can be found here: https://github.com/nalonso2/Tightening-the-
Biological-Constraints-on-Gradient-Based-Predictive-Coding

###  4.1. PC with Separate Feedback Weights

Figure 2. Test errors of PC models with various feedback weights on MNIST with
sigmoid activations.

We begin by testing how well PC networks train with the separate feedback
matrices discussed in section 2.3. Here all networks use sigmoid activation
functions. Rand-PC uses a fixed, random feedback matrix, while the KP-PC uses
a weight matrix trained using the rules discussed in section 2.3. The test
accuracy, averaged over three runs, is shown after each epoch of training
starting after the first training epoch.

Data | Backprop | PC | KP-PC | Rand-PC | PC w/ Div | Rand-PC w/ Div
---|---|---|---|---|---|---
MNIST |  $.019$ $(1x10^{-3})$ |  $.020$ $(1x10^{-3})$ |  $.019$ $(2x10^{-3})$ |  $.021$ $(1x10^{-3})$ |  $.021$ $(8x10^{-4})$ |  $.022$ $(1x10^{-3})$
Fashion-MNIST |  $.107$ $(2x10^{-3})$ |  $.109$ $(3x10^{-3})$ |  $.116$ $(5x10^{-3})$ |  $.112$ $(3x10^{-3})$ |  $.116$ $(4x10^{-3})$ |  $.116$ $(2x10^{-3})$
Table 1. Mean (and standard deviation) of validation errors (%) for PC models
with sigmoid activitations. Three different seeds of each model were trained.
Values shown are the means and standard deviations of the validation errors of
the last three epochs across training runs.

We find the standard PC network, which uses the weight transpose for feedback,
performs as well as backpropagation, which replicates the findings of
(Whittington and Bogacz, 2017). The KP-PC network also performs as well as
backpropagation, while the Rand-PC network only does slightly worse (see table
1). All models achieve a mean accuracy within a standard deviation of $98\%$.
Similar results are found with fashion-MNIST, where these models’ mean test
errors were within two standard deviations of backpropagation.

###  4.2. PC with Constrained Activity Rates

Next, we test how well PC networks work under the constraint that activity
neurons can only take positive firing rate values. As mentioned above, because
firing rates are updated linearly, activity values can become negative even
when initialized to be positive. It is possible then, that preventing
activities from being negative (by passing them through a ReLU function after
each update) may erase gradient information useful for credit assignment, and
this may consequently hurt performance. Here we test how severe this potential
loss of gradient information is when using different activation functions. All
models have sigmoid activations at the output layer, but some models use
sigmoid functions at hidden layers while others use tanh functions at hidden
layers.

Figure 3. Test errors of PC models on MNIST with activities constrained to be
positive (+Act) versus those allowed to go negative.

Figure 3 shows that constraining activities to be positive does not negatively
affect the performance of the network when sigmoid activation functions are
used. However, there is a small but significant drop in performance when Tanh
non-linearities are used. When a small bias is added, as specified in section
3.3, the network using Tanh non-linearities sees no drop in performance.

There are a couple reasons why we see this pattern. When the network is
initialized, the activities are set to the predictions which are passed
through the non-linearity (see 1). The same process is used to generate
predictions at test time. Because sigmoid maps all values to non-negative
numbers, the network, when tested, will have activities with only positive
values. Thus, when ReLUs are applied, no activities are affected. This is not
true when Tanh is applied without a bias, which will map some numbers to
negative values. However, when a bias is added, the Tanh network will not
produce negative predictions.

However, this does not explain why learning in the sigmoid network without
bias is not affected by constrained activity values, since updates to the
activities that occur after the initialization can still push the activities
to negative values. If an update pushed activity values to be negative, then
the negative values will be set to zero and the gradient information they
encode will be lost. This loss of gradient information will presumably
negatively affect learning.

We find through experimentation, however, that the top-down error term (the
second term in equation 2) prevents activities from going too far outside the
range of the predictions. When predictions are passed through a sigmoid, for
example, the top-down error will pull the activities toward the prediction,
which is some positive number between zero and one. This generally prevents
activities from going below zero. The same is not true of predictions passed
through Tanh without a bias, which will not always pull activities toward
positive values. In this case, useful gradient information will be lost when
negative activities are clamped to zero. A bias term, thus, is necessary in
cases where non-linearities allow for negative predictions, but may not be
necessary for non-linearities that map all input to non-negative numbers.

###  4.3. Subtraction versus Division Errors

In section $4.3$, we outlined three error neuron encoding schemes: subtractive
separated encoding, subtractive threshold encoding, and division mismatch
encoding. We also showed how to compute the gradients for the new encoding
schemes. As explained above, the two subtractive error encoding schemes will
look quite different in spiking neurons. However, we showed that both of these
encoding schemes are mathematically equivalent in the firing rate model, and
that the separated subtractive encoding scheme works as well as
backpropagation. Here we test how well the division threshold encoding scheme
compares to backpropagation.

Positive activity values are required for the division encoding scheme (due to
the square-root and logarithm in equations 7 and 8), so we set all models to
have positive activity values. We apply sigmoid activation at hidden layers.
We still find that adding a small bias is necessary to prevent loss of
gradient information, so we add a small bias.

Figure 4. Test errors of PC models on MNIST with division encoding errors.

We can see in figure 4 that the PC models that use division based errors
produce comparable results to backpropagation. The division mismatch model
that uses true gradients achieves a mean test error of $.021$ while the
division model with random feedback achieves a mean test error of $.22$ (see
table 1), which is only slightly worse than the performance of
backpropagation, which had a mean test error of $.019$. The results are
similar for fashion MNIST as well, where both division mismatch models produce
mean test errors with less the one percent difference of the mean test error
of backpropagation (see 1).

##  5\. Discussion

In this paper, we showed that a more biologically constrained version of Rao
and Ballard’s (Rao and Ballard, 1999b) seminal model of predictive coding
performed similarly to backpropagation on supervised learning tasks using
MNIST data. We found this to be true under constraints where 1) separate
feedback weights were used to propagate errors, 2) activity values were
prevented from going negative, and 3) error neuron activities were prevented
from going negative using either division or subtraction based encoding
schemes. We also showed how the gradients for the new encoding schemes could
be computed and incorporated into the model.

These results suggest that it is likely possible for more biologically
constrained models of gradient-based predictive coding to be built using
spiking neural networks. We computed the gradients for the division and
subtraction threshold error encoding schemes, which prevent error neurons from
having negative activity rates. These equations can potentially be used as a
basis for forming new equations for spiking neuron models. We also discussed
how separated subtraction error encoding could be implemented in spiking
neurons. Although spiking models of PC have been previously demonstrated (e.g.
(Wacongne et al., 2012)), as far as we can find, spiking neural models of PC
that utilize gradient-based inference and learning have yet to be developed.

Of course, in spiking neurons the true gradients cannot be computed because
the step functions used as non-linearities are non-differentiable. However,
surrogate gradient methods (Neftci et al., 2019) used to approximate gradients
in spiking neural networks can naturally be incorporated into predictive
coding networks. Such surrogate-gradient spiking PC models could help further
develop the empirical hypothesis that predictive coding is the general
algorithm the brain uses to solve the credit assignment problem. Additionally,
such networks could lead to useful local-learning algorithms that are
compatible with neuromorphic hardware.

It is still an open question whether the cortex is performing some form of
predictive coding. There is good evidence that top-down connections in the
cortex do propagate predictions (Huang and Rao, 2011; Kok and de Lange, 2015;
Walsh et al., 2020). However, the hypothesis that the cortex encodes
prediction errors widely in specialized error neurons (which is a key
implication of PC) is not yet widely accepted within neuroscience. There is
much data that is consistent with the hypothesis that such error neurons exist
in the cortex, but none of it is particularly conclusive (for recent review
see (Walsh et al., 2020)).

One reason why it has been difficult to locate error neurons is that different
error encoding schemes yield different empirical predictions (Walsh et al.,
2020). However, some progress is being made. Spratling (Spratling, 2019), for
example, recently found that division-based error encoding schemes better fit
certain neurophysiological data than the subtractive encoding of the Rao-
Ballard model. Spratling did not compare the physiological data to a
subtractive threshold encoding scheme, like the one we proposed here (equation
6), so it is unclear whether division error encoding fits with the data better
than subtractive schemes generally or only with a particular type of
subtractive encoding present in the Rao-Ballard model. Nonetheless, this study
shows how we can begin to build evidence in favor of one hypothesis over
another.

Our model also illustrates that gradient-based PC is compatible with multiple
different kinds of error encoding schemes. So the general hypothesis that some
form of gradient-based PC is utilized by the cortex does not depend on there
existing a subtractive rather than divisive error encoding. It will, instead,
depend on the way the errors are propagated through the cortex and used to
affect neural activity and learning. It may be that errors and weights are
updated using some other optimization method (e.g. (Spratling et al., 2009)).

More empirical and computational work will be needed to settle these debates.
We hope, however, that the work presented here provides new avenues to further
develop and test neural models of gradient-based predictive coding and to
further develop useful brain-inspired learning algorithms.

##  6\. Conclusion

PC is a theory that is of interest to both cognitive science, neuroscience,
and engineering. It, for one, is able to explain a wide array of
neurophysiological and anatomical data. PC also has potential to provide a
path toward understanding how the brain solves the credit assignment problem.
We also suggested that the local learning rules used within PC models could
lead to useful learning algorithms that are compatible with the constraints
neuromorphic hardware. In this paper, we showed that the standard gradient-
based PC model can learn and perform inference well under tighter biological
constraints. This further supports the position that PC can be developed in a
way that is useful for both neuroscience and neuromorphic computing, and marks
a path along which PC can be further developed.

###### Acknowledgements.

This work was supported by the National Science Foundation under grant 1652159
and 1823366 (EN).

## References

  * (1)
  * Akrout et al. (2019) Mohamed Akrout, Collin Wilson, Peter C Humphreys, Timothy Lillicrap, and Douglas Tweed. 2019\.  Deep learning without weight transport.  _arXiv preprint arXiv:1904.05391_ (2019).
  * Bogacz (2017) Rafal Bogacz. 2017\.  A tutorial on the free-energy framework for modelling perception and learning.  _Journal of mathematical psychology_ 76 (2017), 198–211.
  * Carandini and Heeger (2012) Matteo Carandini and David J Heeger. 2012.  Normalization as a canonical neural computation.  _Nature Reviews Neuroscience_ 13, 1 (2012), 51–62.
  * Davies et al. (2018) Mike Davies, Narayan Srinivasa, Tsung-Han Lin, Gautham Chinya, Yongqiang Cao, Sri Harsha Choday, Georgios Dimou, Prasad Joshi, Nabil Imam, Shweta Jain, et al. 2018\.  Loihi: A neuromorphic manycore processor with on-chip learning.  _Ieee Micro_ 38, 1 (2018), 82–99.
  * Dayan and Abbott (2001) Peter Dayan and Laurence F Abbott. 2001.  _Theoretical neuroscience: computational and mathematical modeling of neural systems_.  Computational Neuroscience Series.
  * Eliasmith and Anderson (2003) Chris Eliasmith and Charles H Anderson. 2003.  _Neural engineering: Computation, representation, and dynamics in neurobiological systems_.  MIT press.
  * Friedmann et al. (2016) Simon Friedmann, Johannes Schemmel, Andreas Grübl, Andreas Hartel, Matthias Hock, and Karlheinz Meier. 2016.  Demonstrating hybrid learning in a flexible neuromorphic hardware system.  _IEEE transactions on biomedical circuits and systems_ 11, 1 (2016), 128–142.
  * Friston and Kiebel (2009) Karl Friston and Stefan Kiebel. 2009.  Predictive coding under the free-energy principle.  _Philosophical Transactions of the Royal Society B: Biological Sciences_ 364, 1521 (2009), 1211–1221.
  * Huang and Rao (2011) Yanping Huang and Rajesh PN Rao. 2011.  Predictive coding.  _Wiley Interdisciplinary Reviews: Cognitive Science_ 2, 5 (2011), 580–593.
  * Keller and Mrsic-Flogel (2018) Georg B Keller and Thomas D Mrsic-Flogel. 2018.  Predictive processing: a canonical cortical computation.  _Neuron_ 100, 2 (2018), 424–435.
  * Kok and de Lange (2015) Peter Kok and Floris P de Lange. 2015.  Predictive coding in sensory cortex.  In _An introduction to model-based cognitive neuroscience_. Springer, 221–244.
  * Kolen and Pollack (1994) John F Kolen and Jordan B Pollack. 1994.  Backpropagation without weight transport. In _Proceedings of 1994 IEEE International Conference on Neural Networks (ICNN’94)_ , Vol. 3. IEEE, 1375–1380.
  * Kunin et al. (2020) Daniel Kunin, Aran Nayebi, Javier Sagastuy-Brena, Surya Ganguli, Jonathan Bloom, and Daniel Yamins. 2020.  Two routes to scalable credit assignment without weight symmetry. In _International Conference on Machine Learning_. PMLR, 5511–5521.
  * Lillicrap et al. (2016) Timothy P Lillicrap, Daniel Cownden, Douglas B Tweed, and Colin J Akerman. 2016.  Random synaptic feedback weights support error backpropagation for deep learning.  _Nature communications_ 7, 1 (2016), 1–10.
  * Millidge et al. (2020) Beren Millidge, Alexander Tschantz, and Christopher L Buckley. 2020\.  Predictive coding approximates backprop along arbitrary computation graphs.  _arXiv preprint arXiv:2006.04182_ (2020).
  * Mostafa et al. (2018) Hesham Mostafa, Vishwajith Ramesh, and Gert Cauwenberghs. 2018\.  Deep supervised learning using local errors.  _Frontiers in neuroscience_ 12 (2018), 608.
  * Neftci et al. (2019) Emre O Neftci, Hesham Mostafa, and Friedemann Zenke. 2019\.  Surrogate gradient learning in spiking neural networks: Bringing the power of gradient-based optimization to spiking neural networks.  _IEEE Signal Processing Magazine_ 36, 6 (2019), 51–63.
  * Qiao et al. (2015) Ning Qiao, Hesham Mostafa, Federico Corradi, Marc Osswald, Fabio Stefanini, Dora Sumislawska, and Giacomo Indiveri. 2015.  A reconfigurable on-line learning spiking neuromorphic processor comprising 256 neurons and 128K synapses.  _Frontiers in neuroscience_ 9 (2015), 141.
  * Rao and Ballard (1999a) Rajesh PN Rao and Dana H Ballard. 1999a.  Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects.  _Nature neuroscience_ 2, 1 (1999), 79–87.
  * Rao and Ballard (1999b) Rajesh PN Rao and Dana H Ballard. 1999b.  Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects.  _Nature neuroscience_ 2, 1 (1999), 79–87.
  * Song et al. (2020) Yuhang Song, Thomas Lukasiewicz, Zhenghua Xu, and Rafal Bogacz. 2020.  Can the brain do backpropagation? exact implementation of backpropagation in predictive coding networks.  _Advances in neural information processing systems_ 33 (2020), 22566\.
  * Spratling (2019) MW Spratling. 2019\.  Fitting predictive coding to the neurophysiological data.  _Brain research_ 1720 (2019), 146313.
  * Spratling (2008) Michael W Spratling. 2008\.  Predictive coding as a model of biased competition in visual attention.  _Vision research_ 48, 12 (2008), 1391–1408.
  * Spratling et al. (2009) Michael W Spratling, Kris De Meyer, and R Kompass. 2009\.  Unsupervised learning of overlapping image components using divisive input modulation.  _Computational intelligence and neuroscience_ 2009 (2009).
  * Wacongne et al. (2012) Catherine Wacongne, Jean-Pierre Changeux, and Stanislas Dehaene. 2012.  A neuronal model of predictive coding accounting for the mismatch negativity.  _Journal of Neuroscience_ 32, 11 (2012), 3665–3678.
  * Walsh et al. (2020) Kevin S Walsh, David P McGovern, Andy Clark, and Redmond G O’Connell. 2020.  Evaluating the neurophysiological evidence for predictive processing as a model of perception.  _Annals of the new York Academy of Sciences_ 1464, 1 (2020), 242\.
  * Whittington and Bogacz (2017) James CR Whittington and Rafal Bogacz. 2017.  An approximation of the error backpropagation algorithm in a predictive coding network with local hebbian synaptic plasticity.  _Neural computation_ 29, 5 (2017), 1229–1262.
  * Xiao et al. (2018) Will Xiao, Honglin Chen, Qianli Liao, and Tomaso Poggio. 2018\.  Biologically-plausible learning algorithms can scale to large datasets.  _arXiv preprint arXiv:1811.03567_ (2018).

##  Appendix A Research Methods

###  A.1. Division Encoding Cost Function

The division encoding error $e^{**}_{l}$ is equal to
$\sqrt{\frac{a_{l}}{f(p_{l})}}$. Under this encoding scheme, $a_{l}=p_{l}$
when $e^{**}_{l}=1$. The goal then is to update weights and activities to
reduce the difference between $e^{**}_{l}$ and 1. To do so, we use the
following cost function:

(8) $C=\sum_{l}^{L}C_{l},\text{ where }C_{l}=\frac{1}{2}\log(e^{**}_{l})^{2}.$

When $e^{**}_{l}=1$ the log of $e^{**}_{l}$ will equal 0 so the cost at level
$l$ will equal 0. Deviations from 1, will lead the cost to increase. We use
the log here, instead of simply subtracting $e^{**}_{l}$ from 1 because it
simplifies the computations of the gradients.

####  A.1.1. Weight Updates for Division Encoding

Here we derive the gradients of C at level $l$ with respect to the forward
weights at level $l$:

$\begin{split}C_{l}&=\frac{1}{2}(\log(e^{**}_{l}))^{2},\\\
&=\frac{1}{2}(\frac{1}{2}(\log(a_{l})-\log(f(p_{l}))))^{2}.\end{split}$

Lets call the term inside the parentheses $u_{l}$, such that
$C_{l}=\frac{1}{2}u_{l}^{2}$. Now we need to compute $\frac{\partial
C_{l}}{\partial W_{l-1}}$, which can be decomposed using the chain rule as
follows.

$\frac{\partial C_{l}}{\partial W_{l}}=\frac{\partial C_{l}}{\partial
u_{l}}\frac{\partial u_{l}}{\partial f(p_{l})}\,\frac{\partial
f(p_{l})}{\partial p_{l}}\frac{\partial p_{l}}{\partial W_{l-1}}.$

Each term can then be computed individually.

$\begin{split}\frac{\partial C_{l}}{\partial
u_{l}}&=u_{l}=\frac{1}{2}(\log(a_{l})-\log(f(p_{l}))),\\\ \frac{\partial
u_{l}}{\partial f(p_{1})}&=-1/f(p_{l}),\\\ \frac{\partial f(p_{l})}{\partial
p_{l}}&=f^{\prime}(p_{l}),\\\ \frac{\partial p_{l}}{\partial
W_{l-1}}&=a_{l-1}^{T}.\end{split}$

Now if we combine terms we get the following weight update

(9) $\begin{split}\frac{\partial C_{l}}{\partial
W_{l}}&=\frac{1}{2}(\log(a_{l})-\log(p_{l}))(\frac{-1}{f(p_{l})})f^{\prime}(p_{l})a_{l-1}^{T},\\\
&=-\log(\sqrt{\frac{a_{l}}{f(p_{l})}})\frac{f^{\prime}(p_{l})}{f(p_{l})}a_{l-1}^{T},\\\
&=-\log(e_{l}^{**})\frac{f^{\prime}(p_{l})}{f(p_{l})}a_{l-1}^{T}.\end{split}$

We can see this update is the outer product of pre-synaptic and post-synaptic
information. In particular, the pre-synaptic activity is multiplied by the
error neuron activities, which are first passed through a non-linearity
(logarithm) and multiplied by information about the predictions at the same
level.

####  A.1.2. Activity Updates for Division Encoding

Here we derive the gradients of the cost at level $l$ w.r.t. the activities.
Like original equation, the gradients w.r.t. $a_{l}$ are derived from the cost
at the same level and the next level (i.e. level $l$ and $l+1$). The gradient
of $C_{l+1}$ w.r.t $a_{l}$ can be decomposed using the chain rule as follows.

$\frac{\partial C_{l+1}}{\partial a_{l}}=\frac{\partial C_{l+1}}{\partial
u_{l+1}}\frac{\partial u_{l+1}}{\partial f(p_{l+1})}\frac{\partial
f(p_{l+1})}{\partial p_{l+1}}\frac{\partial p_{l+1}}{\partial a_{l}}.$

We already computed all of these terms in the last section except for
$\frac{\partial p_{l}}{\partial a_{l-1}}$. The gradient $\frac{\partial
p_{l}}{\partial a_{l-1}}$ is just equal to $W^{T}_{l}$. This gives us the
bottom-up error (the first term in equation 4) for the activity update
$a_{l}$:

$\frac{\partial C_{l+1}}{\partial
a_{l}}=W_{l}^{T}\log(e_{l}^{**})\frac{f^{\prime}(p_{l})}{f(p_{l})}.$

Now we need to compute $\frac{\partial C_{l}}{\partial a_{l}}$. This term can
be decomposed using the chain rule as follows.

$\frac{\partial C_{l}}{\partial a_{l}}=\frac{\partial C_{l}}{\partial
u_{l}}\frac{\partial u_{l}}{\partial a_{l}}.$

We saw in the last section that $\frac{\partial C_{l}}{\partial u_{l}}$ is
equal to $\log(e_{l}^{**})$. $\frac{\partial u_{l}}{\partial a_{l}}$ is
$\frac{1}{a_{l}}$. With these terms we can now compute the top-down error for
the activity update (second term in equation 4):

$\frac{\partial C_{l}}{\partial a_{l}}=\log(e_{l}^{**})\frac{1}{a_{l}}$

Combining the top-down and bottom up error now gets us the activity updates
under the division encoding schema:

(10) $\Delta
a_{l}=W_{l}^{T}\log(e_{l}^{**})\frac{f^{\prime}(p_{l})}{f(p_{l})}-\log(e_{l}^{**})\frac{1}{a_{l}}$