Benchmarking Predictive Coding Networks -- Made Simple

## Abstract

Abstract In this work, we tackle the problems of efficiency and scalability for predictive coding networks in machine learning. To do so, we first propose a library called PCX, whose focus lies on performance and simplicity, and provides a user-friendly, deep-learning oriented interface. Second, we use PCX to implement a large set of benchmarks for the community to use for their experiments. As most works propose their own tasks and architectures, do not compare one against each other, and focus on small-scale tasks, a simple and fast open-source library adopted by the whole community would address all of these concerns. Third, we perform extensive benchmarks using multiple algorithms, setting new state-of-the-art results in multiple tasks and datasets, as well as highlighting limitations inherent to PC that should be addressed. Thanks to the efficiency of PCX, we are able to analyze larger architectures than commonly used, providing baselines to galvanize community efforts towards one of the main open problems in the field: scalability. The code for PCX is available at https://github.com/liukidar/pcax .

## 1 Introduction

The history of predictive coding is long, and spans a large number of disciplines (Friston, ; Spratling, ).
It first appeared as a computational framework in the 50’s, when electronic engineers realized that sending compressed representations of prediction errors in time series data was cheaper than sending the data itself (Elias, ).
A similar algorithm was then used in the neurosciences, first to describe inhibitory signals in the retina (Srinivasan et al., ), and then as a general theory of information processing in different brain regions (Mumford, ).
In 1999, Rao and Ballard () proposed a hierarchical formulation of predictive coding (PC) as a model of visual processing. Recently, researchers realized that this framework could be used to train neural networks using a bio-plausible learning rule (Whittington and Bogacz, ). This has led to different directions of research that either explored interesting properties of PC networks, such as their robustness (Song et al., ; Alonso et al., ) and flexibility (Salvatori et al., ), or proposed variations to improve the performance on specific tasks (Salvatori et al., ).
While interesting and important for the progress of the field, these lines of research have the tendency of not comparing their results against other papers or those of related fields, and to focus on small-scale experiments. The field is hence avoiding what we believe to be the most important open problem: scaling such results to large scale tasks.

There are multiple reasons why such an important problem has been overlooked.
First, it is a hard problem, and it is still unclear why PC is able to perform as well as classical gradient descent with backprop only up to a certain scale, which is of small convolutional models trained to classify the CIFAR10 dataset (Salvatori et al., ). Understanding the reason for this limitation would allow us to develop regularization techniques that stabilize learning, and hence allow better performance on more complex tasks, similarly to what dropout and batch normalization have been for deep learning (Srivastava et al., ; Ioffe and Szegedy, ). Second, the lack of specialized libraries makes PC models extremely slow: a full hyperparameter search on a small convolutional network can take several hours. Third, the lack of a common framework makes reproducibility and iterative contributions impossible, as implementation details or code are rarely provided.
In this work, we make first steps toward addressing these problems with three contributions, that we call *tool*, *benchmarking*, and *analysis*.

#### Tool.

We release an open-source library for accelerated training for predictive coding called PCX. This library runs in JAX (Bradbury et al., ), and offers a user-friendly interface with minimal learning curve through familiar syntax inspired from Pytorch, and extensive tutorials. It is also fully compatible with Equinox (Kidger and Garcia, ), a popular deep-learning-oriented extension of JAX, ensuring reliability, extendability, and compatibility with ongoing research developments. It also supports JAX’s Just-In-Time (JIT) compilation, making it efficient and allowing both easy development and execution of PC networks. We empirically show the gain in efficiency with respect to an existing library.

#### Benchmarking.

We propose a uniform set of tasks, datasets, metrics, and architectures that should be used as a skeleton to test the performance of future variations of PC. The tasks that we propose are the standard ones in computer vision: image classification and generation (with relative metrics). The models that we use, as well as the datasets, are picked according to two criteria: First, to allow researchers to test their algorithm from the easiest model (feedforward network on MNIST) to more complex ones (deep convolutional models), where we failed to get acceptable results, and should hence pave the way for future research; Second, to favor the comparison against related fields in the literature, such as equilibrium propagation (Scellier and Bengio, ). To this end, we have picked some of the models that are consistently used in other research papers. All the source files of the proposed benchmarks, as well as tutorials explaining how to implement them, will be present in the library, to facilitate researchers to use them in their studies.

#### Analysis.

We provide the baselines for future research by performing an extensive comparative study between different hyperparameters and PC algorithms on multiple tasks. We considered standard PC, incremental PC (Salvatori et al., ), PC with Langevin dynamics (Oliviers et al., ; Zahid et al., ), and nudged PC, as done in the Eqprop literature (Scellier and Bengio, ; Scellier et al., ). This is also the first time nudging algorithms were applied in PC models. In terms of quantitative contributions, we get state-of-the-art results for PC on multiple benchmarks and show for the first time that it is able to perform well on more complex datasets, such as CIFAR100 and Tiny Imagenet, where we get results comparable to those of backprop. In image generation tasks, we present experiments on datasets of colored images, going beyond MNIST and FashionMNIST as performed in previous works. We conclude with an analysis on the credit assignment of PC, which tries to shed light on problems that we will need to solve to scale up such models even more. To this end, we believe that future effort should aim towards developing algorithms and models that improve the numbers that we show in this work, as they represent the new state of the art in the field.

## 2 Related Works

#### Rao and Ballard’s PC.

The most related works are those that explore different properties or optimization algorithms of standard PC in the deep learning regime, using formulations inspired by Rao and Ballard’s original work (Rao and Ballard, ). Examples are works that study their associative memory capabilities (Salvatori et al., ; Yoo and Wood, ; Tang et al., , ), their ability to train Bayesian networks (Salvatori et al., , ), and theoretical results that explain, or improve, their optimization process (Millidge et al., , ; Alonso et al., ). Results in this field have allowed to either improve the performance of such models in different tasks, or to study different properties that could benefit from the use of PCNs.

#### Variations of PC.

In the literature, there are multiple variations of PC algorithms, which differ from Rao and Ballard’s original formulation in the way they update the neural activities. Important examples of such variations are biased competition and divisive input modulation (Spratling, ), or the neural generative coding framework (Ororbia and Kifer, ).
The latter is already used in multiple reinforcement learning and control tasks (Ororbia and Mali, ; Ororbia et al., ). For a review on how different PC algorithms evolved through time, from signal processing to neuroscience, we refer to (Spratling, ); for a more recent review specific to machine learning applications, to (Salvatori et al., ). It is also worth mentioning the original literature on PC in the neurosciences, that does not intersect with ours as it is not related to deep neural networks, has evolved from Rao and Ballard’s work into a general theory that models information processing in the brain using probability and variational inference, called the *free energy principle* (Friston, ; Friston and Kiebel, ; Friston, ).

#### Neuroscience-inspired deep learning.

Another line of related works is that of neuroscience methods applied to machine learning, like equilibrium propagation (Scellier and Bengio, ), which is the most similar to PC (Laborieux and Zenke, ; Millidge et al., ). Other methods able to train models of similar sizes are target propagation (Bengio, ; Ernoult et al., ; Millidge et al., ) and SoftHebb (Moraitis et al., ; Journé et al., ). The first two communities, that of targetprop and eqprop, consistently use similar architectures in different research papers to test the performance of their methods. In our benchmarking effort, some of the architectures proposed are the same ones they use, to favor a more direct comparison. There are also methods that differ more from PC, such as forward-only methods (Kohan et al., ; Nøkland, ; Hinton, ), methods that back-propagate the errors using a designated set of weights (Lillicrap et al., ; Launay et al., ), and other Hebbian methods (Moraitis et al., ; Journé et al., ).

## 3 Background and Notation

Predictive coding networks (PCNs) are hierarchical Gaussian generative models with $L$ levels with parameters $\theta=\{\theta_{0},\theta_{1},\theta_{2},...,\theta_{L}\}$, in which each level models a multi-variate distribution parameterized by activation of the preceding level. Let $h_{l}$ be the realization of the vector of random variables $H_{l}$
of level $l$, then we have that the likelihood

$$ $P_{\theta}(h_{0},h_{1},\ldots,h_{L})=P_{\theta_{0}}(h_{0})P_{\theta_{1}}(h_{1} |h_{0})\cdots P_{\theta_{L}}(h_{L}|h_{L-1}).$ $$

For simplicity, we write $P_{\theta_{l}}(h_{l})$ instead of $P_{\theta_{l}}(H_{l}=h_{l})$. We refer to each of the scalar random variables of $H_{l}$ as a neuron. PC assumes that both the prior on $h_{0}$ and the relationships between levels are governed by a normal distribution parameterized as follows:

$$ $\displaystyle P_{\theta_{0}}(h_{0})=\mathcal{N}(h_{0},\mu_{0},\Sigma_{0}),\;\; \;\;\;\;\;\;\,\mu_{0}=\theta_{0},$ $\displaystyle P_{\theta_{l}}(h_{l}|h_{l-1})=\mathcal{N}(h_{l};\mu_{l},\Sigma_{ l}),\;\;\;\mu_{l}=f_{l}(h_{l-1},\theta_{l}),$ $$

where $\theta_{l}$ are the learnable weights parametrizing the transformation $f_{l}$, and $\Sigma_{l}$ is a covariance matrix. As it is standard practice, $\Sigma_{l}$ will be fixed to the identity matrix throughout this work (Whittington and Bogacz, ). If, for example, $\theta_{l}=(W_{l},b_{l})$ and $f_{l}(h_{l-1},\theta_{l})=\sigma_{l}(W_{l}h_{l-1}+b_{l})$, then the neurons in level $l-1$ are connected to neurons in level $l$ via a linear operation, followed by a non linear map, that is the analogous to a fully connected layer in standard deep learning. Intuitively, $\theta$ is the set of learnable weights of the model, while $h=\{h_{0},h_{1},...,h_{L}\}$ is data-point-dependent latent state, containing the abstract representations for the given observations.

#### Training.

In supervised settings, training consists of learning the relationship between given pairs of input-output observations $(x,y)$.
In PC, this is performed by maximizing the joint likelihood of our generative model with the latent vectors $h_{0}$ and $h_{L}$ respectively fixed to the input and label of the provided data-point: $P_{\theta}(h\rvert_{h_{0}=x,h_{L}=y})=P_{\theta}(h_{L}=y,\ldots,h_{1},h_{0}=x)$.
This is achieved by minimizing the so-called variational free energy $\mathcal{F}$ (Friston et al., ):

$$ $\mathcal{F}(h,\theta)=-\ln{P_{\theta}(h)}=-\ln{\left(\mathcal{N}(h_{0}|\mu_{0} )\prod_{l=1}^{L}\mathcal{N}(h_{l};f_{l}(h_{l-1},\theta_{l}))\right)}=\sum_{l=0 }^{L}\frac{1}{2}(h_{l}-\mu_{l})^{2}+k.$ (1) $$

The quantity $\epsilon_{l}=(h_{l}-\mu_{l})$ is often referred to as prediction error of layer $l$, being the difference between the predicted activation $\mu_{l}$ and the current state $h_{l}$. Refer to the appendix, for a full derivation of Eq. (1).
To minimize $\mathcal{F}$, the Expectation-Maximization (EM) (Dempster et al., ) algorithm is used by iteratively optimizing first the state $h$, and then the weights $\theta$ according to the equations

$$ $h^{*}=\underset{h}{\arg\min}\,\mathcal{F}(h,\theta),\;\;\theta^{*}=\underset{ \theta}{\arg\min}\,\mathcal{F}(h^{*},\theta).$ (2) $$

We refer to the first step described by Eq. (2) as inference and to the second as learning phase.
In practice, we do not train on a single pair $(x,y)$ but on a dataset split in mini-batches that are subsequently used to train the model parameters. Furthermore, both inference and learning are approximated via gradient descent on the variational free energy.
In the inference phase, firstly $h$ is initialized to an initial value $h^{(0)}$, and, then, it is optimized for $T$ iterations. Then, during the learning phase we use the newly computed values to perform a single update on the weights $\theta$. The gradients of the variational free energy with respect to both $h$ and $\theta$ are as follows:

$$ $\nabla h_{l}=\frac{\partial\mathcal{F}}{\partial h_{l}}=\frac{1}{2}\left(\frac {\partial\epsilon_{l}^{2}}{\partial h_{l}}+\frac{\partial\epsilon_{l+1}^{2}}{ \partial h_{l}}\right),\;\;\;\;\;\;\nabla\theta_{l}=\frac{\partial\mathcal{F}} {\partial\theta_{l}}=\frac{1}{2}\frac{\partial\epsilon_{l}^{2}}{\partial\theta _{l}}.$ (3) $$

Then, a new batch of data points is provided to the model and the process is repeated until convergence. As highlighted by Eq. (3), each state and each parameter is updated using local information as the gradients depend exclusively on the pre and post-synaptic errors $\epsilon_{l}$ and $\epsilon_{l+1}$. This is the main reason why, in contrast to BP, PC is a local algorithm and is considered more biologically plausible.
In Appendix A, we provide an algorithmic description of the concepts illustrated in these paragraphs, highlighting how each equation is translated to code in PCX.

Figure: Figure 1: Left: Generative and discriminative modes, Right: Pseudocode of one parameter update of PC in supervised learning, as well as an informal description of the different algorithms considered in this work.
Refer to caption: x1.png

#### Evaluation.

This phase is similar to the inference phase, with the difference that we perform it on a test point $\bar{x}$, used to infer the label $\bar{y}$. This is achieved by fixing $h_{0}=\bar{x}$ and compute the most likely value of the latent states $h^{*}\rvert_{h_{0}=\bar{x}}$, again using the state gradients of Eq. (3). We refer to this as discriminative mode. In practice, for discriminiative networks, the values of the latent states computed this way are equivalent to those obtained via a forward pass, that is setting $h_{l}^{(0)}=\mu_{l}^{(0)}$ for every $l\neq 0$, as it corresponds to the global minimum of $\mathcal{F}$ (Frieder and Lukasiewicz, ).

#### Generative Mode.

So far, we have only described how to use PCNs to perform supervised training. However, as we will see in the experimental section, such models can also be used (and were initially developed to be used) to perform unsupervised learning tasks. Given a datapoint $x$, the goal is to use PCNs to compress the information of $x$ into a latent representation, conceptually similar to how variational autoencoders work (Kingma and Welling, ). Such a compression, that should contain all the information needed to generate $x$, is computed by fixing the state vector $h_{L}$ to the data-point, and run inference – that is, we maximize $P_{\theta}(h\rvert_{h_{L}=x})$ via gradient descent on $h$ – until the process has converged. The compressed representation will then be the value of $h_{0}$ at convergence (or, in practice, after $T$ steps). If we are training the model, we then perform a gradient update on the parameters to minimize the variational free energy of Eq.(1), as we do in supervised learning. A sketch of the discriminative and generative ways of training PCNs is represented in Fig. 1(a).

## 4 Experiments and Benchmarks

This section is divided in two areas, that correspond to discriminative and generative inference tasks. The first focuses on classification tasks; the second on unsupervised generation. A sketch illustrating the two modes is provided in Fig. 1. To provide a comprehensive evaluation, we will test our models on multiple computer vision datasets, MNIST, FashionMNIST, CIFAR10/100, CelebA, and Tiny ImageNET; on models of increasing complexity, with both feedforward and convolutional layers; and multiple learning algorithms present in the literature.

#### Algorithms.

We consider various learning algorithms present in the literature: (1) Standard PC, already discussed in the background section; (2) Incremental PC (iPC), a simple and recently proposed modification where the weight parameters are updated alongside the latent variables at every time step; (3) Monte Carlo PC (MCPC), obtained by applying unadjusted Langevin dynamics to the inference process; (4) Positive nudging (PN), where the target used is obtained by a small perturbation of the output towards the original, 1-hot label; (5) Negative nudging (NN), where the target is obtained by a small perturbation *away* from the target, and updating the weights in the opposite direction; (6) Centered nudging (CN), where we alternate epochs of positive and negative nudging.
Among these, PC, iPC, and MCPC will be used for the generative mode, and PC, iPC, PN, and NN, and CN for the discriminative mode. See Fig. 1, and the supplementary material, for a more detailed description.

### 4.1 Discriminative Mode

**Table 1: Test accuracies of the different algorithms on different datasets.**
| % Accuracy | PC-CE | PC-SE | PN | NN | CN | iPC | BP-CE | BP-SE |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| MLP |  |  |  |  |  |  |  |  |
| MNIST | $98.11^{\pm 0.03}$ | $98.26^{\pm 0.04}$ | $98.36^{\pm 0.06}$ | $98.26^{\pm 0.07}$ | $98.23^{\pm 0.09}$ | $\mathbf{98.45^{\pm 0.09}}$ | $98.07^{\pm 0.06}$ | $98.29^{\pm 0.08}$ |
| FashionMNIST | $89.16^{\pm 0.08}$ | $89.58^{\pm 0.13}$ | $89.57^{\pm 0.08}$ | $89.46^{\pm 0.08}$ | $89.56^{\pm 0.05}$ | $\mathbf{89.90^{\pm 0.06}}$ | $89.04^{\pm 0.08}$ | $89.48^{\pm 0.07}$ |
| VGG-5 |  |  |  |  |  |  |  |  |
| CIFAR-10 | $88.06^{\pm 0.13}$ | $87.98^{\pm 0.11}$ | $88.42^{\pm 0.66}$ | $88.83^{\pm 0.04}$ | $\mathbf{89.47^{\pm 0.13}}$ | $85.51^{\pm 0.12}$ | $88.11^{\pm 0.13}$ | $89.43^{\pm 0.12}$ |
| CIFAR-100 (Top-1) | $60.00^{\pm 0.19}$ | $54.08^{\pm 1.66}$ | $64.70^{\pm 0.25}$ | $65.46^{\pm 0.05}$ | $\mathbf{67.19^{\pm 0.24}}$ | $56.07^{\pm 0.16}$ | $60.82^{\pm 0.10}$ | $66.28^{\pm 0.23}$ |
| CIFAR-100 (Top-5) | $84.97^{\pm 0.19}$ | $78.70^{\pm 1.00}$ | $84.74^{\pm 0.38}$ | $85.15^{\pm 0.16}$ | $\mathbf{86.60^{\pm 0.18}}$ | $78.91^{\pm 0.23}$ | $85.84^{\pm 0.14}$ | $85.85^{\pm 0.27}$ |
| Tiny ImageNet (Top-1) | $41.29^{\pm 0.2}$ | $30.28^{\pm 0.2}$ | $34.61^{\pm 0.2}$ | $\mathbf{46.40^{\pm 0.1}}$ | $46.38^{\pm 0.11}$ | $29.94^{\pm 0.47}$ | $43.72^{\pm 0.1}$ | $44.90^{\pm 0.2}$ |
| Tiny ImageNet (Top-5) | $66.68^{\pm 0.09}$ | $57.31^{\pm 0.21}$ | $59.91^{\pm 0.24}$ | $68.50^{\pm 0.18}$ | $69.06^{\pm 0.10}$ | $54.73^{\pm 0.52}$ | $\mathbf{69.23^{\pm 0.23}}$ | $65.26^{\pm 0.37}$ |
| VGG-7 |  |  |  |  |  |  |  |  |
| CIFAR-100 (Top-1) | $56.80^{\pm 0.14}$ | $37.52^{\pm 2.60}$ | $56.56^{\pm 0.13}$ | $59.97^{\pm 0.41}$ | $64.76^{\pm 0.17}$ | $43.99^{\pm 0.30}$ | $59.96^{\pm 0.10}$ | $\mathbf{65.36^{\pm 0.15}}$ |
| CIFAR-100 (Top-5) | $83.00^{\pm 0.09}$ | $66.73^{\pm 2.37}$ | $81.52^{\pm 0.17}$ | $81.50^{\pm 0.41}$ | $84.65^{\pm 0.18}$ | $73.23^{\pm 0.30}$ | $\mathbf{85.61^{\pm 0.10}}$ | $84.41^{\pm 0.26}$ |
| Tiny ImageNet (Top-1) | $41.15^{\pm 0.14}$ | $21.28^{\pm 0.46}$ | $25.53^{\pm 0.77}$ | $39.49^{\pm 2.69}$ | $35.59^{\pm 7.69}$ | $19.76^{\pm 0.15}$ | $45.32^{\pm 0.11}$ | $\mathbf{46.08^{\pm 0.15}}$ |
| Tiny ImageNet (Top-5) | $66.25^{\pm 0.11}$ | $44.92^{\pm 0.27}$ | $50.06^{\pm 0.84}$ | $64.66^{\pm 1.95}$ | $59.63^{\pm 6.00}$ | $40.36^{\pm 0.22}$ | $\mathbf{69.64^{\pm 0.18}}$ | $66.65^{\pm 0.20}$ |

Here, we test the performance of PCNs on image classification tasks. We compare PC against BP, using both Squared Error (SE) and Cross Entropy (CE) loss, by adapting the energy function as described in (Pinchetti et al., ). For the experiments on MNIST and FashionMNIST, we use feedforward models with $3$ hidden layers of $128$ hidden neurons, while for CIFAR10/100 and Tiny ImageNET, we compare VGG-like models. We performed a large hyperparameter search over learning rates, optimizers, activation functions, and algorithm-specific parameters. All the details needed to reproduce the experiments, as well as a discussion about ‘lessons learned’ during such a large search, are in the Appendix B. Results, averaged over $5$ seeds are reported in Tab. 1.

#### Discussion.

The results show that the best performing algorithms, at least on the most complex tasks, are the ones where the target is nudged towards the real label, that are PN and NN. This is in line with previous findings in the Eqprop literature (Scellier et al., ). The recently proposed iPC, on the other hand, performs well on small architectures, as it is the best performing one on MNIST and FashionMNIST, but its performance worsten when it comes to train on large architectures. More broadly, the performance are comparable to those of backprop, except on the largest model. An interesting observation, is that all the best results for PC have been achieved using a VGG5, that has always outperformed the deeper VGG7. Future work should investigate the reason of such a phenomenon, as scaling up to more complex datasets will require the use of much deeper architectures, such as ResNets (He et al., ). In Section 5, we analyze possible causes, as well as comparing the wall-clock time of the different algorithms.

**Table 2: MSE loss for image reconstruction of BP, PC, and iPC on different datasets.**
| MSE ($\times 10^{-3}$) | PC | iPC | BP |
| --- | --- | --- | --- |
| MNIST | $9.25^{\pm 0.00}$ | $9.09^{\pm 0.00}$ | $\mathbf{9.08^{\pm 0.00}}$ |
| FashionMNIST | $10.56^{\pm 0.01}$ | $10.11^{\pm 0.01}$ | $\mathbf{10.04^{\pm 0.00}}$ |

Table: Figure 2: CIFAR10 image reconstruction via autoencoding convolutional networks. In order: original, PC, iPC, BP, BP (half of the parameters).

### 4.2 Generative Mode

We test the performance of PCNs on image generation tasks. We perform three different kinds of experiments: (1) generation from a posterior distribution; (2) generation via sampling from the learned joint distribution; and (3) associative memory retrieval. In the first case, we provide a test image $y$ to a trained model, run inference to compute a compressed representation $\bar{x}$ (stored in the latent vector $h_{0}$ at convergence), and produce a reconstructed $\bar{y}=h_{L}$ by performing a forward pass with $h_{0}=\bar{x}$). The model we consider are three layer networks. As this is an autoencoding task, we compare against autoencoders with three layer encoder/decoder structure (so, six layers in total). In the case of MNIST and FashionMNSIT, we use feedforward layers, in the case of CIFAR10 and CelebA, deconvolutional (and convolutional for the encoder) ones. The results in Tab. 2 and Fig. 2 report comparable performance, with a small advantage for PC compared to BP on the more complex convolutional tasks. In this case, iPC is the best performing algorithm, probably due to the small size of the considered models which allows for better stability. Furthermore, note that the BP architectures have double the amount of parameters, being the PC networks decoder only. If we halve the number of features in each layer of the autoencoder architecture (while keeping the bottleneck dimension unchanged), we get significantly reduced performance for BP (Fig. 2, bottom), achieving a loss of $10.66^{\pm{0.94}}\times 10^{-3}$ on CIFAR10. Details about these and the following experiments are provided in Appendix C.

**Table 3: MSE ($\times 10^{-4}$) of associative memory tasks. Columns indicate the number of hidden neurons while rows shows the training images to memorize. Results over 5 seeds.**
| Noise | 512 | 1024 | 2048 |
| --- | --- | --- | --- |
| $50$ | $6.06^{\pm 0.11}$ | $5.91^{\pm 0.14}$ | $5.95^{\pm 0.06}$ |
| $100$ | $6.99^{\pm 0.19}$ | $6.76^{\pm 0.23}$ | $6.16^{\pm 0.07}$ |
| $250$ | $9.95^{\pm 0.05}$ | $10.14^{\pm 0.06}$ | $8.90^{\pm 0.06}$ |

Figure: Figure 4: Memory recalled images. Top: Original images. Left: Noisy input (guassian noise, $\sigma=0.2$) and reconstruction. Right: Masked input (bottom half removed) and reconstruction.
Refer to caption: x5.png

For the second category of experiments, we tested the capability of PCNs to learn, and sample from, a complex probability distribution.
MCPC extends PC by incorporating Gaussian noise to the activity updates of each neuron. This change enables a PCN to learn and generate samples analogous to a variational autoencoders (VAE). This change shifts the inference of PCNs from a variational approximation to Monte Carlo sampling of the posterior using Langevin dynamics. Data samples can be generated from the learned joint $P_{\theta}(h)$ by leaving all states $h_{l}$ free and performing noisy inference updates. Figure 3 illustrates MCPC’s ability to learn non-linear multimodal distributions using the iris dataset (Pedregosa et al., ) and shows generative samples for MNIST. When comparing MCPC to a VAE on MNIST, both models produced samples of similar quality despite the VAE having twice the number of parameters. MCPC achieved a lower FID score (MCPC: $2.53^{\pm 0.17}$ vs. VAE: $4.19^{\pm 0.38}$), whereas the VAE attained a higher inception score (VAE: $7.91^{\pm 0.03}$ vs. MCPC: $7.13^{\pm 0.10}$).

In the associative memory (AM) experiments, we test how well the model is able to reconstruct an image already present in the training set, after it is provided with an incomplete or corrupted version of it, as done in a previous work (Salvatori et al., ). Fig. 4 show the results obtained by a PCN with 2 hidden layers of 512 neurons given noise or mask corrupted images. In Tab. 3, we study the memory capacity as the number of hidden layers increases. No visual difference between the recall and original images can be observed for MSE up to $0.005$. To evaluate efficiency we then trained a PCN with 5 hidden layers of 512 neurons on 500 TinyImagenet samples, with a batch size of 50 and 50 inference iterations during training. Training takes $0.40\pm 0.005$ seconds per epoch on an Nvidia V100 GPU.

#### Discussion.

The results show that PC is able to perform generative tasks, as well as and associative memory ones using decoder-only architectures. Via inference, PCNs are able to encode complex probability distributions in their latent state which can be used to perform a variety of different tasks, as we have shown. Thus, compared to artificial neural networks, PCNs are more flexible and require only half the parameters to achieve similar performance. This comes at a higher computational cost due to the number of inference steps to perform. Future work should look into this issue, aiming at reducing the inference time by propagating the information more efficiently through the network.

## 5 Analysis and metrics

In this section, we report several metrics that we believe are important to understand the current state and challenges of training networks with PC and compare them with standard models trained with gradient-descent and backprop when suitable. More in detail, we discuss how regularly the energy flows into the model, and how stable training is when changing parameters, initializations, and optimizers. A better understanding of such phenomena would allow us to solve the current problems of PCNs and, hence, scale up to the training of larger models on more complex datasets.

### 5.1 Energy and stability

The first study we perform regards the initialization of the network states $h$, and how this influences the performance of the model. In the literature, they have been either initialized to be equal to *zero*, randomly initialized via a Gaussian prior (Whittington and Bogacz, ), and be initialized via a forward pass. This last technique has been the preferred option in machine learning papers as it sets the errors $\mathcal{\epsilon}_{l\neq L}=0$ at every internal layer of the model. This allows the prediction error to be concentrated in the output layer only, and hence be equivalent to the SE.
To provide a comparison among the three methods, we have trained a $3$-layer feedforward model on FashionMNIST. The results, plotted in Fig. 5(a), show that forward initialization is indeed the better method, although the gap in performance shrinks the more iterations $T$ are performed.

Figure: (a)
Refer to caption: x6.png

#### Energy propagation.

Concentrating the total error of the model, and hence its energy, to the last layer as done when performing forward initialization, makes it hard for the model to then propagate such an energy back to the first layers. As reported in Fig. 5(b), we observe that the energy in the last layer is orders of magnitude larger than the one in the input layer, even after performing several inference steps. However, this behavior raises the question whether better initialization or optimization techniques could result in a more balanced energy distribution and thus better weight updates, as learning in this unbalanced energy regime has been shown problematic for more complex models (Alonso et al., ). An easy way of quickly propagating the energy through the network is to use learning rates equal to $1.0$ for the updates of the states. However, both the results reported in Fig. 5(c), as well as our large experimental analysis of Section 4 show that the best performance was consistently achieved for state learning rates $\gamma$ significantly smaller than $1.0$.

We hypothesize that the current training setup for PCNs favors small state learning rates that are sub-optimal to scale to deeper architectures. Fig. 5(d) shows the energy ratios for different state learning rates: when $\gamma\ll 1.0$, the ratios of energies between layers are small, $\nicefrac{{\epsilon_{l+1}^{2}}}{{\epsilon_{l}^{2}}}\ll 1$. The energy in the the first hidden layer is on average 6 orders of magnitude below the last layer for $\gamma=0.01$. While models trained with large $\gamma$ values achieve better energy propagation, they achieve lower accuracy as shown in Fig 5(c). Note that the decay in performance as function of increasing $\gamma$ is stronger for Adam despite being the overall better optimizer in our experiments. This suggests limitations in the current training techniques and possible direction for future improvements aimed at reducing the energy imbalance between layers. We provide implementation details and results on other datasets in Appendix D.

#### Training stability.

Figure: Figure 6: Updating weights with AdamW becomes unstable for wide layers as the accuracy plummets to random guessing for progressively smaller state learning rates as the network’s width increases. Contrarely to using SGD, the optimal state learning rate depends on the width of the layers.
Refer to caption: x10.png

We observed a further link between the weight optimizer and the structure of a PCN that might hinder the scalability of PC, that is, the influence of the hidden dimension on the performance of the model.
To better study this, we trained feedforward PCNs with different hidden dimensions, state learning rates $\gamma$ and optimizers, and reported the results in Fig. 6. The results show that, when using Adam, the width strongly affects the values of the learning rate $\gamma$ for which the training process is stable. Interestingly, this phenomenon does not appear when using both the SGD optimizer, nor on standard networks trained with backprop. This behavioral difference with BP is unexpected and suggests the need for better optimization strategies for PCNs, as Adam was still the best choice in our experiments, but could be a bottleneck for larger architectures.

### 5.2 Properties of predictive coding networks

With PCX, it is straightforward to inspect and analyze several properties of PCNs. Here, we use $\mathcal{F}$ to differentiate between in-distribution (ID) and out-of-distribution (OOD) due to a semantic distribution shift Liu et al. (), as well as to compute the likelihood of a datasets Grathwohl et al. (). This can occur when samples are drawn from different, unseen classes, such as FashionMNIST samples under an MNIST setup Hendrycks and Gimpel ().
To understand the confidence of the predictions of a PCN, we compare the distribution of the probability $p(x,\hat{y};\theta)$, for ID and OOD samples against these to the distribution of softmax values of the PCN classifier, and compute their negative log-likelihoods (NLL), according to the assumptions stated in Section 3, that is,

$$ $\mathcal{F}=-\ln{p(x,y;\theta)}\implies p(x,y;\theta)=e^{-\mathcal{F}}.$ (4) $$

Our results in Fig. 7(a) demonstrate that a trained PCN classifier can effectively (1) assess OOD samples out-of-the-box, without requiring specific training for that purpose Yang et al. (), and (2) produce energy scores for ID and OOD samples that initially correlate with softmax values prior to the optimization of the states variables, $h$. However, after optimizing the states for $T$ inference steps, the scores for ID and OOD samples become decorrelated, especially for samples with lower softmax values as shown in Fig. 7(b).
To corroborate this observation, we also present ROC curves for the most challenging samples, including only the lowest $25\%$ of the scores. As shown in Fig.7(c), the probability (i.e., energy-based) scores provide a more reliable assessment of whether samples are OOD. Experiment details and results on other datasets are provided in in Appendix E.

Figure: (a)
Refer to caption: x11.png

**Table 4: Comparison of the training times of BP against PC on different architectures and datasets.**
| Epoch time (seconds) | BP | PC (ours) | PC (Song) |
| --- | --- | --- | --- |
| MLP - FashionMNIST | $1.82^{\pm 0.01}$ | $1.94^{\pm 0.07}$ | $5.94^{\pm 0.55}$ |
| AlexNet - CIFAR-10 | $1.04^{\pm 0.08}$ | $3.86^{\pm 0.06}$ | $17.93^{\pm 0.37}$ |
| VGG-5 - CIFAR-100 | $1.61^{\pm 0.04}$ | $5.33^{\pm 0.02}$ | $13.49^{\pm 0.05}$ |
| VGG-7 - Tiny ImageNet | $7.59^{\pm 0.63}$ | $54.60^{\pm 0.10}$ | $137.58^{\pm 0.08}$ |

## 6 Computational Resources and Implementations Details

PCX is developed on top on JAX, focusing on performance and versatility. We plan to further improve its efficiency and expand its capabilities to support new developments in the field of backprop-free training techniques that align with the PC core principles. In relation to current alternatives, however, PCX is a very competitive choice. We measured the wall-clock time of our PCNs implementation against another existing open-source library (Song, ) used in many PC works (Song et al., ; Salvatori et al., , ; Tang et al., ), as well as comparing it with equivalent BP-trained networks (developed also with PCX for a fair comparison). Tab. 4 reports the measured time per epoch, averaged over 5 trials, using a A100 GPU. Despite being a sub-optimal architecture in term of classification performance, we also tested on AlexNet (Krizhevsky et al., ), showing that our library can efficiently train models of more than 100 million parameters (having AlexNet $\approx 160$ million parameters). We also outperform alternative methods such as Eqprop: using the same architecture on CIFAR100, the authors report that one epochs takes $\approx 110$ seconds, while we take $\approx 5.5$ on the same hardware (Scellier et al., ). However, we stress that this is not an apple-to-apple comparison, as the authors are more concerned with simulations on analog circuits, rather than achieving optimal GPU usage.

#### Limitations.

The efficiency of PCX could be further increased by fully parallelizing all the computations occuring within a PCN. In fact, in its current state, JIT seems to be unable to parallelize the executions of the layers; a problem that can be addressed with the JAX primitive vmap, but only in the unpractical case where all the layers have the same dimension. To test how different hyperparameters of the model influence the training speed, we have taken a feedforward model, and trained it multiple times, each time increasing a specific hyperparameter by a multiplicative factor. The results, reported in Fig. 8, show that the two parameters that increase the training time are the number of layers $L$ (when the calculations are not explicitly vmapped), and the number of steps $T$. Ideally, only $T$ should affect the training time as inference is an inherently sequential process that cannot be parallelized, but this is not the case, as the time scales linearly with‘F the amount of layers. Details are reported in Appendix F.

## 7 Discussion

The main contribution of this work is the introduction and open-source release that can be used to perform deep learning tasks using PCNs. Its efficiency relies on JAX’s Just-In-Time compilation and carefully structured primitives built to take advantage of it. A second advantage of our library is its intuitive setup, tailored to users already familiar with other deep learning frameworks such as PyTorch. Together with the large number of tutorials we release will make it easy for new users to train networks using PC.

We have also performed a large-scale comparison study on image classification and image generation tasks, unifying under the same computational framework multiple optimization algorithms for PCNs present in the literature. In addition, we have tried multiple, popular, optimizers and training techniques, as well as an extremely large choice of hyperparameters. For CIFAR100 only, for example, we have conducted thousands of individual experiments, that have been used to obtain new state of the art results, as well as provide insights on what works and what does not, that will be useful in the future to researchers tackling deep learning problem with PCNs.

## Acknowledgments

This work was supported by Medical Research Council grant MC_UU_00003/1 awarded to RB.

## Appendix

The code for PCX is available at [https://github.com/liukidar/pcax](https://github.com/liukidar/pcax). The experiments conducted in the main body and in the appendix can be found in the examples folder.

Here we provide the details on how experiments were conducted and results obtained. We opt for a more descriptive approach to convey the fundamental concepts, and leave all details for reproducibility in the provided code, as well as in the next sections. There, each section will link to the exact directory corresponding to the described experiments.

## Appendix A PCX – A Brief Introduction

In this section, we illustrate the core ideas of PCX by describing the main building blocks necessary to train and evaluate a feedforward classifier in predictive coding. For more detailed and complete explanations, please refer to the tutorial notebooks in the examples folder of the library.

In Section 3, we defined PCNs as models with parameters $\theta=\{\theta_{0},\ldots,\theta_{L}\}$ and state $h=\{h_{0},\ldots,h_{L}\}$. In PCX, we divide a model in components of two main categories: layers (i.e., the traditional deep-learning transformations such as ’Linear’ or ’Conv2D’) and vodes (i.e., vectorized nodes that store the array of neurons representing state $h_{l}$). A PCN is defined as follows:
{minted}python
import jax.nn as jnn
import pcx.predictive_coding as pxc
import pcx.nn as pxnn

class MLP(pcx.EnergyModule):
def __init__(self, in_dim, h_dim, out_dim):
self.layers = [
pxnn.Linear(in_dim, h_dim),
pxnn.Linear(h_dim, h_dim),
pxnn.Linear(h_dim, out_dim)
]

self.vodes = [
pxc.Vode((dim,)) for dim in (h_dim, h_dim, out_dim)
]

def __call__(self, x, y = None):
for layer, vode in zip(self.layers, self.vodes):
u = jnn.leaky_relu(layer(x))
x = vode(u)

if y is not None:
self.vodes[-1].set("h", y)

return u
In the __call__ method, we forward the input $x$ through the network. Note that every time we call a vode, we are effectively storing in it the activation $u_{l}$ (so that we can later compute the energy $\epsilon_{l}^{2}$ associated to the vode) and returning its state $h_{l}$ (i.e., x = vode(u) corresponds to vode.set("u", u); x = vode.get("h")). During training, the label $y$ is provided to the model and fixed to the last vode by overwriting its state $h$. Note that, since both during training and evaluation the state of the first vode would be fixed to the input $x$, we avoid defining it (i.e., we avoid computing $P_{\theta_{0}}(h_{0})$ since it would be constant), and directly forward $x$ to the first layer transformation.

The class pxc.EnergyModule provides a .energy() function that computes the variational free energy $\mathcal{F}$ as per Eq. (1). We can compute the state and parameters gradients as per Eqs. (3) by calling pxf.value_and_grad, a wrap around the homonymous JAX function. Having defined two optimizers, optim_w and optim_h, for parameters and state respectively, we can define training on a pair $(x,y)$ as following:
{minted}python
import pcx.utils as pxu
import pcx.functional as pxf

def energy(x, y, *, model):
model(x, y)
return model.energy()

grad_h = pxf.value_and_grad(
pxu.Mask(pxc.VodeParam, [False, True])
)(energy)

grad_w = pxf.value_and_grad(
pxu.Mask(pxc.LayerParam, [False, True])
)(energy)

def train(T, x, y, *, model, optim_h, optim_w):
model.train()

# Initialization
with pxu.step(model, pxc.STATUS.INIT, clear_params=pxc.VodeParam.Cache):
model(x)

# Inference steps
for i in range(T):
with pxu.step(model, clear_params=pxc.VodeParam.Cache):
_, g_h = grad_h(x, y, model=model)
optim_h.step(model, g_h["model"], True)

# Learning step
with pxu.step(model, clear_params=pxc.VodeParam.Cache):
_, g_w = grad_w(x, y, model=model)
optim_w.step(model, g_w["model"])
A few notes on the above code:

- •
JAX [Bradbury et al., ] is a functional library, PCX is not. Modules in PCX are PyTrees, using the same philosophy as another popular JAX library, equinox [Kidger and Garcia, ], with which PCX modules are fully compatible. However, their state is managed by PCX so that each parameter transformation is automatically tracked. The user can opt in for this behavior by passing arguments as keyword argmunets (such as in the above example). Positional function parameters, instead are ignored by PCX and it is the user’s duty to track their state as done in JAX or equinox.
- •
pxf.value_and_grad allows to specify a Mask object to identify which parameters to target with the given transformation. In the case above, we first compute the gradient of $\mathcal{F}$ with respect of the state (VodeParam) and, then, of the weights (LayerParam) of the model.
- •
In the train function, we use pxu.step to set the model status to pxc.STATUS.INIT to perform the state initialization. In PCX, forward initialization is the default method, however other ones can be easily specified. pxu.step is also used to clear the PCN’s cache which is used to store intermediate values such as the activations $u_{l}$.
- •
The actual examples in the library are on mini-batches of data, so all transformations above are vmapped in the actual experiments.

For the evaluation function, being in discriminative mode, we simply perform a forward pass through the PCN which sets $\epsilon_{l}=0$ for all layers.
{minted}python
def eval(x, *, model):
with pxu.step(model, pxc.STATUS.INIT, clear_params=pxc.VodeParam.Cache):
return model(x)

## Appendix B Discriminative experiments

#### Model.

We conducted experiments on three models: MLP, VGG-5, and VGG-7. The detailed architectures of these models are presented in Table 5.

**Table 5: Detailed Architectures of base models**
|  | MLP | VGG-5 | VGG-7 |
| --- | --- | --- | --- |
| Channel Sizes | [128, 128] | [128, 256, 512, 512] | [128, 128, 256, 256, 512, 512] |
| Kernel Sizes | - | [3, 3, 3, 3] | [3, 3, 3, 3, 3, 3] |
| Strides | - | [1, 1, 1, 1] | [1, 1, 1, 1, 1, 1] |
| Paddings | - | [1, 1, 1, 0] | [1, 1, 1, 0, 1, 0] |
| Pool window | - | 2 × 2 | 2 × 2 |
| Pool stride | - | 2 | 2 |

For each model, we conducted experiments with seven different algorithms:

- 1.
Standard PC with Cross-Entropy Loss (PC-CE) / Squared Error Loss (PC-SE): already discussed in the background section.
- 2.
PC with Positive Nudging (PC-PN):
Unlike standard Predictive Coding with Squared Error Loss (PC-SE), where the output is clamped to the target, we “nudge” the output towards the target in PC with nudging. This is achieved by fixing the representation $h$ of last layer $h_{L}$ to $\mu_{L}+\beta(y-\mu_{L})$, where $\mu_{L}$ is the predicted activation of the last layer after forward initialisation, $y$ is the target, and $\beta\in(0,1)$ is a scalar parameter that controls the strength of nudging. Note that when $\beta=1$, PC with nudging is equivalent to the standard PC.
During training procedure, as the model output gradually approaches to the target, we employ a strategy of increasing $\beta$. At the end of each epoch, the value of $\beta$ is incremented by a fixed learning rate $\beta_{ir}$. When $\beta$ becomes greater than or equal to 1, we set it to 1. This strategy allows the model to learn and explore in the early stages of training, while gradually transitioning to standard PC in the later stages.
- 3.
PC with Negative Nudging (PC-NN):
In this algorithm, we do the opposite of positive nudging: we push the output away from the target. Therefore, we fix the representation $h$ of the last layer to $\mu_{L}-\beta(y-\mu_{L})$. We use the same strategy of dynamically increasing $\beta$. When $\beta$ becomes greater than or equal to 1, we set it to -1.
In the learning stage, to ensure that the direction of the weight update is consistent with the target (since we fixed $h_{L}$ to the opposite direction), we invert the weight update: $\theta_{l}\leftarrow\theta_{l}-\Delta\theta_{l}$
where $\Delta\theta_{l}$ defined in the Eq. 3.
- 4.
PC with Center Nudging (PC-CN):
Center Nudging [Scellier et al., ] is used in equilibrium propagation to improve and stabilize performance compared to both positive and negative nudging, and it is obtained as an average of the gradients produced by the two methods. Here, we approximate this behavior by randomly alternating between epochs in which we train with either negative or positive nudging. In this way, the training model can benefit from both methods without any extra computational cost.
- 5.
Incremental PC (iPC), a simple and recently proposed modification where the weight parameters are updated alongside the latent variables at every time step [Salvatori et al., ].
- 6.
Standard Backpropagation with Cross-Entropy Loss (BP-CE) / Squared Error Loss (BP-SE):
the most popular way to do the credit assignment in the neural networks. The model is trained by computing the gradients of the loss function with the weights of the network using the chain rule.

#### Experiments.

The benchmark results of MLP are obtained with MNIST and Fashion-MNIST, the results of VGG-5 are obtained with CIFAR-10, CIFAR-100 and Tiny ImageNet, the results of VGG-7 are obtained with CIFAR-100 and Tiny ImageNet. The data is normalized as in Table 6.

**Table 6: Data normalization**
|  | Mean ($\mu$) | Std ($\sigma$) |
| --- | --- | --- |
| MNIST | 0.5 | 0.5 |
| Fashion-MNIST | 0.5 | 0.5 |
| CIFAR-10 | [0.4914, 0.4822, 0.4465] | [0.2023, 0.1994, 0.2010] |
| CIFAR-100 | [0.5071, 0.4867, 0.4408] | [0.2675, 0.2565, 0.2761] |
| Tiny ImageNet | [0.485, 0.456, 0.406] | [0.229, 0.224, 0.225] |

For data augmentation on the training sets of CIFAR-10, CIFAR-100, and Tiny ImageNet, we apply random horizontal flipping with a probability of 50%. Additionally, we employ random cropping with different settings for each dataset. For CIFAR-10 and CIFAR-100, images are randomly cropped to 32×32 resolution with a padding of 4 pixels on each side. In the case of Tiny ImageNet, random cropping is performed to obtain 56×56 resolution images without any padding. And on the testing set of Tiny ImageNet, we use center cropping to extract 56×56 resolution images, also without padding, since the original resolution of Tiny ImageNet is 64x64.

The model hyperparameters are determined using the search space shown in Table 7. The results presented in Table 1 were obtained using 5 seeds with the optimal hyperparameters, which are stored in the YAML files located in the "examples/s4_1_discriminative_mode/" subdirectories of the PCX library.

As for the optimizer and scheduler, we use mini-batch gradient descent (SGD) with momentum as the optimizer for the $h$, and we utilize AdamW with weight decay as the optimizer for the $\theta$. Additionally, we apply a warmup-cosine-annealing scheduler without restart for the learning rates of $\theta$. We also tried SGD with momentum for the weights $\theta$, but we did not perform complete hyperaparemter searches on all combinations of architectures and datasets as its performance was suboptimal compared to AdamW in all the cases tested.

**Table 7: Hyperparameters search configuration**
| Parameter | PC | iPC | BP |
| --- | --- | --- | --- |
| Epoch (MLP) | 25 |  |  |
| Epoch (VGG-5 and VGG-7) | 50 |  |  |
| Batch Size | 128 |  |  |
| Activation | [leaky relu, gelu, hard tanh] | [leaky relu, gelu, hard tanh, relu] |  |
| $\beta$ | [0.0, 1.0], 0.05^1 | - | - |
| $\beta_{ir}$ | [0.02, 0.0] | - | - |
| $lr_{h}$ | (1e-2, 5e-1)^2 | (1e-2, 1.0)^2 | - |
| $lr_{\theta}$ | (1e-5, 3e-4)^2 | (3e-5, 3e-4)^2 |  |
| $momentum_{h}$ | [0.0, 1.0], 0.05^1 | - |  |
| $weightdecay_{\theta}$ | (1e-5, 1e-2)^2 | (1e-5, 1e-1)^2 | (1e-5, 1e-2)^2 |
| T (MLP and VGG-5) | [4,5,6,7,8] | - |  |
| T (VGG-7) | [8,9,10,11,12] | - |  |

#### Results.

All the results presented in this study were obtained using forward initialization, a technique that initializes the model’s parameters by performing a forward pass on a zero tensor with the same shape as the input data. Besides, in our experiments, we limited the range of $T$ to ensure a fair comparison with BP in terms of training times. Higher $T$ correspond to a greater number of optimization rounds of $h$, which can lead to improved model performance but also increased computational costs and longer training durations. To maintain comparability with BP, we restricted our searching space of T that resulted in training times similar to those observed in BP-based training.

Momentum helps significantly. In Figure 9, we present the accuracy of the VGG-7 model trained on CIFAR-100 using different momentum values, both without nudging(Figure 9(a)) and with nudging(Figure 9(b)). It is evident from Figure 9 that selecting an appropriate momentum value can substantially improve model accuracy.
By comparing Figures 9(a) and 9(b), we can observe that different training algorithms have different optimal momentum values. The optimal momentum for training with nudging is generally higher than that for training without nudging. Furthermore, the optimal momentum for negative nudging is larger than that for positive nudging. These differences in optimal momentum values highlight the importance of carefully tuning the momentum hyperparameter based on the specific training algorithm and nudging method employed. For reference, the optimal model parameters and momentum values for various tasks and models can be found in the example/discriminative_experiments folder of the PCX library.

Figure: (a)
Refer to caption: x15.png

Activation function also plays a crucial role in improving model accuracy. For models using Cross-Entropy Loss, the “HardTanh” activation function is a better choice. In the case of models using Squared Error Loss without nudging, the “LeakyReLU” activation function tends to perform better. When using Positive Nudging, the optimal activation function varies depending on the model architecture. For Negative Nudging, the “GeLU” activation function is the most suitable choice.

Nudging improves performance. Fig. 11 illustrates the relationship between the learning rate of $h$ and accuracy with or without nudging. From the plot, we can observe that when nudging is not used (red dots), the model achieves better results at lower learning rates. However, when nudging is employed (purple and blue dots), regardless of whether it is positive nudging or negative nudging, the model can attain better accuracy at higher learning rates compared to the case without nudging. Additionally, Fig. 9(b) shows the relationship between momentum and accuracy. We can see that after applying nudging, the model can achieve better results at higher momentum values. We believe this is the reason why nudging can improve performance. The ability to use higher learning rates and momentum values without sacrificing accuracy is a significant advantage of nudging, as it can lead to faster convergence and improved generalization performance.

Figure: Figure 10: Different training runs for VGG-7 on TinyImageNet using CN. Most of the seeds diverge, resulting in poor average performance.
Refer to caption: x17.png

Training instability and model size.
Similarly to what we show in Section 5.1, we noticed that largest architectures present significantly more training instability. iPC is stable and produces optimal results only for the smallest, fully connected, architectures, while PC achieves its peak performance on VGG-5. In particular, we noticed that VGG-7 is able to reach similar highest accuracy (e.g., $\approx 43\%$ on Tiny ImageNet) but, on average, performs notably worse, as most seeds result in a diverged model (Fig. 10).

Figure: (a) VGG-7
Refer to caption: x18.png

## Appendix C Generative experiments

### C.1 Autoencoder

An Autoencoder is a network that learns how to compress a high-dimensional input into a much smaller dimensional space, called the bottleneck dimension or the hidden dimension, as accurately as possible. Thus, a backpropagation-based Autoencoder consists of two parts: an encoder, that compresses the input from the original high-dimensional space into the bottleneck dimension, and a decoder, that reconstructs the original input from the bottleneck dimension. A mean-squared error (MSE) between the original and the reconstructed input is used as a loss to train the Autoencoder network in an unsupervised manner.

Predictive Coding (PC) alleviates the need in the encoder part of an Autoencoder. Specifically, only the decoder part of an Autoencoder is used, with a PC layer acting as the bottleneck dimension and as an input to the decoder. Moreover, PC layers are inserted after each layer of the decoder.

Figure: Figure 12: Left: An Autoencoder implemented with backpropagation consists of both an encoder and a decoder. The encoder compresses the input data into the bottleneck dimension, and the decoder restores the original image. Right: An Autoencoder implemented with Predictive Coding. The state of the first PC layer is the bottleneck dimension. The state of the last PC layer is the original input, and the predicted state of the last PC layer is the predicted input. Inference steps update the bottleneck dimension to make it a good compressed representation.
Refer to caption: extracted/5702260/Figures/appendix/autoencoder-bp-vs-pc.png

**Table 8: Hyperparameters and search spaces for deconvolution-based autoencoders**
| Parameter | PC | iPC | BP |
| --- | --- | --- | --- |
| Number of layers | 3 | conv layers: 3 deconv layers: 3 |  |
| Internal state dimension | 4x4 |  |  |
| Internal state channels | 8 |  |  |
| Kernel size | [3, 4, 5, 7] |  |  |
| Activation function | [relu, leaky_relu, gelu, tanh, hard_tanh] |  |  |
| Batch size | 200 |  |  |
| Epochs | 30 |  |  |
| T | 20 | - |  |
| Optim $h$ | SGD+momentum | - |  |
| $lr_{h}$ | (1e-2, 5e-1)^2 | (1e-2, 1.0)^2 | - |
| $momentum_{h}$ | [0.0, 0.95] | - |  |
| Optim $\theta$ | AdamW |  |  |
| $lr_{\theta}$ | 3e-5, 1e-3^2 |  |  |
| $weightdecay_{\theta}$ | (1e-5, 1e-2)^2 | (1e-5, 1e-1)^2 | (1e-5, 1e-2)^2 |

**Table 9: Hyperparameters and search spaces for linear-based autoencoders**
| Parameter | PC | iPC | BP |
| --- | --- | --- | --- |
| Number of layers | 3 | encoder: 3 decoder: 3 |  |
| Internal state dimension | 64 |  |  |
| Activation function | [relu, leaky_relu, gelu, tanh, hard_tanh] |  |  |
| Batch size | 200 |  |  |
| Epochs | 30 |  |  |
| T | 20 | - |  |
| Optim $h$ | SGD+momentum | - |  |
| $lr_{h}$ | (1e-2, 5e-1)^2 | (1e-2, 1.0)^2 | - |
| $momentum_{h}$ | [0.0, 0.95] | - |  |
| Optim $\theta$ | AdamW |  |  |
| $lr_{\theta}$ | (3e-5, 1e-3)^2 |  |  |
| $weightdecay_{\theta}$ | (1e-5, 1e-2)^2 | (1e-5, 1e-1)^2 | (1e-5, 1e-2)^2 |

A PC-based Autoencoder works as follows:

- 1.
The energy function of the last PC layer is set to MSE upon its creation. In PCX, the squared error is the default energy function. The squared error is then summed across all dimensions in the input and averaged over the batch, that approximates the MSE up to a multiplication constant.
- 2.
The current state of the last PC Layer $L$, $h_{L}$, is fixed to the original input data, which means that $h_{L}$ is not changed during inference steps.
- 3.
Since the energy of the last layer $L$ now encodes the MSE loss between the predicted image $\mu_{L}$ and the original input stored as $h_{L}$, the inference steps will update the current states $h_{l}$ of all PC layers but the last one, including the one that represents the bottleneck dimension, to minimize this MSE loss.
- 4.
Once the inference steps are done, the state of the bottleneck dimension PC layer will converge to the compressed representation of the original input.

### C.2 MCPC

Model.
Monte Carlo predictive coding (MCPC) is a version of predictive coding that can be used for generative learning. MCPC differs from PC by its noisy neural dynamics. Unlike PC where the neural activity converges to a mode of the free-energy, the neural activity of MCPC performs noisy gradient descent which is used for Monte Carlo sampling. When an input is provided, the noisy neural activity samples the posterior distribution of the generative model given the sensory input. When no input is provided the neural activity samples the generative model encoded in the model parameters.
Specifically, the neural dynamics of MCPC leverage the following Langevin dynamics:

$$ $\Delta h_{l}=-\gamma\nabla_{h_{l}}\mathcal{F}_{h_{l}}(h,\theta)+\sqrt{2\gamma}N$ (5) $$

where N is a Gaussian random variable with variance $\sigma_{mcpc}^{2}$.
These neural dynamics can be extended to 2nd-order Langevin dynamics for faster sampling:

$$ $\displaystyle\Delta h_{l}=\gamma r_{l}$ (6) $\displaystyle\Delta r_{l}=\gamma\nabla_{h_{l}}\mathcal{F}(h,\theta)-\gamma(1-m )r_{l}+\sqrt{2(1-m)\gamma}N$ (7) $$

where $m$ is a momentum constant.

An MCPC model is trained following a Monte Carlo expectation maximisation scheme which iterates over the following two steps: (i) MCPC’s neural activity samples the model’s posterior distribution for the given data, and (ii) the model parameters are updated to increase the model log-likelihood under the samples of the posterior. In practice, we run MCPC inference for a limited number of steps after which we update the model parameters with a single sample of the posterior similarly to how model parameters are updated in variational auto encoders.

After training, samples of a trained model are generated by leaving all neurons unclamped and recording the activity of input neurons (the neurons clamped to data during training). The activity is recorded after a limited number of activity update steps. This process is repeated for each data sample.

MCPC’s implementation in PCX utilizes a noisy SGD optimizer for the state $h$. Compared to PC than uses an SGD or Adam optimizer, MCPC incorporates an optimizer that merges the addition of noise to the model’s gradients with an SGD optimizer.
The variance of the noise added to the gradients needs to be carefully crafted to scale appropriately with the learning rate and the momentum as shown in equations (5 - 7).

Experiments.
All the MCPC experiments use feedforward models with Squared Error (SE) loss. The SE loss of the state layer $h_{L}$ is also scaled by a variance parameter $\sigma^{2}_{h_{L}}$. This additional parameter is introduced to prevent the Gaussian layer $h_{L}$ from having a variance much larger than the variance of the data which would prevent learning. Moreover, for unconditional learning and generation, the layer $h_{0}$ is left unclamped during both training and generation. In contrast, for the conditional learning task on MNIST, the layer $h_{0}$ is clamped to labels during training and generation.

For the iris dataset, we train a model with layer dimensions [2 x 64 x 2], tanh activation function and default parameter values (state learning rate $\gamma$=0.01, state momentum = 0.9 , noise state variance $\sigma_{mcpc}^{2}=1$, parameter learning rate $lr_{\theta}$, parameter decay = 0.0001, Adam parameter optimizer, layer variance $\sigma^{2}_{h_{L}}$ = 0.01 and a batch size of 150). We use 500 state update steps during learning and 10000 for generation.

For the unconditional learning task on MNIST, we train models with layer dimensions [30 x 256 x 256 x 256 x 784]. The model hyperparameters for MCPC and VAE were determined using the hyperparameter search shown in table 10 to optimize the FID and the inception score separately. Refer to the code for exact optimal parameter values. We use 1000 state update steps during learning and 10000 for generation.

For the conditional learning task on MNIST, we train models with layer dimensions [2 x 256 x 256 x 256 x 784]. The labels used in this task, clamped to $h_{0}$, specify whether an image corresponds to an even or odd number. The model hyperparameters are determined using the search space shown in table 10. We use 1000 state update steps during learning and 10000 for generation.

Results.
Figure 13 shows samples generated by the trained models for hyperparameters that maximize the inception score.

**Table 10: Bayes hyperparameter search configuration for MCPC and VAE (where applicable) on MNIST.**
| Parameter | Value |
| --- | --- |
| activation | {ReLU, Silu, Tanh, Leaky-ReLU, Hard-Tanh} |
| $\gamma$ | log-uniform(0.0001, 0.05) |
| momentum | {0.0, 0.9} |
| $\sigma^{2}_{mcpc}$ | {1.0, 0.3, 0.01, 0.001} |
| $lr_{\theta}$ | log-uniform(0.0001, 0.1) |
| parameter decay | {0.0, 0.1, 0.01, 0.001, 0.0001} |
| $\sigma^{2}_{h_{L}}$ | log-uniform(0.03, 1.0) |
| batch size | {150, 300, 600, 900} |

Figure: Figure 13: Samples generated by trained models that optimize the inception score under the unconditional and conditional learning regimes.
Refer to caption: x20.png

### C.3 Associative memories

This section describes the experimental setup of associative memory tasks.

#### Model.

A generative PCN is first trained on $n$ images sampled from the Tiny ImageNet dataset until its parameters have converged. Then, a corrupted version of the training images is presented to the sensory layer of the model ($h_{L}$) and we run inference $\nabla h_{l}$ on all layers, including the sensory layer, until convergence. Note that in masked experiments, the intact top half of the images is kept fixed during inference. Intuitively, suppose the model has minimized its free energy with its sensory layer fixed at each of the $n$ training examples during training. In that case, it has formed attractors defined by these training examples and would thus tend to “refine" the corrupted images to fall back into the energy attractors.

#### Experiments.

Here, the benchmark results are obtained with Tiny ImageNet, corrupted with either Gaussian noise with 0.2 standard deviation, or a mask on the bottom half of the images (examples shown in Fig. 4). We vary the model size and number of training examples to memorize, to study the capacity of the models. Specifically, we use a generative PCN with architecture $[512,d,d,12288]$ where $d=[512,1024,2048]$ ($12288$ being the flattened Tiny ImageNet images) and varied $n=[50,100,250]$. We performed a hyperparameter search for each $d$ and $n$ on the parameter learning rate $lr_{\theta}\in\{1\times 10^{-4}+k\cdot 5\times 10^{-5}\mid k\in\mathbb{Z},0
\leq n\leq 18\}$, the state learning rate $\gamma\in\{0.1+k\cdot 0.05\mid k\in\mathbb{Z},0\leq n\leq 18\}$, training inference steps $T_{\textnormal{train}}\in[20,50,100]$ and recall inference steps $T_{\textnormal{recall}}\in[50000,100000]$. We fix the activation function of the model to Tanh, and the number of training epochs to $500$ and a batch size of $50$. The results in Table 3 are obtained with $5$ seeds with the searched optimal hyperparameters, which are stored in the hps.json file under the examples/s4_2_generative_mode/associative_memory folder in the PCX library.

## Appendix D Energy and Stability

This section describes the experimental setup of Section 5.1, provides replications on other datasets and ablations.

*The code can be found in*: examples/s5_analysis_and_metrics/energy_and_stability.

### D.1 Energy propagation

We test a grid of models on multiple datasets to examine the energy propagation in the models. We test on the FashionMNIST, Two Moons, and, Two Circles datasets. The Two Circles dataset is particularly interesting, as poor energy distribution intuitively results in a linear inductive bias (we primarily learn a one-layer network). This linear inductive bias harms the performance on Two Circles (linear model accuracy $\approx 50$%) more than FashionMNIST ($\approx 83$%) and Two Moons ($\approx 86$%).

#### Experimental Setup.

We train a grid of feedforward PCNs with 2 hidden layers. We train on three datasets: FahionMNIST (as reported in the main body) and additionally Two Moons and Two Circles. For all models, we train for 8 epochs with $T=8$ inference steps. States are optimized with SGD and forward initialization. The grid is formed over weight learning rate $lr_{\theta}\in\{1\times 10^{-5},1\times 10^{-4},\dots,1\}$, state learning rate $\gamma\in\{1\times 10^{-3},3\times 10^{-3},1\times 10^{-2},3\times 10^{-2},1
\times 10^{-1},3\times 10^{-1},1\}$, activation functions $f\in\{\text{LeakyReLU},\text{HardTanh}\}$ (the former is unbounded the latter is bounded), optimization with AdamW or SGD with momentum $m\in\{0.0,0.5,0.9,0.95\}$ and hidden widths of $\{512,1024,2048,4096\}$ for FashionMNIST and $\{128,256,512,1024\}$ for Two Moons and Two Cricles. We replicate all experiments on 3 seeds for FashionMNIST and 10 seeds for the other datasets.

#### Results.

Fig. 5(b) in the main paper shows the average energy across the last batch at the end of training for the best performing model on the grid. Fig. 5(c) compares SGD with momentum $0.9$ and AdamW. It is obtained for activation function “HardTanh” and a width of 1024. We replicate this figure for the other combinations of activation functions and widths below in Fig 14. We observe that across all conditions, small to medium state learning rates are generally preferred by SGD, while AdamW has a stronger preference to smaller state learning rates. Given the uneven distribution of energies across layers, Adam, in particular, may not scale to deeper architectures. We further, observe a larger variance in performance for Adam, especially for wider layers, which we discuss in paragraph “Training Instability“ in Sec. 5.1 and below.
Fig. 5(d) is based on all models trained with AdamW. Many models with high state learning rates diverge, we only plot models achieving accuracy $>0.5$.

Figure: Figure 14: Model accuracies for a range of combinations of activation functions and model widths. Adam perfers small learning rates and tends to be less stable than SGD. Obtained on FashionMNIST.
Refer to caption: x21.png

Figure: (a)
Refer to caption: x22.png

Figure: (a)
Refer to caption: x25.png

Below we present the results of experiments on the Two Moons and Two Circles datasets. Fig. 15(b), 15(a), and 15(c) replicate figures Fig. 5(b),  5(c), and 5(d) for Two Moons, and Fig. 16(b), 16(a), and 16(c) for Two Circles. Results are very similar to FashionMNIST: The energy is concentrated in the last layer, even after $T$ inference steps. However, in the example for Two Circles, we actually observe a training effect for earlier layers: While the energy increases first due to error propagation (still orders of magnitude below later layers), the energy is reduced afterwards. Energy ratios are consistenly indicating poor energy propagation for state learning rates $\gamma$, that perform well. As predicted the variance in results is significantly larger for Two Circles, especially for small state learning rates.

### D.2 Training Stability

We test a grid of PCNs to analyze the interaction between model width, state learning rates and weight optimizers.

#### Experimental Setup.

We train models on FashionMNIST (as reported above) and Two Moons. We train feedforward PCNs (2 hidden layers) with “LeakyReLU” activations over a grid of parameters. All models are trained over 8 epochs. The widths of the hidden layers are $\{32,64,\ldots,4096\}$. State variables are trained for $T=8$ steps with SGD and learning rates $\gamma\in\{1\times 10^{-5},3\times 10^{-5},\ldots,0.3\}$. The weights are updated through SGDor the Adam optimizer with a learning rate of 0.01 for FashionMNIST and 0.03 for Two Moons. Both optimizers uses 0.9 momentum for weights. We further train baseline BP models with the same hyperparameters. For FashionMNIST we replicate each run over 3 random initializations, for Two Moons over 10.

#### Results.

We replicate Fig. 6 (FashionMNIST) here for the Two Moons dataset, see Fig. 17. We observe effects for Two Moons that are analog to FashionMNIST as presented above: The stability of optimization strongly depends on the width of the hidden layers for Adam. This effect is not observed for SGD on either dataset. This further supports the our conclusion in Sec. 5.1: While Adam is the better optimizer, this interaction effect (width $\times$ $\gamma$) can hinder the scaling of PCNs with Adam. Optimization methods for PCNs require further attention from the research community.

Figure: Figure 17: The instability of optimization with Adam given architectural choices can be observed for Two Moons.
Refer to caption: x28.png

#### Ablation.

We further provide an ablation on FashionMNIST. In the experiments above, the hidden layer width is altered, introducing changes in the *absolute size* of the hidden layers (i.e. number of neurons), but also changing the *relative size* of the hidden layers in the network, as input and output layers remain the same size across all experiments. Hence, we provide another experiment on FashionMNIST, where we increase the image size and augment the label vector with 0 0s, such that the width of all layers is equal. All other experimental variables remain as described above. The results are shown in Fig. 18 and follow the trend observed in Fig. 6 and 17: We find that there exists an interaction between the optimization and the width of the network as described above. Hence, accounting for relative changes in layer width does *not* sufficiently explain the problem and we conclude that the *absolute* size of the layers plays a role in the stability of optimization with AdamW.

### D.3 Integrating Skip Connections into VGG19 to Enhance PC Networks for CIFAR10 Classification

#### Skip connections.

We investigate the integration of skip connections into the VGG19 architecture to enhance its performance on the CIFAR10 image classification task, showing a significant increase in maximum achieved test accuracy from 25.32% to 73.95%. The vanishing gradient problem, a notable challenge in deep Predictive Coding (PC) models, becomes pronounced with increased network depth, hindering error transmission to earlier layers and impacting learning efficacy. To address this, we introduce skip connections that allow gradients to bypass multiple layers, enhancing gradient flow and overall learning performance.

**Table 11: Hyperparameter configuration and best accuracy for VGG19 with and without skip connections on CIFAR10**
| Parameter | Range | Best Value |
| --- | --- | --- |
|  | With Skip Connections |  |
| Epochs | 30 | 30 |
| Batch size | 128 | 128 |
| Activation functions | {GELU, Leaky ReLU} | Leaky ReLU |
| Optimizer for network parameters - Learning rate | {5e-2, 1e-1, 5e-1} | 0.5 |
| Optimizer for network parameters - Momentum | {0.0, 0.5, 0.9, 0.99} | 0.5 |
| Optimizer for weight parameters - Learning rate | 1e-4 | 1e-4 |
| Optimizer for weight parameters - Weight decay | {5e-4, 1e-4, 5e-5} | 5e-4 |
| Number of inference steps (T) | {24, 36} | 24 |
| Best Accuracy |  | 73.95% |
|  | Without Skip Connections |  |
| Epochs | 30 | 30 |
| Batch size | 128 | 128 |
| Activation functions | {GELU, Leaky ReLU} | GELU (default) |
| Optimizer for network parameters - Learning rate | {5e-2, 1e-1, 5e-1} | 0.1 |
| Optimizer for network parameters - Momentum | {0.0, 0.5, 0.9, 0.99} | 0.99 |
| Optimizer for weight parameters - Learning rate | 1e-4 | 1e-4 |
| Optimizer for weight parameters - Weight decay | {5e-4, 1e-4, 5e-5} | 1e-4 |
| Number of inference steps (T) | {24, 36} | 24 |
| Best Accuracy |  | 25.32% |

Figure: Figure 19: Performance comparison of VGG19 with and without skip connections on the CIFAR-10 dataset over 30 epochs. The plot shows the mean test accuracy along with the shaded area representing the variability across three different seeds.
Refer to caption: x30.png

#### Results

Our modified VGG19 model includes a skip connection from an early layer within the feature extraction stage, with the output flattened and adjusted using a linear layer before being reintegrated during the classification stage. The model underwent rigorous training and evaluation on the CIFAR10 dataset, employing standard preprocessing techniques like normalization and data augmentation (horizontal flips and rotations). Detailed hyperparameter tuning revealed optimal configurations for both models, with and without skip connections, exploring various optimizers, learning rates, momentum values, and weight decay settings, significantly enhancing the model performance with skip connections as summarized in Table 11.Figure 19 shows the test accuracy progression over 30 epochs for the VGG19 model with and without skip connections on the CIFAR10 dataset, using three different seed values and identical hyperparameters for both simulations. All experiments and scripts used in these experiments can be found in the examples/s5_analysis_and_metrics/energy_and_stability/skip_connections folder of the PCX library.

## Appendix E Properties of predictive coding networks

This section describes the experimental setup of Section 5.2 and displays the utility of using the free energy of a PCN classifier to differentiate between in-distribution (ID) and out-of-distribution (OOD) data Liu et al. []. We show how one can compute the negative log-likelihood of various datasets Grathwohl et al. [] under the PCN. We further provide analyses on the relationship between maximum softmax values and energy values before convergence and after convergence at the state optimum. We compare results across multiple datasets to corroborate our results as well as to show how PCNs can be used for OOD detection out of the box based on a single trained PCN classifier for which we study the receiver operating characteristic (ROC) curve based on different percentiles of the softmax and energy scores.

### E.1 Free energy and out-of-distribution data.

#### Experimental Setup.

We train a PCN classifier on MNIST using a feedforward PCNs with 3 hidden layers each of size $H=512$ with “GELU” activation and cross entropy loss in the output layer. We train the model until test error convergence using early stopping at epoch 75. During training the state variables are optimized for $T=10$ steps with SGD and state learning rate $\gamma=0.01$ without momentum. The weights are optimized using the SGD optimizer with a momentum of $m_{\theta}=0.9$ and the weight learning rate is chosen as $lr_{\theta}=0.01$. During test-time inference, we optimize the state variables until convergence for $T=100$. To understand the confidence of a PCN’s predictions, we compare the distribution of energy for ID and OOD samples against the distribution of the softmax scores that the classifier generates. We compute negative log-likelihoods for ID and OOD samples under the PCN classifier via:

$$ $\mathcal{F}=-\ln{p(x,y;\theta)}\implies p(x,y;\theta)=e^{-\mathcal{F}},$ (8) $$

We conduct the experiments on MNIST as the in-distribution (ID) dataset and we compare it against various out-of-distribution datasets such as notMNIST, KMNIST, EMNIST (letters) as well as FashionMNIST.

#### Results.

In the following we briefly interpret the additional results on the basis of experiments supported by various figures

Figure: Figure 20: Energy distributions before and after state optimization.
Refer to caption: x31.png

In Fig. 20 we see how the energy is distributed at test-time before and after state optimization. We can see, that all OOD datasets have significantly larger initial energies as well as final energies compared to the ID dataset (MNIST).

Figure: Figure 21: Energy histograms against ID data before and after state optimization.
Refer to caption: x32.png

In Fig. 21 we then show how each energy distribution for the OOD dataset compares against the energy of the in-distribution dataset by overlaying the histograms of the energies before and after state optimization. We can see that by plotting the histograms, a pattern emerges, namely, that a majority of the OOD data samples do not overlap with ID data samples, which supports the idea that energy can be used for OOD detection.

Figure: Figure 22: Softmax histograms overlapped with ID dataset.
Refer to caption: x33.png

Next in Fig. 22 we show how this pattern might look like when comparing the softmax scores of ID against OOD datasets. One can see, that the softmax scores are less informative for determining if samples are OOD as can be seen by the bigger overlap in the range of softmax values that ID and OOD samples have in common.

Figure: Figure 23: Non-linear relationship between energy and softmax scores.
Refer to caption: x34.png

In Fig. 23 we further study the relationship between softmax scores and energy values before and after state convergence. The plot shows that while the energy and softmax scores are strongly correlated before inference, a non-linear relationship is evident after convergence, especially for smaller values where the model is more uncertain. This indicates, that softmax scores and energy values do not fully agree on which samples we should have less confidence in.

Figure: Figure 24: Energy and NLL for various OOD datasets before and after inference.
Refer to caption: x35.png

In Fig. 24 we show how the energy distributions for all datasets look like before and after inference. Each box plot represents a different scenario and a different dataset. In addition, we compute the NLL of each dataset and display it as part of the box plot labels. We observe that across all OOD datasets, the initial and final energy values are significantly higher than the MNIST (ID) dataset. Furthermore, we can see that the variance of the energy scores is smaller for the in-distribution data as can be seen by the fact, that there are no outlier samples for MNIST beyond the whiskers of the box plot. Finally, the NLL values for each scenario confirm this observation, with the likelihood of the MNIST data being significantly higher than that of the OOD distributions.

Figure: Figure 25: Performing OOD detection with PCN energy and classifier softmax scores.
Refer to caption: x36.png

Finally, in Fig. 25 we show how the PCN can be used to classify samples as belonging to the ID or some OOD data. We use the PCN classifier’s energy to perform OOD detection and we show that the ROC curves for energy-based detection are superior to ROC curves created via softmax scores. This observation becomes even clearer, when looking at the most challenging samples by picking the 25% percentile of the scores and energies, in effect the samples, that the PCN model is least confident about as reflected by small energy or softmax values.

## Appendix F Computational Resources

Fig. 8 was obtained by taking a small feedforward PCN made by 2 layers of 64 neurons each and training it on batches of 32 elements (generated as random noise so to avoid any overhead due to loading training data to the GPU) for $T=8$ steps. Then, each parameter was scaled independently to measure its effect on the total training time. Each model obtained this way was trained for 5 epochs and the mean time was reported. In all our timing measurements, we skip the first epoch to avoid including the JIT compilation time. Results were obtained on a GTX TITAN X, showing that parallelization is potentially achievable also on consumer GPUs.

## Peer Reviews

confidence: 4
contribution: 4
initial_rating: 8
presentation: 3
questions: 1. Why only JAX? Why not Pytorch? It seems like this goes against the spirit of the inclusiveness of your work. 
2. On Scalability Solutions: While the paper does an excellent job of identifying scalability issues in PCNs, do the authors have any preliminary ideas or suggestions for future work that could address these challenges?

3. Efficiency vs. Backpropagation: Given that the proposed PCX library is still slower than traditional backpropagation for complex models, are there any plans or ongoing work to bridge this efficiency gap?
3. Spelling, grammar and clarity points:

Here are the sections that could benefit from improved clarity, fluency, or grammar:

1. Line 110-111:
"different directions of research that either explored interesting properties of PC networks, such as their robustness and flexibility (Song et al., 2024; Alonso et al., 2022) , or proposed variations to"
- Extra space before comma after parenthesis
- Somewhat awkward phrasing

2. Lines 180-182: 
"In practice, we do not train on a single pair (x,y) but on a dataset split in mini-batches that are subsequently used to train the model parameters."
- Should be "split into mini-batches"

3. Line 229: 
"While interesting and important for the progress of the field,"
- Awkward start to sentence, could be rephrased for better flow

4. Line 270:
"We test the performance of PCNs on image generation tasks."
- Somewhat abrupt transition, could use better connection to previous section

5. Line 394:
"We observed a further link between the weight optimizer and the structure of a PCN that might hinder the scalability of PC"
- Run-on sentence that could be clearer if split or restructured

6. Line 443:
"PCX extensively relies on Just-In-Time compilation."
- Inconsistent capitalization of "Just-In-Time" compared to other instances of "just-in-time"

7. Line 469:
"In fact, in its current state, JIT is unable to parallelize the executions of the layers;"
- "executions" should be singular "execution"

8. Line 476:
"we have lied new foundations"
- Should be "laid new foundations"

9. Line 485:
"When this condition is released,"
- Should be "When this condition is relaxed," or "When this limitation is removed,"

10. Throughout the paper:
- Inconsistent spacing around mathematical symbols and equations
- Some inconsistent use of commas in complex sentences
- Occasional mixing of American and British English spellings

Figure 1: 

You want a figure to stand on its own without any context from the paper ideally. What is F here? Can we improve the clarity and informativeness of the caption and the figure itself?

These issues are relatively minor and don't significantly impact the paper's readability or scientific contribution, but addressing them would improve the overall polish of the manuscript.
soundness: 4
strengths: - Addresses Core Issues in PC Research: The paper tackles fundamental shortcomings in the subfield, such as inconsistent comparisons, overly simple benchmarks, and a lack of rigorous, standardized practices.

- Unified Benchmarking Framework: Introduces a comprehensive and standardized benchmark for PCN models, incorporating major datasets like CIFAR-100, Tiny ImageNet, FashionMNIST, and MNIST to ensure consistent and meaningful comparisons.

- Implementation of Open-Source Tool (PCX): Offers a JAX-based library that enhances research efficiency and accessibility, boasting up to 50x speed improvements via JIT compilation and seamless compatibility with modern ML tools.

- Thorough Experimental Evaluation: Conducts extensive tests on key PCN architectures and variants, backed by large-scale hyperparameter studies, establishing a solid baseline for future advancements.

- Deep Analysis of Training Mechanics: Provides valuable insights into energy propagation, initialization strategies, and stability across optimizers and architectures, which can guide targeted future research.

- Community-Centric Contribution: The open-source nature, complete with detailed documentation and tutorials, lowers the barrier to entry and encourages collaboration, setting a strong foundation for collective progress.

- Motivation for Future Research: Clearly identifies current limitations, such as scalability and energy flow issues, laying the groundwork for targeted solutions and deeper exploration.

- Competitive Performance on Benchmarks: Shows that PCNs can achieve performance on par with BP-trained models for small to medium-sized architectures, establishing a base for scaling up and further innovation.

Insightful 'weaknesses' highlighted by the paper that should not be perceived as a weakness of the paper itself imho.

Identification of Scalability Challenges: The paper highlights scalability issues inherent to PCNs but does not propose new solutions. While this reflects well on the depth of the analysis, it may be seen as a limitation in terms of contribution.
Performance Bottlenecks in PCNs: The analysis shows that while the proposed library offers performance improvements over prior implementations, it remains slower than backpropagation for more complex models. This insight underscores the current challenges facing PCNs but could be perceived as a weakness in terms of efficiency.
Sequential Processing Constraints: The authors identify that sequential layer updates are a significant bottleneck, impacting parallelism and efficiency. This is a critical finding for the field, though it might be seen as a shortcoming in practical use.
Theoretical Insight Depth: The paper provides an analysis of training mechanics, but it could go further in offering deep theoretical explanations for issues like energy propagation and training stability. This gap invites further exploration but could be seen as a weakness by readers looking for comprehensive theoretical coverage.
Implementation Constraints and Parallelization: The paper highlights that training time scales linearly with the number of layers due to current parallelization limitations. This insight is crucial for understanding the limitations of PCNs but might be perceived as a practical constraint in the work presented.
summary: The authors undertake the very hard task of trying to streamline, standardise and robustify the scientific process in the subfield of Predictive Coding networks. They identify that the subfield suffers from a variety of issues such as poor comparisons between models, too simple baselines or benchmarks, and overall a lack of rigour in how potential proposed techniques are studied, compared, contrasted and analysed. This leads to the issue of making important issues hard to identify, and even harder to get a good critical mass of people working on them for resolution. 

The authors try to improve on this issues by offering: 

1.  A unified benchmark for predictive coding network models that is more expansive than others before it, utilising CIFAR100, Tiny ImageNet, FashionMNIST and MNIST. 
2. A library that allows quick and easy application of the benchmark on a target approach. 
3. A thorough evaluation of a number of key architectures and PC approaches on their benchmark, that can serve both a source of insights for the current state of the field and a good starting baseline for future approaches. 

Overall the proposal of the authors is an important and valuable one, and while the writing is OK, it is at times sloppy, with fluency and spelling mistakes, as well as inconsistent use of various terms. 

However overall a very good paper.
weaknesses: Weaknesses of the Paper:
- Writing and Terminology: The paper has inconsistencies in terminology and occasional grammar and spelling issues, along with unclear phrasing that affects readability and comprehension.
- Comparative Analysis Limitations: The paper could be strengthened by including more comparisons with other biologically inspired approaches and better situating its findings within the broader machine learning context.
- Practical Implications: The paper lacks a thorough exploration of how its findings translate to practical, real-world applications beyond academic benchmarks.
- Diversity of domains/tasks/modalities while a lot more than anything other out there for PC research, is still very limited compared to much more general benchmarks out there for general deep learning evaluation. Adding more domains, tasks and modalities would be very useful to having rigorous and insightful works done by people using it. 
- Architectures applied do not include any resnet, densenet or transformer architectures -- yes depth is an issue I understand, but using a shallow variant of any of these would be very useful.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. Can we test some models beyond VGG and MLP as mentioned in weakness? 
2. For the results in table 1, the result of VGG-5 is better VGG-7, does that mean PC works better when the number of layers is small. I think we need more data points here.
3. Also for table 1, can you include the same datasets for both VGG-5 and VGG-7, currently, the missing data points (CIFAR-10) makes it both hard to read the table and hard to make more in-depth comparison between these two cases.
4. I think it's better to provide the URLs at the beginning of the paper instead of at the end.

I will consider raising the scores if my questions can be addressed.
soundness: 3
strengths: 1. The benchmark results included in this work is very extensive which includes most of the popular baselines.
2. The design and usage of the contributed library are well documented with provided repos (I didn't run the code in this repo).
3. The analysis provided in the paper is well explained such as showing different results from different methods which makes the paper easy to read
4. With the contributed library, the paper is able to scale up the field to larger datasets which greatly improved the applicability of the methods in the field.
5. The work also shared some insights for sub-optimal cases of PCNs which can further guide the development of new methods.
summary: The paper contributes a comprehensive library/framework that unifies the developing and testing for predictive coding networks (PCNs). Traditionally, the field suffers from lacking of unified metrics and open-source implementations, the contributed library well addresses the problem. Furthermore, with the contributed library, the paper is able to scale existing works onto larger datasets and achieves new SOTA on all the datasets. The paper also includes the popular algorithms in the field and benchmarked them extensively.
weaknesses: 1. Some figures can be improved, e.g. for figure 2, the method name can be draw above the images for better representation.
2. The main models used in the paper is MLP and VGG, since the library has been developed, we should probably test a few more model types such as resnet (shallow is OK if deep ones cannot be test) or vit. This could probably real more insights on the overall states of PCN algorithms.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 2
questions: Can you motivate the use of VGG-5 and VGG-7 models in the benchmarks? In terms of model sizes, they are roughly 3.8M and 4.5M, respectively, so they’re quite close to each other.
soundness: 2
strengths: ### Originality
The main contribution of the paper is the PCX library which can be used to implement predictive coding networks, and the proposed benchmarks, which can be used to evaluate the scalability of new predictive coding networks.

### Quality 
The PCX library, owing to Just-In-Time compilation, is much faster than the baseline library. The reported performance (training times) of the PCX library shines against the baseline by Song et al.

### Clartity
The tutorials are of good quality. The Jupyter Notebook tutorials released with the paper are easy to set up and follow.

### Significance
The paper is well motivated.
1. Providing the community with the tools and benchmarks to compare their methods against others’ is a meaningful endeavour.
2. The benchmarks are designed to address scalability of PCNs.
3. The benchmarks are not limited to only discriminative tasks like classification, but also include generative tasks such as autoencoding and sampling from a learned joint distribution.

The following insights into PCNs were very interesting:
1. The AdamW optimizer performs poorly as the hidden layer size increases. This was also demonstrated well in the tutorial notebook.
2. Skip connections have a dramatic effect on the performance of PCNs.

The evaluation is extensive. A fair number of existing PCN methods were evaluated (standard PC, positive/negative/centered nudging PC, incremental PC, Monte Carlo PC).
summary: This paper addresses the trend in predictive coding networks (PCNs) literature of not comparing results against each other. To this end, the authors present a library that facilitates implementing and testing new PCNs, and they also propose several new benchmarks for evaluation. These benchmarks include larger networks and more complex tasks than those in previous works, and span discriminative as well as generative tasks. The paper evaluates several existing methods, including standard backpropagation, on these newly proposed benchmarks and analyzes the results.
weaknesses: I gave a score of 3 for Contribution but only 1 for Soundness. This is because I believe the PCX library is a valuable contribution to the community; however, I have reservations about some of the conclusions and claims made in the paper.

Consider the Energy propagation subsection of Section 5.1 (Energy and Stability). L377 says, “Note that the decay in performance as function of increasing γ is stronger for Adam despite being the overall better optimizer in our experiments. This suggests that possible directions for future improvements should aim at reducing the energy imbalance between layers.”. It is unclear to me how this conclusion that energy imbalance must be reduced was arrived at. I understand energy imbalance to be different from energy ratios in that the imbalance is minimized when the ratio is 1, i.e., both the layers have the same energy.
In Figure 5 (right), we see that the learning rate which minimizes the energy imbalances (i.e, make the ratio as close to 1 as possible) is 0.3. In Figure 5 (center right), we see that the performance of AdamW optimizer for a learning rate of 0.3 is much worse than 1e-3, 3e-3, and 1e-2, all of which have much higher energy imbalances. The same pattern can be observed in Figures 13 and 14. How do we conclude from this behaviour of AdamW that the energy imbalance must be minimzed? This also affects the claim in the abstract that the paper “clearly highlights what the current limitations of PCNs are”.

Next, in Section 4.1 (Discriminative Mode), L243 says that "The results show that the best performing algorithms, at least on the most complex tasks, are the ones where the target is nudged towards the real label, that are PN and NN.". The results for Discriminative mode are reported in Table 1. Assuming, based on L067, that complex tasks refer to CIFAR100 and Tiny ImageNet, I am unable to find data in the table that substantiates this claim. For VGG-5 models, centered nudging (CN) dominates both positive nudging (PN) and negative nudging (NN) for all the tasks except Tiny ImageNet (Top-1).
For VGG-7 models, all PC methods are beaten by either of the two variants of Back Propagation (BP-CE and BP-SE). Even considering only the PC methods, PN is consistently beaten by PC-CE, as is NN, with the exception of CIFAR-100 (Top-1).

Finally, in the Discussion subsection of Section 4.2 (Generative Mode), line 339 claims that “Thus, compared to artificial neural networks, PCNs are more flexible and require only half the parameters to achieve similar performance”. The fact that PCNs require only half as many parameters as an autoencoder neural network is evident because they can infer the hidden states of layers that have not been fixed (either to the input image for compression or to the compressed representation for generation). However, the paper demonstrates that this holds true only for small models, as the generative experiments were conducted exclusively on relatively small models (3 layers for autoencoding and generation via sampling, and 2 hidden layers for associative memory experiments). Conducting at least one experiment on a larger model, such as a VGG-7-scale model, could have strengthened this claim.

I am open to reconsidering my scores if these concerns are addressed.

The following is feedback and minor nitpicks about the presentation of the paper. Addressing these will improve the quality of the paper, but will not affect the scores.

1. The abbreviation “BP” is used in the paper without introducing it. The first reference to it is in the Appendix.
2. Figure 1 (b). Shouldn’t it be “hL is fixed to y”?
3. Table 3 can be made more intuitive to read by mentioning in the caption that it is a comparison of the two modes of associative memory tasks - Noise and Mask.
4. H and theta should be subscripts in L141
5. Citations are missing for several works (which, admittedly, a lot of us take for granted). Examples are CIFAR-10, MNIST, Fashion-MNIST, Tiny ImageNet, SGD, Adam, AdamW, VGG models.
6. Citations for each of the methods in the Algorithms subsection of Section 4 would make it easier for the readers to refer to the papers, even though they are all cited before.
7. Figure 5 would be better if it were indexed as Fig 5a to Fig 5d, instead of Fig 5 left, center-left, center-right and right.
8. The citation style in L239 is wrong.
9. L347 “How regularly the energy flows into the model”. Improve phrasing.
10. All references to Figure 5 right and center-right are inverted. See, e.g., L373 and L377.
11. Grammatical errors in a few places, such as L267, L355, L467

Meta-review: The paper makes several contributions to predictive coding networks in machine learning in the form of a library to enable and facilitate comparisons, an attempt at unification of benchmarks, as well as rigorous empirical evaluation. These latter points address several identified pain points and open aspects in existing literature. As such, all reviewers agree that the paper is valuable, necessary, and will likely have impact in the community. The AC agrees with these points and believes the paper presents a timely contribution to advance the field. The raised weaknesses were mostly addressed in the discussion period. Overall, the AC recommends to accept the paper and agrees with the reviewers that the paper’s contributions are extensive and the library helpful to the community.

Decision: Accept (Spotlight)