\chapter{Temporal Predictive Coding}\label{chap:time}

\section{Introduction}
This study is concerned with extending the theory of predictive coding to the problem of processing dynamically changing sequences of sensory inputs arriving over time. Predictive coding, which originated from an algorithm for compression in information theory \citep{spratling2017review}, was initially developed and applied to the analysis of the brain by \citet{srinivasan1982predictive} and \citet{mumford1992computational} and then formalized into a general computational model of the cortex by  \citet{rao1999predictive}. The core hypothesis behind predictive coding is that the brain computes predictions of its observed input, and compares these predictions to the actually received input. The difference between the two is called the prediction error and quantifies how incorrect the brain's prediction was. Predictive coding proposes that the brain then adjusts its neural activities and synaptic strengths to minimize prediction errors which ultimately results in more accurate predictions \citep{friston2005theory,clark2015surfing}. Thus, solely by minimizing prediction errors, the brain is forced to learn a general \emph{world model} which can generate accurate predictions of its incoming sensory input \citep{friston2003learning,friston2005theory}. Moreover, these prediction errors are computed \emph{locally} between the local input to the neuron and the predictions it receives. This means that learning in predictive coding model requires only local information and can be accomplished in most cases with purely Hebbian synaptic plasticity
\citep{buckley2017free,bogacz2017tutorial,millidge2020relaxing}.

As has been introduced in Chapter \ref{chap:intro}, predictive coding has become an influential theoretical model for understanding cortical functions, and recent progresses have also established its capability in machine learning tasks and relationship to backpropagation. However, most of these works are concerned with inputs that are independently and identically distributed (i.i.d.) samples from some datasets and are presented in batches to the predictive coding model in random order. The visual input to the brain, on the other hand, is a continually changing sensory stream conveying a sequence of individual images with much correlation and rich structure embedded in the timing of the sequence elements. Therefore, to better describe the information processing in the brain, predictive coding models must take into account one more crucial element: \emph{time}.

% In their original study, \citet{rao1999predictive} trained predictive coding networks to generate static images of natural scenes, and demonstrated that the network learned receptive fields with similarities to those found in V1, as well as reproduced extra-classical receptive fields and end-stopping. Since then it has been demonstrated that predictive coding networks can explain many intriguing phenomena such as repetition suppression \citep{auksztulewicz2016repetition}, bistable perception \citep{hohwy2008predictive,weilnhammer2017predictive}, illusory motions \citep{watanabe2018illusory}, retinal stabilization \citep{millidge2019fixational} and even potentially attentional modulation of perception \citep{feldman2010attention,kanai2015cerebral}. Moreover, there has been much work matching the underlying neurophysiology of cortical microcircuits to the fundamental computations required by the predictive coding algorithm \citep{keller2018predictive,bastos2012canonical,millidge2021predictive}, thus providing a potential low-level basis for the implementation of predictive coding in neural circuitry.

% Alongside the aforementioned works that have successfully reproduced many neurophysiological phenomena,
% %with predictive coding, 
% recent progresses in predictive coding have been concerned with machine learning tasks such as the classification or generation of static images \citep{song2020can,millidge2020predictive, ororbia2022neural, sun2020predictive}, and multiple lines of research have investigated the relationship between predictive coding and backpropagation, the driving force behind modern machine learning systems \citep{whittington2017approximation,song2020can}. 

There are several established algorithms in statistics and machine learning for sequence processing over time, but they typically require very complex computations that would be difficult for biological circuits to perform. Nevertheless, it is useful to consider them as a reference against which the performance of biologically plausible models can be assessed. When there is a linear relationship between the current and future states (and noise is assumed to be Gaussian), the optimal temporal predictions can be achieved by the Kalman filter \citep{kalman1960new}. For more complex nonlinear problems, one can employ recurrent neural networks \citep{jordan1990attractor}, which contain recurrent\footnote{Unlike recurrent PCNs in the previous chapter, the ``recurrent'' connections discussed in this chapter will involve time.} connections that maintain and update an internal hidden state over time. While these recurrent networks, and more advanced successors such as the Long-Short-Term-Memory (LSTM) \citep{hochreiter1997long} can be very expressive and powerful, they are typically trained with backpropagation through time algorithm (BPTT) \citep{williams1989learning,williams1992simple,werbos1988generalization} which requires storing a history of all computations through a sequence and then ``unrolling'' it sequentially backwards through time to make updates. This algorithm is biologically implausible, because the brain can only receive inputs in a sequential stream, and must be able to process them online, i.e. as the inputs are received, and seems unlikely to be able to unfold a precise sequence of computations in reverse. In this work, we focus our comparison to these statistical and machine learning models on the Kalman filter, an online algorithm that processes data sequentially, as biological systems appear to do.

Within the field of predictive coding, there are a few tracks of research that have considered incorporating the temporal dimension. Earlier works \citep{rao1997correlates, rao1997dynamic} employed Kalman filtering to model the visual processing of dynamical sequences.
%which has similar computational steps to that of predictive coding. 
Instead of using fixed transformation matrices which are assumed to be known as in Kalman filtering, these works introduced learning rules to model synaptic plasticity in neural circuits performing filtering. \citet{friston2008hierarchical} have proposed the notion of extending predictive coding to use \emph{generalized coordinates} \citep{friston2008DEM,friston2010generalised} which model a dynamical state by including a set of temporal derivatives in the state vector and making predictions of these derivatives along with the current state. A few recent studies have also developed predictive coding models that perform well in various tasks involving temporal dependencies, such as those commonly examined in machine learning \citep{kutschireiter2017nonlinear, ororbia2020continual, jiang2022dynamic}. 
However, the works mentioned above have departed from the simple and flexible architecture of classical predictive coding for static inputs \citep{rao1999predictive} in order to take into account the temporal dimension. For example, the pioneering work that introduced Kalman filtering into predictive coding \citep{rao1997correlates,rao1997dynamic} has not described how the computations of Kalman filtering could be implemented in biological networks. Mapping of models employing generalized coordinates \citep{friston2008hierarchical} on a neural circuit would require explicit hard-coding of the expected temporal dependencies between different dynamical orders. Other models include specialised network features to aid the temporal processing, such as hyper-networks \citep{jiang2022dynamic}, multiple sub-networks \citep{kutschireiter2017nonlinear} or complex connectivity and neuron types \citep{ororbia2020continual}.

In this study, we propose a simple predictive coding network that also incorporates the temporal dimension, which we call \emph{temporal predictive coding} (tPC). This model generates predictions not only about the current inputs but also about its own \emph{future} neural responses, which is achieved by recurrent connections between neurons to transmit the prediction of one time-step to the next. 
The study makes four main contributions: First, we propose a predictive coding model that addresses the problem of temporal prediction, while inheriting from static predictive coding \citep{rao1999predictive} the simple and biologically plausible neural network implementations employing only local connectivity and Hebbian plasticity rules. Second, when the model is linear, we show that our model is a close approximation of the Kalman filtering model analytically, and has empirical performance comparable to Kalman filtering in benchmark filtering tasks while being computationally cheaper and requiring only biologically plausible operations. Also, unlike the Kalman filter, our model can learn the parameters of the system online.
%This also connects our model to earlier linear models that directly used Kalman filtering to model the dynamical processing of sensory inputs in the brain \citep{rao1997correlates,rao1997dynamic}. 
{Third, when trained with natural moving stimuli, we find that the model develops localized, Gabor-like receptive fields similar to those observed by Rao and Ballard \citep{rao1999predictive}, and more importantly, the receptive fields are also motion-sensitive, a property unique to neurons responding to dynamic stimuli \citep{baker1990spatial}.}
Overall, our model provides a possible computational mechanism underlying the cortical processing of dynamic inputs based on predictive coding, suggesting that the brain may learn to represent both static and continuous sensory observations using a single computational framework.

\section{Models}


The structure of the exposition in this section is as follows. Firstly, we present the underlying graphical structure of our proposed generative tPC network (tPCN) i.e., a Hidden Markov Model (HMM). Next, we show that, with Gaussian assumptions on the HMM and certain parametric assumptions on the nonlinear generative processes, we can derive an objective function that is similar to the original predictive coding network \citep{rao1999predictive} but taking into account the temporal dimension. We then demonstrate that neural dynamics and plasticity rules can be derived via the minimization of the objective function by gradient descent, and a corresponding training algorithm is presented. Finally, we show that the proposed computations and algorithms afford biologically plausible neural implementations in several different cases and discuss how they may be mapped to neural circuitry in the cortex. Since this study contains much mathematical notation, for ease of reading, we have collected it in Table \ref{notation_table} for quick reference \footnote{Note that in this Chapter, we use $\bfx$ to denote hidden or latent state and $\bfy$ to denote observations. This is to echo the notational tradition in control theory as we will compare tPC with Kalman filtering. This notation is different from elsewhere in this thesis.}.

\begin{table}[!t]
\begin{center}
\caption{
{\bf Table of mathematical notation used in the study.} {Vector and matrix variables are defined with their dimensions. Otherwise, the variables are scalars.}}
\begin{tabular}{|c c|} 
 \hline
 Notation & Meaning  \\ [0.5ex] 
 \hline\hline
 $\mathbf{y}$ & Observation {($\mathbb{R}^{d_y}$)}\\ 
 \hline
 $\mathbf{x}$ & Latent state {($\mathbb{R}^{d_x}$)}\\
 \hline
 $\mathbf{u}$ & Control input {($\mathbb{R}^{d_u}$)}\\
 \hline
 $k$ & Observation ``frame'' \\
 \hline
 $t$ & Inference time \\
  \hline
 $\bfW$ & Dynamics matrix {($\mathbb{R}^{d_x \times d_x}$)}\\ 
 \hline
 $\bfB$ & Control matrix {($\mathbb{R}^{d_x \times d_u}$)}\\
 \hline
 $\bfF$ & Observation matrix {($\mathbb{R}^{d_y \times d_x}$)}\\
 \hline
 $\bm{\Sigma}_x$ & Process noise {($\mathbb{R}^{d_x}$)}\\
 \hline
 $\bm{\Sigma}_y$ & Observation noise {($\mathbb{R}^{d_y}$)}\\
 \hline
 $f$ & Nonlinear function \\
 \hline
 $\hat{\mathbf{x}}$ & Inferred latent state {($\mathbb{R}^{d_x}$)}\\
 \hline
 % $g$ & Nonlinear observation function \\
 % \hline
   $q(\cdot)$ & Variational distribution \\
   \hline
 $\mathcal{F}$ & Variational free energy \\
  \hline
 $\bm{\bm{\epsilon}}^\bfx$ & Dynamics prediction error {($\mathbb{R}^{d_x}$)}\\ 
 \hline
 $\bm{{\epsilon}}^\bfy$ & Observation prediction error {($\mathbb{R}^{d_y}$)}\\
 \hline
 $\eta$ & learning rate \\
 \hline
 $\Delta t$ & Inference step size \\
 \hline
  $\bm{\Sigma}_k$ & Posterior variance {($\mathbb{R}^{d_x \times d_x}$)}\\ 
 \hline
 $\bfK$ & Kalman Gain matrix {($\mathbb{R}^{d_x \times d_y}$)}\\
  \hline
 % $\KL[\cdot || \cdot ]$ & KL divergence \\
 %   \hline
 % $\phi$ & parameters of variational distribution \\
 % [1ex] 
 % \hline
\end{tabular}
\label{notation_table}
\end{center}
\end{table}

\subsection{Model foundations: HMM and free energy}

%We situate our explanation of 
The conceptual level of our model is grounded within the Bayesian Brain paradigm \citep{knill2004bayesian,pouget2013probabilistic,friston2005theory,friston2012free}. Specifically, we assume that the problem of perception is fundamentally an \emph{inference} problem, where there exists some real world ``out there'', from which we only receive noisy and distorted sensory input. We assume that the task of perception is to untangle and counteract the noise in order to reconstruct the real (but hidden) state of the world given only our sensory observations. Thus, mathematically, we can represent the problem of perception as trying to infer a series of latent states of the world $\mathbf{x}_k$ ($k=1,2,\dots$)\footnote{Throughout this thesis, $k$ will be used to denote time steps in the sensory domain, and $t$ will be used for the rate of change of neural activities as in $dx/dt$.} from their corresponding sensory observations $\mathbf{y}_k$ ($k=1,2,\dots$). 
%Classical, or `static' predictive coding assumes that the observations $\mathbf{y}_k$ are independent and thus does not provide a model of the underlying dynamics generating the observations i.e., at each time step $k$, we independently infer the hidden $\mathbf{x}_k$ from $\mathbf{y}_k$, ignoring the past and the future \citep{bogacz2017tutorial}. To take into account the dynamical and temporal aspects of the observations, 
We assume that the underlying graphical structure of our tPC network is an HMM, where the hidden states $\mathbf{x}_k$ follow a Markov chain. That is, the current hidden state of the world only depends upon the previous hidden state. 
% The HMM also assumes that
Also, the current observation is generated by the current state of the world only, with no dependence on past noisy observations. Fig  \ref{fig:gen_model_schematic} shows the generative process of the tPC where, for generality, we have also added ``control'' inputs $\mathbf{u}$ which can be thought of as known inputs to the system at every time step (such inputs are included in the Kalman filter that we will later use as a benchmark). This control input is useful to model systems where there are known external forces acting on the system (such as an agent's own actions) which we do not necessarily want to model simply as part of the environmental dynamics.

\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{figs/chap4/hmm.pdf}
    \caption{{\bf Graphical model of the generative process assumed by temporal predictive coding.} $\mathbf{x}_k$ correspond to hidden states, $\mathbf{y}_k$ to observations, and $\mathbf{u}_k$ to control inputs. Circles denote latent variables, squares denote observations, and arrows denote conditional dependence of the variables (the absence of an arrow indicates conditional independence).}
    \label{fig:gen_model_schematic}
\end{figure}

From the generative model in Fig \ref{fig:gen_model_schematic}, we can also write out specific equations for the dynamics of the states and observations in what is called the state-space representation:
\begin{eqnarray}\label{eq:general_transition}
    \mathbf{x}_k = \bfW f(\mathbf{x}_{k-1}) + \bfB\mathbf{u}_{k} + \bm{\omega}_x,
\end{eqnarray}
\begin{eqnarray}\label{eq:general_emission}
    \mathbf{y}_k = \bfF f(\mathbf{x}_k) + \bm{\omega}_y
\end{eqnarray}
where $\bfW$ is the matrix that transitions the previous hidden state $\mathbf{x}_{k-1}$ to $\mathbf{x}_{k}$, $\bfB$ is the matrix governing how the control input $\mathbf{u}_k$ affects the current hidden state, and $\bfF$ is the ``emission'' matrix that determines how the observation $\mathbf{y}_k$ is generated from the hidden state. $f$ is a function transforming $\mathbf{x}_{k}$ that may be nonlinear. The above state-space representation also includes sources of noise $\bm{\omega}_x$ and $\bm{\omega}_y$. We will assume a white Gaussian noise model such that
$\bm{\omega}_x \sim \mathcal{N}(0, \bm{\Sigma}_x)$ and $\bm{\omega}_y \sim \mathcal{N}(0, \bm{\Sigma}_y)$ are zero-mean Gaussian random variables with covariance matrices $\bm{\Sigma}_x$, $\bm{\Sigma}_y$ \footnote{In the control theory literature, the process noise $\bm{\Sigma}_x$ is often denoted as $\mathbf{Q}$ and the observation noise $\bm{\Sigma}_y$ as $\mathbf{R}$. We instead follow the convention that has arisen within the predictive coding literature, which we also believe is more straightforward, by calling them $\bm{\Sigma}_x$ and $\bm{\Sigma}_y$, since it makes explicit that these variances are simply the variances of the two distributions composing the generative model.
%(likelihood and transition distributions respectively).
}. Therefore, $\mathbf{x}_k$ and $\mathbf{y}_k$ can be considered as random variables that follow the Gaussian distributions:
\begin{eqnarray}\label{eq:gaussian-latent}
    \mathbf{x}_k | \mathbf{x}_{k-1},\mathbf{u}_k \sim \mathcal{N}(\bfW f(\mathbf{x}_{k-1}) + \bfB\mathbf{u}_{k}, \bm{\Sigma}_x),
\end{eqnarray}
\begin{eqnarray}\label{eq:gaussian-obs}
    \mathbf{y}_k | \mathbf{x}_k \sim \mathcal{N}(\bfF f(\mathbf{x}_k), \bm{\Sigma}_y).
\end{eqnarray}

The objective of our model is to obtain an estimate $\hat{\mathbf{x}}_k$ of the current $\mathbf{x}_k$, given the previously estimated $\hat{\mathbf{x}}_{k-1}$, current noisy observation $\mathbf{y}_k$ and control input $\mathbf{u}_k$. This objective can be expressed as estimating the posterior distribution $p(\mathbf{x}_k | \mathbf{y}_k, \hat{\mathbf{x}}_{k-1}, \mathbf{u}_k)$.
%, which takes an analogous form of the objective used in static predictive coding and more generally, free energy principle \citep{friston2005theory,bogacz2017tutorial}. 
%Similar to the static case
If function $f$ is nonlinear, an analytic expression for the posterior cannot be found, and thus variational inference \citep{neal1998view,beal2003variational,ghahramani2000graphical} is utilized to find an approximate of the posterior. Similar to what was introduced in Chapter \ref{chap:intro} for the static case, we assume that there exists an approximate posterior $q(\mathbf{x}_k)$ and we seek an approximate posterior that is as close as possible to the ``true'' posterior $p(\mathbf{x}_k | \mathbf{y}_k, \hat{\mathbf{x}}_{k-1}, \mathbf{u}_k)$, and variational inference finds such a $q(\mathbf{x}_k)$ by minimizing an upper bound on the divergence between the true and approximate posteriors called the variational free energy. If we make an additional simplifying assumption that $q(\mathbf{x}_k)$ follows a Gaussian distribution with its density highly concentrated at its mean,  the free energy $\mathcal{F}_k$ becomes approximately equal to \citep{buckley2017free,bogacz2017tutorial}:
\begin{eqnarray}\label{eq:tPC-FE}
    \begin{aligned}
    \mathcal{F}_k &= - \log p(\mathbf{y}_k, \mathbf{x}_k | \hat{\mathbf{x}}_{k-1}, \mathbf{u}_k) \\ & = - \log p(\mathbf{y}_k| \mathbf{x}_k) p(\mathbf{x}_k | \hat{\mathbf{x}}_{k-1}, \mathbf{u}_k)
\end{aligned}
\end{eqnarray}
and it is also sufficient to estimate the mode of the approximate posterior (which is the same as its mean) instead of estimating the whole distribution with this assumption. Furthermore, with the Gaussian assumptions underlying Eqs. \ref{eq:gaussian-latent} and \ref{eq:gaussian-obs}, the free energy becomes (we choose to omit the other terms in the multivariate Gaussian density as they do not affect the optimization over $\mathbf{x}_k$ and $\bfW$, $\bfB$, $\bfF$):
% \LabelLines{Omit}{\Omit}
\begin{eqnarray}\label{eq:tpc-objective}
\begin{aligned}
    \mathcal{F}_k &= \frac{1}{2}(\mathbf{y}_k - \bfF f(\mathbf{x}_k))^{\top}\bm{\Sigma}_y^{-1} (\mathbf{y}_k - \bfF f(\mathbf{x}_k)) \\&+ \frac{1}{2}(\mathbf{x}_k - \bfW f(\hat{\mathbf{x}}_{k-1}) - \bfB \mathbf{u}_k)^{\top} \bm{\Sigma}_x^{-1} (\mathbf{x}_k - \bfW f(\hat{\mathbf{x}}_{k-1})- \bfB \mathbf{u}_k) + C
\end{aligned}
\end{eqnarray}

Importantly, we can express this objective as the sum of squared \emph{prediction errors} weighted by their inverse covariances (which are often called precisions in the predictive coding literature). In this model, there are two kinds of prediction errors -- ``sensory'' prediction errors which are the difference between the observation and predicted observation $\mathbf{y}_k - \bfF f(\mathbf{x}_k)$ and ``temporal'' prediction errors which are the difference between the inferred current state and the current state predicted from the previous state $\mathbf{x}_k - \bfW f(\hat{\mathbf{x}}_{k-1})- \bfB \mathbf{u}_k$. Thus, by finding an estimate $\hat{\mathbf{x}}_k = \argmin_{\mathbf{x}_k} \mathcal{F}_k$, we effectively minimize the squared sum of these prediction errors while the precision matrices serve to weight the importance of the sensory and temporal prediction errors in accordance with their intrinsic variance (so highly variable prediction errors are weighted less). After the minimization finishes, the estimated $\hat{\mathbf{x}}_k$ can then be used to estimate the hidden state at the next step $k+1$.

\subsection{{Inference and learning algorithm}}

With the objective function in Eq. \ref{eq:tpc-objective}, it is then possible to derive an iterative algorithm to perform its minimization via gradient descent. Similar to static predictive coding \citep{rao1999predictive}, the gradient descent is performed on two sets of values: the hidden states $\mathbf{x}_k$ and the weight parameter matrices $\bfW$, $\bfB$ and $\bfF$. As we will show in the next subsection, the former can be implemented as neural responses and the latter can be implemented as synaptic connections in a neural circuit in a similar way to static predictive coding. For the hidden state $\mathbf{x}_k$, its dynamics follow:
\begin{eqnarray}\label{eq:x-GD}
    \tau\frac{d \mathbf{x}_k}{dt} = - \frac{\partial \mathcal{F}_k}{\partial \mathbf{x}_k},
\end{eqnarray}
{where $\tau$ is the time constant of the neurons}. At convergence, we can say that the equilibrium value $\hat{\mathbf{x}}_k$ represents the optimal inference about the mean of the true hidden state of the world. It is worth mentioning that we now introduced two different indices $k$ and $t$ for distinct time scales: for discrete steps $k$ at which the observations arise and continuous real time $t$ in which computations are made within the model, and we will discuss the relationship between them in detail in the next subsection.

% For simplicity at the initial part of the work we will make an additional assumption that the likelihood and transition functions are \emph{linear} which corresponds to the assumption that the external world dynamically evolves according to a linear process and that the observations we receive are a linear map of the hidden state, as defined in Eqs. \ref{eq:ABC}
%such that we replace the likelihood function $g$ with the matrix $C$ and the transition function $f$ with the matrices $A$ and $B$ 
% (later on in the work, when we consider nonlinear inference, this assumption will be relaxed). 
%We can also assume a `control input' to the system $u$ where the effect of $u$ on the system is parametrized by the matrix $B$.  
To derive the exact expression for the inferential dynamics, we can define the precision-weighted state and observation prediction errors as $\bm{\epsilon}^{\bfx}$ and $\bm{\epsilon}^{\bfy}$ respectively as follows:
\begin{eqnarray}\label{eq:err-y}
    \bm{\epsilon}^{\bfy} =  \bm{\Sigma}_y^{-1} ( \mathbf{y}_k - \bfF f(\mathbf{x}_k) ),
\end{eqnarray}
\begin{eqnarray}\label{eq:err-x}
    \bm{\epsilon}^{\bfx} = \bm{\Sigma}_x^{-1} ( \mathbf{x}_k - \bfW f(\hat{\mathbf{x}}_{k-1}) - \bfB \mathbf{u}_k ).
\end{eqnarray}

We can then write Eq. \ref{eq:x-GD} as:
\begin{eqnarray}\label{eq:dynamics_x}
    \tau\frac{d \mathbf{x}_k}{dt} = - \frac{\partial \mathcal{F}_k}{\partial \mathbf{x}_k} = - \bm{\epsilon}^{\bfx} + f'(\mathbf{x}_k) \odot \bfF^{\top} \bm{\epsilon}^{\bfy},
\end{eqnarray}
where $\odot$ denotes the element-wise product between vectors. The dynamics have contributions from both the sensory and the temporal prediction errors, and the contribution of the sensory prediction error is ``mapped backwards'' through the transpose matrix $\bfF^{\top}$. 

Similarly, if we assume that the $\bfW$, $\bfB$ and $\bfF$ parameter matrices are learnable, we can derive update rules also following gradient descent on $\mathcal{F}_k$:
\begin{eqnarray}\label{eq:A_C_updates}
    \begin{aligned}
    \Delta \bfW &= - \eta \frac{\partial \mathcal{F}_k}{\partial \bfW} = \eta \bm{\epsilon}^{\bfx} f(\hat{\mathbf{x}}_{k-1})^{\top}, \\
    \Delta \bfB &= - \eta \frac{\partial \mathcal{F}_k}{\partial \bfB} = \eta \bm{\epsilon}^{\bfx}\mathbf{u}_k^{\top}, \\
    \Delta \bfF &= - \eta \frac{\partial \mathcal{F}_k}{\partial \bfF} = \eta \bm{\epsilon}^{\bfy} f(\mathbf{x}_k)^{\top}.
    \end{aligned}
\end{eqnarray}

In the above equations, $\eta$ denotes a scalar learning rate. As we will see below, the iterative updates will correspond to local Hebbian plasticity in the neural implementation of the model.

Typically, if the $\bfW$, $\bfB$, and $\bfF$ matrices are learned, then the matrices are updated by a single step according to Eqs. \ref{eq:A_C_updates} after the $\mathbf{x}_k$ has already converged and using the equilibrium values $\hat{\mathbf{x}}_k$. This is because it is often assumed that these variables represent more slowly changing variables in the real world. Moreover, in the neural implementation, these matrices are often assumed to be implemented by synaptic strengths which typically change slowly while the $\hat{\mathbf{x}}$ variables are typically mapped to neural firing rates, which change quickly.  

The process of learning a tPC model is shown in Algorithm~\ref{alg:tPC}. The algorithm assumes that sensory observations $\mathbf{y}_k$ arrive at discrete steps $k$. At each step $k$, we first initialize the $\mathbf{x}_k$ values with the previous estimated $\hat{\mathbf{x}}_{k-1}$ value from the last step.
%passed through the dynamics function. 
Then, we iterate Eq.\ref{eq:dynamics_x} until convergence for a given observation $\mathbf{y}_k$. Upon convergence our $\hat{\mathbf{x}}_k$ becomes our best estimate of the true state of the world. Given this $\hat{\mathbf{x}}_k$, we can update the $\bfW$, $\bfB$, and $\bfF$ matrices using Eqs.\ref{eq:A_C_updates}. 

\begin{algorithm}
\caption{Single training epoch for temporal predictive coding}\label{alg:tPC}
\begin{algorithmic}[1]
\State $K$: Discrete steps of observations
% $T$: Continuous inference time until $\mathbf{x}_k$ converges\\
\For{$k=1,...,K$}
    \State Initialize $\mathbf{x}_k$ with previously inferred $\hat{\mathbf{x}}_{k-1}$ 
    \While{$\mathbf{x}_k$ not converged} 
        \State Perform inference by updating $\mathbf{x}_k$ (Eq~\ref{eq:dynamics_x})
    \EndWhile
    \State Update weight matrices $\Delta \bfW$, $\Delta \bfB$, $\Delta \bfF$ (Eq~\ref{eq:A_C_updates}) using inferred $\hat{\mathbf{x}}_k$
\EndFor
\end{algorithmic}
\end{algorithm}


\subsection{Neural circuit implementation}

The update rules and dynamics we have derived in Eqs.\ref{eq:dynamics_x} and \ref{eq:A_C_updates} can be mapped to a recurrent predictive coding network architecture with biologically plausible Hebbian learning rules. In this section, we present two examples of networks implementing the algorithm exactly, and then a simplified network that approximates the algorithm.

Fig \ref{fig:recurrent_multistep}A presents a potential example of how the update rules we have derived can be implemented in a neural network with only local and Hebbian plasticity, which is similar to the standard predictive coding network \citep{rao1999predictive}. In this network, observations $\mathbf{y}_k$ enter at the lowest level, and cause sensory prediction errors $\bm{\epsilon}^{\bfy}$ as the observations are met by top-down predictions. These prediction errors are explicitly represented by the activities of  ``prediction error neurons''. These prediction error neurons receive top-down inhibitory connections from ``value neurons'' in the layer above which, through their firing rates, represent the inferred posterior values $\hat{\mathbf{x}}_k$. 
Similarly, at the layer above there are additional prediction error neurons that represent the difference between the current activity and that predicted based on the previous inference mapped through the dynamics function (Eq.\ref{eq:err-x}), and we assume that there are dedicated neurons that ``memorize'' latent activities $\hat{\mathbf{x}}_{k-1}$ inferred at the previous discrete time step, which is used to make a prediction of the current latent activities and reloaded at the end of inference at each time step.

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth]{figs/chap4/Fig2.pdf}
    \caption{{\bf Possible neural implementations of temporal predictive coding.} A: Potential neural circuit implementing the iterative recurrent predictive coding algorithm. For simplicity, we have depicted each neural layer as possessing only two neurons. 
    B: Version of the model where the prediction errors are represented by the difference in membrane potential in soma and at apical dendrites (depicted as {ellipses}).
    C: Neural circuitry required to implement the single-iteration predictive coding algorithms. 
    %The key difference to the previous discrete multistep algorithm is that there no longer needs to be 
    This model no longer includes a separate set of neurons explicitly storing the estimate of the previous timestep, but instead, the temporal prediction errors are computed naturally through recurrent connections. For simplicity, we omitted the control inputs $\bfB\mathbf{u}_k$, which can be implemented in a similar way to the recurrent inputs $\bfW\hat{\mathbf{x}}_{k-1}$ to the error neurons or apical dendrites.
    }
     \label{fig:recurrent_multistep}
\end{figure}

There are several important aspects to note about this model. First, all the required update rules can be implemented using purely local information. The dynamics of the hidden state estimates $\mathbf{x}_k$ (Eq.\ref{eq:dynamics_x}) can be reproduced locally since the value neurons receive inhibitory connections from the prediction error neurons $\bm{\epsilon}^{\bfx}$ at their layer as well as excitatory connections from the prediction error neurons $\bm{\epsilon}^{\bfy}$ at the layer below. Similarly, the prediction errors $\bm{\epsilon}^{\bfy}$ can be computed according to Eq.\ref{eq:err-y} as the corresponding prediction error neurons receive excitatory input $\mathbf{y}_k$ and inhibition from neurons encoding $\mathbf{x}_k$. 
% achieved through additional inhibitory connections between prediction error neurons (not shown for simplicity) as described by \citep{bogacz2017tutorial}, and it was recently demonstrated that the inverse of the covariance matrix can also be reparameterized as a recurrent weight matrix to circumvent the implausibility of the inversion \citep{tang2023recurrent}. 
The prediction error $\bm{\epsilon}^{\bfx}$ can be computed analogously.
The update rules for the $\bfW$, $\bfB$, and $\bfF$ weight matrices (Eqs.\ref{eq:A_C_updates}) are also precisely Hebbian, since they are outer products between the prediction errors and the value neurons of the layer above which, crucially, are also precisely the pre and post-synaptic activities of the neurons where the synapses implementing these weight matrices are located. {Moreover, we empirically demonstrate in the Results sections that scaling by the inverse covariance matrices $\bm{\Sigma}_x^{-1}$ and $\bm{\Sigma}_y^{-1}$ could be encoded in the learned $\bfW$ and $\bfF$ matrices, similarly as it has been shown in static predictive coding models (Chapter \ref{chap:am}; \citet{tang2023recurrent}). Thus, the tPC model can represent the covariance matrices implicitly in its synaptic weights, without needing to implement them in explicit synaptic weights.} 
% \newcommand{\sometexts}{
    {The nonlinear function $f(\cdot)$ can be implemented in this circuit following the way specified in \citet{whittington2017approximation}, through local inhibitory neurons.}
% }
% \LabelLines{sometexts}{\sometexts{}}
%The local and Hebbian synaptic plasticity achievable with predictive coding networks like this is elegant and rarely obtained in other network architectures \citep{crick1989recent,bengio2017stdp,scellier2017equilibrium,sacramento2018dendritic,millidge2020predictive}.%The expressions for prediction errors in Eqs.\ref{prediction_errors_definition} and \ref{eq:error_x} are scaled by inverses of covariance matrices. However, if the covariance matrices are diagonal, such scaling would correspond to a multiplication by a scalar value ($\bm{\Sigma}_{i,i}^{-1}$) which could be encoded in the gain of prediction error neurons. %For simplicity, in the simulations of this work we consider only cases of diagonal covariance matrices, and a possible neural implementation with non-diagonal matrices is described by \citep{bogacz2017tutorial}.

The network shown in Fig \ref{fig:recurrent_multistep}A follows a standard predictive coding architecture, but it could be simplified because the prediction error neurons encoding $\bm{\epsilon}^{\bfx}$ only project to corresponding neurons encoding $\bfx$, and we could thus borrow the idea of dendritic computing similar to the model of \citet{sacramento2018dendritic}. In particular, substituting Eq.\ref{eq:err-x} into Eq. \ref{eq:dynamics_x}, we obtain:
\begin{eqnarray}
    \tau\frac{d \mathbf{x}_k}{dt} = - \mathbf{x}_k + \bfW f(\hat{\mathbf{x}}_{k-1}) + \bfB \mathbf{u}_k + f'(\mathbf{x}_k) \odot \bfF^{\top} \bm{\epsilon}^{\bfy} ,
    \label{eq:dendritic}
\end{eqnarray}
where, for clarity of explanation, we omitted the precision $\bm{\Sigma}_x^{-1}$ that can be encoded in $\bfW$. By writing the dynamical equation in this way, we assume that there is no building block within the model that encodes the error explicitly; rather, the apical dendrite will encode the inputs $\bfW f(\hat{\mathbf{x}}_{k-1}) + \bfB \mathbf{u}_k$ and send the signal to the soma of the pyramidal neuron. This dendritic signal excites the soma and drives the inferential dynamics (Eq. \ref{eq:dendritic}), together with the decay $- \mathbf{x}_k$ that is intrinsic to the soma and feedback signals from the observation layer. The corresponding neural implementation is shown in Fig \ref{fig:recurrent_multistep}B. Although the architecture of the network becomes simpler, learning parameters $\bfW$ and $\bfB$ is less straightforward because the prediction error $\bm{\epsilon}^{\bfx}$ is not explicitly represented in the activity of any neurons in the network. 
Nevertheless, the prediction error $\bm{\epsilon}^{\bfx}$ is equal to the difference between the neural activity and the membrane potential in the apical dendrite, and it has been proposed that such difference drives plasticity of synapses on the apical dendrite \citep{urbanczik2014learning}. 
Since both the neural activity and the membrane potential are encoded within the same neuron, it is plausible that their difference could be computed within the neuron (e.g. the information on the neural activity could be brought to the synapses on apical dendrites via backpropagating action potentials).  
Such a signal could then drive local synaptic plasticity to learn the $\bfW$ and $\bfB$ matrices. 

While simulating the model, we update the state estimates by numerically solving Eq. \ref{eq:dendritic} using the Euler method, i.e. we calculate the state estimates for every interval $\Delta t$:
\begin{eqnarray}\label{eq:inf-discrete}
    \mathbf{x}_k(t+\Delta t) = \mathbf{x}_k(t) + \frac{\Delta t}{\tau} \left[ - \mathbf{x}_k(t) + \bfW f(\hat{\mathbf{x}}_{k-1}) + \bfB \mathbf{u}_k  + f'(\mathbf{x}_k(t))\odot \bfF^{\top} \bm{\epsilon}^{\bfy} \right] .
    \label{eq:discrete}
\end{eqnarray}
{where $\frac{\Delta t}{\tau}$ is effectively the step size of the discretized inference that we tune in our simulations.} The above expression highlights that the algorithm has two nested timescales -- firstly there is the ``external'' timescale which is where sensory inputs $\mathbf{y}_k$ are received in a sequence of steps we index by subscript $k$. Then, for each external step, there is an internal inference of the hidden state that is numerically implemented as a set of recurrent iterations within that external step, which we denote as $t$. In such a nested framework, the implementations in Figs \ref{fig:recurrent_multistep}A and \ref{fig:recurrent_multistep}B
% \newcommand{\timescale}{
    {need to}
    store and hold fixed the state estimate of the previous step while the iterative inference of the state estimate for the current step is ongoing.
% }
% \LabelLines{timescale}{\timescale{}} 
Specifically, the estimate from the previous step $\hat{\mathbf{x}}_{k-1}$ needs to be held fixed throughout the iterative procedure while the actual current values $\mathbf{x}_k(t)$ vary in order to find a balance $\hat{\mathbf{x}}_k$ between the demands of matching the prediction from the last step and also the current observation (Algorithm \ref{alg:tPC}). Once the iterations are complete for a step, the new value of $\hat{\mathbf{x}}_k$ needs to be loaded into the memory and stored as the last step for the next set of iterations. 
In situations where the observations are separated in time, it is known that neurons are able to store the representation of stimuli presented a few seconds earlier in their activity \citep{rainer1998selective}.
In the case where sensory input arrives continuously, the two timescales coincide, and there may be no time for inference between steps of sensory input. In that latter case, the algorithm can be adapted to remove the issue of nested timescales without unduly harming filtering performance by simply using a \emph{single} iteration of the internal inference for each external step. This means that there is effectively no inner loop of the algorithm anymore, since the inner loop consists of just a single iteration. This makes the algorithm fully online in the sense that it receives a new sensory input for every step. 
In particular, note that by equating time indices $\Delta k =\Delta t$, we obtain:
\begin{eqnarray}
    \hat{\mathbf{x}}_{k} = \mathbf{x}_k + \frac{\Delta t}{\tau} \left[ - \mathbf{x}_k + \bfW f(\hat{\mathbf{x}}_{k-1}) + \bfB \mathbf{u}_k  + f'(\mathbf{x}_k)\odot \bfF^{\top} \bm{\epsilon}^{\bfy} \right].
\end{eqnarray}

A diagram of a potential neural circuit implementing this single-step algorithm is presented in Fig \ref{fig:recurrent_multistep}C. This network no longer includes neurons storing past inferences.
Instead, the temporal prediction errors are computed solely using recurrent connections labelled $\bfW$, which are now assumed to introduce a temporal delay of one step. 

The advantage of this approach is that it eschews the challenge of storing and loading the memory of the last step; instead, this memory can be dynamically maintained across a single external step simply through recurrent connections via their intrinsic synaptic delays. The disadvantage of this approach is that the update rules of the algorithm were derived as gradient descent on an objective function, and this approach is equivalent to taking only a single gradient step for each example. Clearly, in many cases, such an algorithm simply would not work because a single step is nowhere near enough to approach the optimum. However, there are two features of the problem that ameliorate much of this difficulty in practice. 
The first is that, when $f$ is a linear function, the objective is actually convex, and thus the loss landscape is extremely well-behaved. This allows for the use of relatively high integration steps to move large distances in a single, or a few steps, without fear of overshooting the optimum or running into divergences. The other factor is due to the nature of the external world: typically, visual scenes change relatively slowly on a microsecond-by-microsecond level, and thus the optimal estimated hidden state in a single step is likely close to the optimum hidden state for the next. In this case, since we initialize the inference of each step with the optimum of the last step, this will usually be close to the optimum for the current step as well, thus meaning that the algorithm simply does not have to make many iterations to achieve the optimum since it already starts close by. 
In the next section, we show that, in practice, on standard tracking and filtering tasks, these two factors can often simplify the inference problem enough that this single iteration approach often works successfully, although it usually does not perform quite as well as multi-step methods.

\section{Results}

The results of this work are partitioned into a theoretical section and experimental sections. In the theoretical section, we examine the relationship between the tPC model and Kalman filtering and demonstrate that the tPC network, under certain assumptions, is equivalent to a Kalman filter with a fixed posterior variance.
%of $\bm{\Sigma}_k=0$. 
In simulations, we demonstrate that, despite not correctly representing the posterior variance, the tPC network nevertheless exhibits strong and robust tracking performance on filtering tasks while also being capable of online system identification through the learning of the $\bfW$ and $\bfF$ matrices, unlike the Kalman filter. 
{Moreover, we show that when tPC is trained with natural movies, the simulated neurons in the latent layer develop motion-sensitive receptive fields resembling those of neurons in the primary visual cortex.}


\subsection{Relationship to Kalman filtering}

Here we show that both Kalman filtering \citep{kalman1960new} and tPC can be derived as special cases of the \emph{Bayesian filtering} problem. 
%However, since they make different assumptions on the prior information at each time step, predictive coding differs from Kalman filtering theoretically. 
Bayesian filtering concerns the problem of inferring the \emph{sequence} of hidden causes $\mathbf{x}_1,...,\mathbf{x}_K$ of the observations $\mathbf{y}_1,...,\mathbf{y}_K$. This problem can be effectively factorised into a sequence of online inference problems i.e., inferring the hidden state $\mathbf{x}_k$ at time step $k$, given the whole history of observations $\mathbf{y}_1,...,\mathbf{y}_k$ \citep{jazwinski2007stochastic,stengel1986stochastic}. For simplicity of notation, we denote $\mathbf{x}_{1:k}=\mathbf{x}_1,...,\mathbf{x}_k$ and $\mathbf{y}_{1:k}=\mathbf{y}_1,...,\mathbf{y}_k$. The Bayesian filtering problem can thus be formulated as inferring the following posterior distribution:
\begin{eqnarray}
    p (\mathbf{x}_k | \mathbf{y}_{1:k}).
\end{eqnarray}

We show in Appendix \ref{app:chap4} that this posterior distribution is proportional to:
\begin{eqnarray}\label{eq:bayesian-filtering}
    p (\mathbf{x}_k | \mathbf{y}_{1:k}) \propto p(\mathbf{y}_k |\mathbf{x}_k)\int p(\mathbf{x}_k|\mathbf{x}_{k-1})p(\mathbf{x}_{k-1} | \mathbf{y}_{1:k-1}) d\mathbf{x}_{k-1}
\end{eqnarray}
where the integral is effectively the marginal distribution $p(\mathbf{x}_k|\mathbf{y}_{1:k-1})$ of the joint $p(\mathbf{x}_k, \mathbf{x}_{k-1} | \mathbf{y}_{1:k-1})$ and can be considered as the prior on $\mathbf{x}_k$. Notice that the term $p(\mathbf{x}_{k-1} | \mathbf{y}_{1:k-1})$ is exactly the posterior inferred from the previous time step $k-1$, making Bayesian filtering a recursive method \citep{chen2003bayesian}. As a special case of Bayesian filtering, Kalman filtering assumes that the conditional distributions $p(\mathbf{y}_k |\mathbf{x}_k)$ and $p(\mathbf{x}_k|\mathbf{x}_{k-1})$ can be parameterized \emph{linearly} as follows:
\begin{eqnarray}\label{eq:linear-assumption}
    \mathbf{y}_k |\mathbf{x}_k \sim \mathcal{N}(\bfF\mathbf{x}_{k},\bm{\Sigma}_y); \ \ \mathbf{x}_k | \mathbf{x}_{k-1} \sim \mathcal{N}(\bfW\mathbf{x}_{k-1}+\bfB\mathbf{u}_k,\bm{\Sigma}_x).
\end{eqnarray}

Further, it assumes that the posterior estimated at the previous step $k-1$ follows:
\begin{eqnarray}\label{eq:KF-gaussian}
    \mathbf{x}_{k-1} | \mathbf{y}_{1:k-1} \sim \mathcal{N}(\hat{\mathbf{x}}_{k-1}, \bm{\Sigma}_{k-1})
\end{eqnarray}
where $\hat{\mathbf{x}}_{k-1}$ is the MAP estimate from the previous step $k-1$ and is the mode (or mean) of the Gaussian posterior with covariance $\bm{\Sigma}_{k-1}$ at step $k-1$. Under the above Gaussian assumptions, the prior $p(\mathbf{x}_k|\mathbf{y}_{1:k-1})$ on $\mathbf{x}_k$ (i.e., the integral in Eq. \ref{eq:bayesian-filtering}) can be written as \citep{bishop2006pattern}:
\begin{eqnarray}
    p(\mathbf{x}_k|\mathbf{y}_{1:k-1}) = \mathcal{N}(\bfW\hat{\mathbf{x}}_{k-1}+\bfB\mathbf{u}_k, \bfW\bm{\Sigma}_{k-1}\bfW^{\top}+\bm{\Sigma}_x).
\end{eqnarray}

Kalman filtering then performs maximum a posteriori (MAP) to find $\hat{\mathbf{x}}_k$:
\begin{eqnarray}\label{eq:KF-objective}
    \begin{aligned}
        \hat{\mathbf{x}}_k 
        &= \argmax_{\mathbf{x}_k} \log p (\mathbf{x}_k | \mathbf{y}_{1:k}) \\
        &= \argmin_{\mathbf{x}_k} (\mathbf{y}_k - \bfF{\mathbf{x}}_k)^{\top}\bm{\Sigma}_y^{-1}(\mathbf{y}_k - \bfF{x}_k) \\&+ (\mathbf{x}_k - \bfW\hat{\mathbf{x}}_{k-1} - \bfB\mathbf{u}_k)^{\top} (\bfW\bm{\Sigma}_{k-1}\bfW^{\top}+\bm{\Sigma}_x)^{-1}({\mathbf{x}}_k - \bfW\hat{\mathbf{x}}_{k-1}- \bfB\mathbf{u}_k).
    \end{aligned}
\end{eqnarray}

Since all the transformation functions in this optimization problem are linear, an analytical expression for $\hat{\mathbf{x}}_k$ can be derived, which will result in the well-known algorithm for Kalman filtering i.e., the ``projection'':
\begin{eqnarray}\label{eq:KF-projection}
\begin{aligned}
    \hat{\mathbf{x}}^-_k &= \bfW \hat{\mathbf{x}}_{k-1} + \bfB \mathbf{u}_{k} \\
    \bm{\Sigma}^-_k &= \bfW \bm{\Sigma}_{k-1} \bfW^{\top} + \bm{\Sigma}_x 
\end{aligned}
\end{eqnarray}
and the ``correction'' step where we then incorporate the new information we have received from the environment to correct our estimates:
\begin{eqnarray}\label{eq:KF-correction}
\begin{aligned}
    \hat{\mathbf{x}}_k &= \hat{\mathbf{x}}^-_k + \bfK(\mathbf{y}_k - \bfF \hat{\mathbf{x}}^-_k) \\
    \bm{\Sigma}_k &= (\bfI - \bfK \bfF)\bm{\Sigma}^-_k \\
    \bfK &= \bm{\Sigma}^-_k \bfF^{\top} \big[ \bfF \bm{\Sigma}^-_k \bfF^{\top} + \bm{\Sigma}_y]^{-1}
\end{aligned}
\end{eqnarray}
where $\bfK$ is known as the Kalman Gain matrix and is central to the Kalman filter update rules for the estimated mean and variance in the correction step \citep{welch1995introduction}. $\hat{\mathbf{x}}_k$ and $\bm{\Sigma}_k$ are then our estimated mean and covariance of the posterior Gaussian distribution. The derivation of the projection and correction rules can be found in prior works \citep{chen2003bayesian,millidge2021neural}, while we also provide the derivation in Appendix \ref{app:chap4} as to compare to the update rules of our tPC model, which is demonstrated below.

Our tPC model also aims to solve the Bayesian filtering problem in Eq. \ref{eq:bayesian-filtering}, and it can also make the linear and Gaussian assumptions underlying Eq. \ref{eq:linear-assumption}. Notice that Eq. \ref{eq:linear-assumption} is identical to Eqs. \ref{eq:gaussian-latent} and \ref{eq:gaussian-obs}, but with a linear $f$. tPC differs from Kalman filtering by making a different assumption on the distribution on the previous-step posterior $p(\mathbf{x}_{k-1} | \mathbf{y}_{1:k-1})$. Instead of assuming it as a Gaussian distribution in Eq. \ref{eq:KF-gaussian}, it assumes:
\begin{eqnarray}\label{eq:PC-dirac}
    \mathbf{x}_{k-1} | \mathbf{y}_{1:k-1} \sim \delta(\mathbf{x}_{k-1} - \hat{\mathbf{x}}_{k-1})
\end{eqnarray}
where $\delta(\mathbf{x}_{k-1} - \hat{\mathbf{x}}_{k-1})$ denotes a Dirac distribution with its density concentrated at $\hat{\mathbf{x}}_{k-1}$. The prior on $\mathbf{x}_k$ for predictive coding thus becomes:
\begin{eqnarray}
    p(\mathbf{x}_k|\mathbf{y}_{1:k-1}) = p(\mathbf{x}_k | \hat{\mathbf{x}}_{k-1}) = \mathcal{N}(\bfW\hat{\mathbf{x}}_{k-1}+\bfB\mathbf{u}_k, \bm{\Sigma}_x)
\end{eqnarray}
since the density of $p(\mathbf{x}_{k-1} | \mathbf{y}_{1:k-1})$ is concentrated at $\hat{\mathbf{x}}_{k-1}$. The MAP estimation of $\hat{\mathbf{x}}_k$ performed by predictive coding is thus:
\begin{eqnarray}\label{eq:PC-objective}
    \begin{aligned}
        \hat{\mathbf{x}}_k 
        &= \argmax_{\mathbf{x}_k} \log p (\mathbf{x}_k | \mathbf{y}_{1:k}) \\
        &= \argmin_{\mathbf{x}_k} (\mathbf{y}_k - \bfF{\bfx}_k)^{\top}\bm{\Sigma}_y^{-1}(\mathbf{y}_k - \bfF{\bfx}_k) \\&+ (\bfx_k - \bfW\hat{\mathbf{x}}_{k-1} - \bfB\mathbf{u}_k)^{\top} \bm{\Sigma}_x^{-1}(\bfx_k - \bfW\hat{\mathbf{x}}_{k-1}- \bfB\mathbf{u}_k)
    \end{aligned}
\end{eqnarray}
which is the free energy in Eq. \ref{eq:tPC-FE} with a linear $f$. The above derivation thus provides another interpretation of the tPC model, i.e. a special case of Bayesian filtering that assumes observation and state noise are the same—and known. Again, as all transformations in the optimization objective are linear, we can derive an analytical expression for $\hat{\mathbf{x}}_k$:
\begin{eqnarray}\label{eq:PC-update}
\begin{aligned}
    \hat{\mathbf{x}}^-_k &= \bfW \hat{\mathbf{x}}_{k-1} + \bfB \mathbf{u}_{k} \\
    \hat{\mathbf{x}}_k &= \hat{\mathbf{x}}_k^- + \bfK(\mathbf{y}_k - \bfF \hat{\mathbf{x}}_k^-) \\
    \bfK &= \bm{\Sigma}_x \bfF^{\top} \big[ \bfF \bm{\Sigma}_x \bfF^{\top} + \bm{\Sigma}_y]^{-1}
\end{aligned}
\end{eqnarray}
which is similar to the projection and correction steps for Kalman filtering (Eqs. \ref{eq:KF-projection} and \ref{eq:KF-correction}), but with $\bm{\Sigma}_{k-1}$ assumed to be $0$, i.e., tPC does not propagate the uncertainty estimations. The derivation of these equations can be found in Appendix \ref{app:chap4}.

It is also worth mentioning that in the tPC model, although we did not specify the estimation of uncertainty $\bm{\Sigma}_k$ at each step, not \emph{propagating} the uncertainty $\bm{\Sigma}_{k-1}$ from the previous step is not equivalent to not \emph{estimating} it. In fact, as was shown in \citep{baltieri2021kalman}, even when we choose to estimate $\bm{\Sigma}_k$ in the tPC model, it is still not propagated. Therefore, our MAP estimation $\hat{\mathbf{x}}_k$ will not be affected.

To summarize, by assuming a linear $f$, here we have shown that both Kalman filtering and tPC are special cases of the Bayesian filtering problem, while the key difference is that tPC ignores the uncertainty estimated at the previous step when estimating the most likely hidden cause at the current time step, whereas Kalman filtering always estimates this uncertainty. However, as we will show in the experimental results section, in the benchmark tracking tasks, tPC performs on par with Kalman filtering, albeit not estimating the posterior uncertainty. Importantly, there are several advantages of tPC over Kalman filtering as a model of dynamical processing in the brain. Firstly, the projection and correction steps of Kalman filtering require complicated matrix algebra and are challenging to compute in neural circuitry, especially the Kalman Gain matrix $\bfK$. On the other hand, although we can derive analytical results for the estimates of tPC as well, these estimates can be obtained via the iterative methods mentioned above, which afford plausible circuitry implementations (Fig \ref{fig:recurrent_multistep}). Secondly, the iterative nature of our tPC model also makes it adaptable to nonlinear $f$, where there are no analytical solutions to the Bayesian filtering problem. In contrast, extending the Kalman Filter to nonlinear systems is challenging and standard methods such as the extended Kalman filter \citep{ruck1992comparative} work by linearizing around the nonlinearity and thus require knowledge of the Jacobian of the nonlinearity at every state, which is also challenging to implement in neural circuits. Finally, the Kalman filter assumes knowledge of the correct $\bfW$, $\bfB$, and $\bfF$ matrices while these must presumably be learned from sensory observations in the brain. 

\subsection{Performance in linear filtering problems}
Here we first present results for the linear tPC model on a simple tracking task of the kind to which the Kalman filter is commonly applied in industry.
% , and which also shares a close relationship with those the brain must solve -- such as tracking a moving object with visual saccades during smooth pursuit. 
Here, the goal is simply to infer the unknown hidden state (position, velocity, acceleration) of an object that is undergoing an unknown acceleration. We receive noisy observations of the position, state, and acceleration of the object which are mapped through a random $\bfF$ matrix with additional observation noise. We use a random $\bfF$ matrix for the observation mapping to simulate and test the most difficult scenario where the observations are entirely scrambled. 

Mathematically, the generative process of this task can be represented according to a linear state space model. We assume that the position and velocity of the object follow the usual laws of Newtonian physics, and there is a persistent acceleration which is affected by the control input, giving us the following $\bfW$ and $\bfB$ matrices,
\begin{eqnarray}
    \label{eq:A}
    \bfW = \begin{bmatrix}
    1 & \Delta k & \frac{1}{2}\Delta k^2 \\
    0 & 1 & \Delta k \\
    0 & 0 & 1
    \end{bmatrix}, \
    \bfB = \begin{bmatrix}
    0 & 0 & 1
    \end{bmatrix}.
\end{eqnarray}

In Eq. \ref{eq:A}, $\Delta k$ denotes the duration of the interval between successive sensory observations (we used $\Delta k=0.001$).
Additionally, we draw a fixed $\bfF$ matrix from a random Gaussian distribution $\bfF \sim \mathcal{N}(0, 1)$, {and modeled the control inputs as $\mathbf{u}_k=e^{-0.01k}$}. The process and observation noise $\bm{\Sigma}_x$ and $\bm{\Sigma}_y$ were set to identity matrices. We then generate the true latent states $\mathbf{x}_k$ and the noisy observations $\mathbf{y}_k$ using Eqs. \ref{eq:general_transition} and \ref{eq:general_emission}, initialized with a zero vector when $k=0$. The performance of the models is then measured as the mean squared error (MSE) between the estimated $\hat{\mathbf{x}}_k$ and true $\mathbf{x}_k$ across all observation time steps. Fig \ref{tracking_performance_fig}A shows an example of the true system state that we generated across $1000$ time steps, and Fig \ref{tracking_performance_fig}B shows its corresponding noisy observations. As can be seen, the projected observations are completely scrambled by matrix $\bfF$,
%and $\bm{\Sigma}_y$, 
making it a challenging task for the models to retrieve the true system states.

\begin{figure}[t]
\centering
    \includegraphics[width=0.95\linewidth]{figs/chap4/Fig3.pdf}
    \caption{{\bf The tracking task and the impact of inference step size and the number of inference steps on performance.} A. The dynamics of the true hidden state are represented as a $3$-dimensional vector at each time step, with entries corresponding to position ($x_1$), velocity ($x_2$) and acceleration ($x_3$). B. The projected noisy observations from the true system state in A. C: Estimates of the acceleration with different models, zoomed in at the interval between $560$ and $600$ time steps. D: MSE difference between tPC and Kalman filter, with varying numbers of inference steps and step sizes for predictive coding. PC stands for temporal predictive coding and KF stands for Kalman filter. {All values are with arbitrary units (a.u.).}}
    \label{tracking_performance_fig}
\end{figure}

We then investigate tracking performance using tPC compared to the Kalman filter for both $5$ steps of inference between observations ($\Delta k = 5 \Delta t$) and $1$ step ($\Delta k = \Delta t$). Since the problem is linear, we also investigate the performance of a tPC when its inference dynamics have reached the equilibrium, using the equilibrium condition derived in Eq. \ref{eq:PC-update}. Since tracking performance is {visually indistinguishable} when zoomed out over $1000$ timesteps, in Fig \ref{tracking_performance_fig}C, we plot the estimates of acceleration ($x_3$) on the $40$ timesteps between $560$ and $600$ steps in. It can be seen from Fig \ref{tracking_performance_fig}C that both the tPC model with $5$ inference steps and that with fully converged inferential dynamics could achieve comparable performance to the Kalman filter. 
Interestingly, the estimates of tPC tend to be closer to those of the Kalman filtering, rather than the true values.
Although the tPC model with a single inference step (corresponding to the neural implementation in Fig \ref{fig:recurrent_multistep}C) has worse tracking performance, it is able to estimate a smoothed version of the trajectory of the system state. 
% This allows us to see more clearly how the predictive coding and Kalman filtering estimates are performing. 
% Interestingly, we see that overall the predictive coding estimate appears largely to be less variable than the Kalman filter estimate and is effectively a smoothed version of it. 
We hypothesize that the smoothed estimate with a single inference step is likely due to the fact that the tPC model does not completely converge in $1$ iteration, and so does not completely optimize its estimate on every timestep, with the effect that the estimate is less sensitive to new information and effectively averages over recent experiences rather than optimally solving each one independently. 
% We also see that both the $5$ and the $1$ gradient step largely track the underlying true state correctly, although a single gradient step performs slightly worse and, in the case of the acceleration diverges somewhat from the true signal.
A similar performance comparison is obtained on the position ($x_1$) and velocity ($x_2$) and is shown in Appendix \ref{app:chap4}.

To quantify the effect of the number of inference iterations and integration step size $\frac{\Delta t}{\tau}$ upon performance, in Fig \ref{tracking_performance_fig}D we plotted the MSE difference between predictive coding models with various inference iterations and inference step sizes and the Kalman filter, which is the optimal solution to the tracking problem. The MSE is calculated as the mean squared difference between the estimated system state $\hat{\bfx}_k$ and true state $\mathbf{x}_k$, averaged across time steps and trials. We find that with a small number of inference steps, the performance of the tPC model is worse, indicating that additional steps of inference aid the tracking
performance of the algorithm. Moreover, although increasing the step size will initially improve the performance, the MSE will start to increase if the step size is too large. 
%We believe that the relatively high step sizes allow the recurrent predictive coding network to converge quickly, but can cause overshooting and oscillatory behaviour once the network is close to the optimum. However, due to the strict real-time constraints of computation in the brain, we believe that it is likely that the brain uses large learning rates to quickly get usable inferences rather than small step sizes to obtain very refined inferences slowly. 
It can also be seen that with more inference steps and appropriate step sizes, the performance of tPC will be able to approximate that of the (optimal) Kalman filter.

\subsection{Learning the synaptic weights}

In the previous investigations, we fixed the parameters $\bfW$, $\bfB$ and $\bfF$ to the true values and only performed the inference dynamics. However, in many cases, simply inferring the hidden states of the world is not enough because we cannot assume that we know a-priori the structure of the dynamics or observation functions of the world. That is, in most real-world situations, the $\bfW$, $\bfB$, and $\bfF$ matrices are unknown. Instead, we must \emph{learn} these matrices from observations. In the tPC model, we have a natural Hebbian plasticity-based learning rule which we can use to learn these matrices directly (Eq. \ref{eq:A_C_updates}). Here, we investigate how learning these parameters affects the performance of the tPC model. Specifically, we compare three different ways of setting the values for $\bfW$ and $\bfF$: 1) fixing them to the true values used for generating the data; 2) learning them using Eq. \ref{eq:A_C_updates} and Algorithm \ref{alg:tPC}; 3) fixing them to random values. We then examine the performance of these models on two levels, the latent state level ($\mathbf{x}$) and the observation level ($\mathbf{y}$), by measuring how well the model estimates the activities on both levels. It is worth noting that the observation estimates are calculated by performing a forward pass at each time step i.e., $\hat{\mathbf{y}}_k = \bfF(\bfW\hat{\mathbf{x}}_{k-1}+\bfB\mathbf{u}_k)$, where the values of $\bfW$ and $\bfF$ are obtained at each time step $k$ via the aforementioned three approaches. The results are shown in Fig \ref{linear_learning_fig}. 

\begin{figure}[t]
    \centering
    \includegraphics[width=0.9\linewidth]{figs/chap4/Fig4.pdf}
    \caption{{\bf Effects of learning parameters $\bfW$ and $\bfF$.} {A, B}: Estimation of the state and observation trajectories respectively by different models. `True', `learned' and `Random' denote the predictive coding model with true, learned and random $\bfW$ and $\bfF$ respectively. Only the first dimension of the latent and observation is shown for simplicity. The other two dimensions have similar performance. {C, D}: MSE of the predictions on the hidden and observation levels respectively. Boxplots were obtained with $40$ trials for each model. {Both $x$ and $y$ are with arbitrary units (a.u.).}} 
    \label{linear_learning_fig}
\end{figure}

For this set of results, we use $20$ inference steps with a step size $0.2$ to get the optimal performance for the tPC models, based on Fig \ref{tracking_performance_fig}D. Fig \ref{linear_learning_fig}A shows that, while our tPC estimates the latent state well with true $\bfW$ and $\bfF$, when asked to learn the parameters, the model fails to accurately estimate the latent state. 
Quantitatively, as shown in Fig \ref{linear_learning_fig}C, the estimation MSE of the learning model on the state level is similar to tPC with totally random parameters, which is much higher than those of the Kalman filter and tPC with true parameters. On the other hand, however, we find that the model learning $\bfW$ and $\bfF$ can accurately estimate the observations even with the incorrectly estimated latent state (Fig \ref{linear_learning_fig}B). This effect arises because the problem of inferring the true hidden state from the data is fundamentally under-determined. There are many possible hidden states that, given a flexible learned mapping, could result in an identical predicted observation. %As such, absent any regularization, the 
Importantly, despite inferring a different representation of the hidden state, the network is able to learn $\bfW$ and $\bfF$ that correctly predicts the incoming observations (Fig \ref{linear_learning_fig}D).

\subsection{Learning the noise covariance matrices}
For all our experiments above, we have used $\bm{\Sigma}_x = \bm{\Sigma}_y = \bfI$ when generating the training data, where $\bfI$ is the identity matrix. However, the noise covariance underlying a natural dynamic process may not always be identity. Although earlier works have proposed to encode the noise precision matrices $\bm{\Sigma}_x^{-1}$ and $\bm{\Sigma}_y^{-1}$ into additional connections explicitly \citet{bogacz2017tutorial}, this approach would introduce extra complexity into the neural implementation of tPC.
%(Fig~\ref{fig:recurrent_multistep}). 
On the other hand, it has been shown that in the static case, the noise precision matrix can be implicitly encoded in recurrent connections similar to the $\bfW$ matrix in our tPC model \citep{tang2023recurrent}, without needing to represent the precision matrix explicitly as in \citep{bogacz2017tutorial}. Therefore, here we investigate whether $\bfW$ and $\bfF$ can encode the precision matrices of the process and observation noise respectively after learning. We used the same $\bfW$, $\bfB$ and $\bfF$ matrices for data generation as before, but set the noise covariance matrices as:
\begin{eqnarray}
    \bm{\Sigma}_x = \bm{\Sigma}_y = \begin{bmatrix}
        10 & 0 & 0 \\
        0 & 1 & 0 \\
        0 & 0 & 1 \\
    \end{bmatrix} \text{or } \bm{\Sigma}_x = \bm{\Sigma}_y = \begin{bmatrix}
        10 & 2 & 0.5 \\
        2 & 1 & 0.4 \\
        0.5 & 0.4 & 1 \\
    \end{bmatrix}.
\end{eqnarray}

We refer to the first case as ``non-identity diagonal'' and the second case as ``positive definite''. The choice of these covariance matrices is arbitrary, although we intentionally make one diagonal entry larger than the others to examine tPC's capability of learning such a unique structure. The training hyper-parameters for this experiment are identical to the ones used for Fig~\ref{linear_learning_fig}, with which the learning for $\bfW$ and $\bfF$ could also converge. As Fig~\ref{fig:learning-precision}A shows, when the noise covariance matrices are non-identity diagonal and positive definite, the recurrent matrix $\bfW$ is able to encode the large diagonal entry $10$ in $\bm{\Sigma}_x$ into its diagonal entries. Likewise, the $\bfF$ matrix develops stronger weights to account for the larger variance of the inputs, suggesting that our tPC model can learn the noise covariance $\bm{\Sigma}_x$ and $\bm{\Sigma}_y$ into its recurrent and feedforward weights, without accounting for them in the energy function and learning/inference dynamics explicitly. 
%It is also worth noting that the learned weights are also wrong in these cases, further demonstrating the under-determined nature of the model.

\begin{figure}[t]
    \centering
    \includegraphics[width=\linewidth]{figs/chap4/Fig5.pdf}
    \caption{{\bf Performance with non-identity noise covariance.} A: True and learned $\bfW$ and $\bfF$ matrices with different underlying noise covariance matrices. B, C: MSE of the predictions on the hidden and observation levels with different noise covariance matrices. Error bars obtained with $40$ trials.} 
    
    \label{fig:learning-precision}
\end{figure}

We then conducted a quantitative analysis of the impact of non-identity noise covariance in Fig~\ref{fig:learning-precision}B and \ref{fig:learning-precision}C, similar to Fig~\ref{linear_learning_fig}. Here, both ``Kalman Filter'' and ``True'' use the correct $\bfW$, $\bfF$, $\bm{\Sigma}_x$ and $\bm{\Sigma}_y$, although the ``True'' tPC model performs inference to get the hidden states using Eqs.~\ref{eq:err-y},\ref{eq:err-x} and \ref{eq:dynamics_x}. The ``learned'' model has no access to all the matrices and has to learn both the generative process and the noise covariance with its weights. The results are similar to Fig~\ref{linear_learning_fig}: the model that learns $\bfW$ and $\bfF$ will fail to learn the correct hidden states but its observation estimates are on par with Kalman filtering, even when the noise covariance matrices are non-identity. Interestingly, we also observed slightly degraded latent estimation performance of the ``True'' model with non-identity noise covariance. We hypothesize that this is due to the large diagonal entry $10$ affecting the inference dynamics, which can be solved using coordinate-wise inference step sizes.

Overall, these results suggest that the synaptic plasticity of tPC can encode the noise covariance in the generative process, without representing them explicitly in the dynamics and circuit implementations. It is also interesting to investigate the exact theoretical relationship between the weights and noise covariance, similar to the analytical relationship in \citet{tang2023recurrent}, and we intend to investigate it in future explorations.

\subsection{Training tPC with natural movies}

Thus far, we have examined tPC in low-dimensional examples to understand its computational properties.
In this section, we demonstrate that this model can also be applied to high-dimensional stimuli and provide a plausible account of how the biological visual system develops representations of dynamical inputs. To do so, we trained the tPC model on a dataset of patches extracted from movies of natural scenes. Each of the movies consisted of 50 frames of $200 \times 200$ pixels that we spatially bandpass filtered by a retina-like centre-surround filter \citep{olshausen1997sparse}. We extracted $16 \times 16$ patches from the movies with the highest motion to form a dataset of $2000$ moving patches of $50$ frames. We also augmented the dataset with a left-right flipped version of itself to ensure a rich variety of motion directions. Two examples of the first $10$ frames of the movies can be seen in Fig~\ref{fig:natural_movie}A. 
We then trained a tPC model with $320$ hidden neurons with this dataset. Specifically, due to the redundancy of information in natural scenes and neural connectivity and energy constraints, we trained a tPC model with $l1$ sparsity-constraint latent activities and forward weights $\bfF$. Therefore, the objective function we used in this set of experiments is:
\begin{eqnarray}\label{eq:sparse}
    \mathcal{F}_k = \Vert \mathbf{y}_k - \bfF f(\mathbf{x}_k) \Vert_2^2 + \Vert \mathbf{x}_k - \bfW f(\hat{\mathbf{x}}_{k-1}) \Vert_2^2 + \lambda_x |\mathbf{x}_k|_1 + \lambda_{F} \sum_{i,j} |F_{ij}|
\end{eqnarray}
where $|\cdot|_1$ denotes vector $1$-norm, and we used a linear $f(\cdot)$ in these experiments. We also assume there is no control input to the model, so $\bfB=0$. The sparsity constraints are similar to the classical sparse coding model \citep{olshausen1997sparse} and recent temporal prediction neural network models of dynamical inputs \citep{singer2018sensory, singer2023hierarchical}. To describe the cortical processing of natural scenes more accurately, for the experiments with these natural movies, we also used a time constant $\tau=10ms$ estimated for visual cortical neurons \citep{nowak2003electrophysiological}. 
% That is, the discretized inferential dynamics in Eq.~\ref{eq:inf-discrete} becomes:
% \begin{eqnarray}\label{eq:dynamic-tau}
%     x_k(t+\Delta t) = x_k(t) + \frac{\Delta t}{\tau} \left[ - \epsilon_x  + f'(x_k(t))\odot C^T \epsilon_y \right]
% \end{eqnarray}
% where we used $\tau=10ms$ according to experimental data from visual cortex \citep{nowak2003electrophysiological}.
The natural movies used in this experiment were recorded with frame rates in the range of $20$Hz to $30$Hz. In our simulations, we set $\Delta t$ to $0.1ms$ (and therefore the inference step size $\frac{\Delta t}{\tau}=0.01$) and perform inference for $330$ iterations to (approximately) match the $30$Hz frame rate, although we did not observe any qualitative difference when we used more inference iterations to match the lower frame rates.

We then trained our adjusted model with the natural movies. In Fig~\ref{fig:natural_movie}B we show the (sparsity-constraint) forward weight $\bfF$ (of size $256 \times 320$), which develops Gabor-like and localized projective fields similar to those observed in mammalian primary visual cortex \citep{hubel1962receptive}. This finding is unsurprising as such a sparsity-constraint generative model is similar to the original predictive coding \citep{rao1999predictive} and sparse coding models \citep{olshausen1997sparse}, which also reproduced localized Gabor filters. 

\begin{figure}[t]
    \centering
    \includegraphics[width=0.9\linewidth]{figs/chap4/Fig7.pdf}
    \caption{{\bf Representations developed by the model when trained with patches from movies of dynamic natural scenes.} A: First $10$ frames of $2$ example training movies used in our experiments. Patches extracted from movies obtained at websites pexels.com, pixabay.com and commons.wikimedia.org. B: The projective fields $\bfF$ developed Gabor-like filters after training. C: Space-time receptive fields developed by hidden neurons of the tPC model.}
    \label{fig:natural_movie}
\end{figure}

To study the dynamical properties of tPC, we then performed a reverse correlation analysis \citep{dayan2005theoretical} on the hidden neurons of the trained tPC. 
Specifically, after training, we supply a sequence of white noise stimuli to the model, let the hidden neurons relax according to Eq.~\ref{eq:inf-discrete} to develop latent representations of the stimuli, and average the stimuli giving more weights to those producing the highest response of a neuron. In particular, for each time step $k$ and each neuron, we multiply the hidden activity $\mathbf{x}_k$ with each of the white noise stimuli from $\mathbf{y}_{k-4}$ to $\mathbf{y}_k$ and sum up the products. The weighted sum is the spatio-temporal receptive field (STRF) of this neuron \citep{dayan2005theoretical}. Formally, the STRF of the $i$th hidden neuron is defined as:
\begin{eqnarray}
    \text{STRF}_i = \frac{1}{T-5} \sum_{k=5}^T x_{i,k} y_{k-4:k},
\end{eqnarray}
where $x_{i,k}$ denotes and activity of the hidden neuron $i$ at time step $k$ and $y_{k-4:k}$ denotes the $5$ frames of white noise preceding and including time step $k$. Some examples of the STRFs developed by the hidden neurons of the tPC are shown in Fig~\ref{fig:natural_movie}C, where each row denotes a neuron and each column a time step. Two important properties emerge: 1) The STRFs are temporally asymmetric i.e., the power of the receptive fields decays back in time. This is consistent with observations in real neurons \citep{ohzawa1996encoding} and earlier computational models \citep{singer2018sensory, singer2023hierarchical}; 2) Importantly, many neurons also develop motion-sensitive STRFs, i.e., the spatial regions each neuron is responsive to shift location in the visual field over history steps, moving in a fixed direction (right-ward for neurons $160$ and $298$, left-ward for neurons $14$ and $32$ in the shown examples). Such motion sensitivity has been observed in real neurons in early visual areas \citep{baker1990spatial,livingstone1987psychophysical} and demonstrated in earlier sparse coding models \citep{bogacz2000emergence,olshausen2002sparse}. 
%However, these previous models were not really `temporal'; instead, they treated time as a spatial dimension. On the other hand, our tPC model is fully temporal and dynamically processes sequential inputs. 
Overall, these results demonstrate that the tPC model with realistic time constants matching cortical neurons can reproduce neural representations observed in the visual areas, providing a possible computational mechanism underlying the learning of such representations.
% \end{revised}

% \subsection{Extending tPC to nonlinear tasks}

% {Here we examine our tPC model in nonlinear tasks. In order to make the task more challenging, we train the model on a simulated motion of a pendulum, in} which the generative process does not explicitly follow that of tPC.
% Fig \ref{pendulum-full}A shows a free-body diagram of the pendulum that we simulated in this experiment, demonstrating the mass, length $L$, and force vectors acting on the system.
% We describe the state of the pendulum by the angle of the pendulum $\theta_1$ and its angular velocity $\theta_2$.
% According to Newton's Second Law of Motion, the angle of the pendulum evolves according to:  
% % system using the rotation motion form:
% % \begin{align}
% %        \label{rotational_motion}
% %        \tau = I \alpha
% % \end{align}
% % which is just a slightly modified version of Newton's second Law of Motion ($F = ma$). Here $\tau$ is the torque and is given by the perpendicular component force on the mass, $I$ is the rotational inertia and is equal to $mL^2$ where $m$ is the mass hanging from a string with length $L$, and $\alpha$ is the angular acceleration given by the second time derivative of $\theta$, the angle between the gravity and force vectors. Thus, Eq. \ref{rotational_motion} can be re-written as:
% \begin{eqnarray}
%     \begin{aligned}
%        \label{rotational_motion_2}
%        \ddot{\theta}_1 = -\frac{g}{L}\sin(\theta_1) 
%        %&= mL^2 \ddot{\theta}\\
%       %^\Rightarrow mL^2 \ddot{\theta} + -%mgL\sin(\theta) &= 0 \\
%       %\label{rotational_motion_friction} 
%     %   \Rightarrow mL^2 \ddot{\theta} + bL^2 \dot{\theta} + -mgL\sin(\theta) &= 0
% \end{aligned}
% \end{eqnarray}
% where $g$ is gravity. 
% % Note that we did not include a damping parameter in the equations above as we will not use any damping in our simulations. The second-order ODE described here was reduced to two first-order equations by introducing $\theta_1$ and $\theta_2$ to represent the angular displacement and angular velocity, respectively, as:
% % \begin{align}
% %        \label{angular_displacement}
% %        \theta_1 &= \theta\\
% %       \label{angular_velocity}
% %        \theta_1 &= \dot{\theta}
% % \end{align}
% We can express this motion in a set of first-order equations:
% \begin{eqnarray}
%     \begin{aligned}
%        \label{first_order_1}
%        \dot{\theta}_1 &= \theta_2\\
%       \label{angular_velocity}
%        \dot{\theta}_2 &= - \frac{g}{L} \sin(\theta_1).
% \end{aligned}
% \end{eqnarray}

% \begin{figure}
%     \centering
%     \includegraphics[width=0.9\linewidth]{figs/chap4/figure6.pdf}
%     \caption{{\bf Simulations of the pendulum.} A: A free-body diagram of a simple pendulum that has a mass $m$ attached to a string with length $L$. Also shown are the forces applied to the mass. The restoring force $-mg \sin \theta$ is a net force toward the equilibrium position. B: A phase portrait of the pendulum simulation showing the result of our linear versus nonlinear models prediction for the ground-truth data. The vector field (i.e. set of small arrows) was created by computing the derivatives of $\frac{d \theta_1}{dt}$ and $\frac{d \theta_2}{dt}$ at $t = 0$ on a grid of 30 points over the range of $-\pi$ to $+\pi$ and -4 to +4 for $\theta_1$ and $\theta_2$, respectively. C: The barplot shows the difference between the mean prediction errors of the linear model versus the nonlinear model from 100 simulations with varying noise profiles. The mean errors are significantly different ($p << 0.001$).}
%     \label{pendulum-full}
% \end{figure}

% In the simulation, we used the following parameters: $g$ = 9.81 m/s\textsuperscript{2}, $L$ = 3.0 m.
% %, and $m$ = 0.1 kg. 
% We then simulated the system as an initial value problem by numerically integrating the equations using the explicit Runge-Kutta method for 2500 seconds with $\Delta k = 0.1$ second time steps and initial values of $\theta_1$ = 1.8 rad and $\theta_2$ = 2.2 rad/s. These values were chosen to simulate the pendulum motion with a large amplitude of oscillation to shift the system into a more nonlinear regime.
% The time series presented to the models $y_k$ were created by adding these numerical solutions and Gaussian noise with a zero mean and standard deviation of 0.1. 

% We trained both nonlinear (with a hyperbolic tangent nonlinearity $f(\cdot)$) and linear tPC models to predict these time series. Parameters of $\bfF$ and $\bfW$ were initialized {to the identity matrix and the zero matrix}, respectively. The learning was performed {in the same fashion} as in the previous section.

% Fig \ref{pendulum-full}B shows the results from the pendulum simulations using the phase portrait of the system.  The solutions of the ground-truth simulation and our nonlinear model prediction are plotted on the vector field for the final 80 seconds. Even though both models performed relatively well by correctly predicting the behaviour of the pendulum motion, the linear model performed worse when the pendulum reached its highest angular displacements (see the noisy prediction around the extremities of the $\theta_1$ axis in Fig \ref{pendulum-full}B). On the other hand, the nonlinear model stayed relatively stable throughout the entire prediction. 
% Fig \ref{pendulum-full}C shows the mean squared error of prediction averaged over 100 separate simulations, and the nonlinear model consistently outperformed the linear one. 
% In Appendix \ref{app:chap4} we also examined a nonlinear task where the data generation process follows the generative model in Eqs.~\ref{eq:general_emission} and \ref{eq:general_transition} and obtained similar results.

\section{Discussion}

\subsection{Summary}
In this chapter, we have analysed the recurrent predictive coding architecture for temporal prediction and filtering. This task is important because processing time-varying sequences of inputs and using them to infer dynamically changing hidden states of the world is often considered a core task of the sensory regions of the brain \citep{doya2007bayesian}.
%, much more so than simple Bayesian inference on static datasets. 
As such, it is likely that these regions have an architecture heavily specialised for performing such filtering tasks. Here we have shown that the filtering problem can be tackled with a simple and biologically plausible algorithm with a straightforward implementation in neural circuitry.

We have derived our tPC model from first principles as a variational filtering algorithm, providing a clear algorithmic derivation in terms of gradient descents on the resulting free energy functionals. We have also proposed a direct implementation in neural circuitry which relies on only local information transmission as well as purely Hebbian plasticity that inherits the simple neural implementation of static predictive coding networks. Furthermore, we have also demonstrated that in the linear case, the algorithm is closely related to Kalman filtering, and is capable of robustly solving filtering and tracking tasks at a level close to the optimal linear solution. Moreover, and unlike the Kalman filter, we have demonstrated that our model can perform online \emph{learning} of the parameters using Hebbian plasticity, which works rapidly and effectively in predicting the correct observations. {Importantly, when trained on natural stimuli and constrained by sparse weights and activities, the tPC model develops motion-sensitive Gabor-filters of the visual scene, which is consistent with representations developed in the visual cortex.}

\subsection{Neural implementation of time}
There are several interesting questions regarding the biological plausibility of the multi-step neural implementation in tPC. Initially these schemes, while they arise directly from the gradient descent derivation, appear biologically implausible for two reasons. The first is the issue of storage. Iterative schemes require the initial conditions (state estimate $\hat{\bfx}_{k-1}$) to be held fixed throughout multiple iterations, and this means that this state must be stored somewhere accessible to be utilized multiple times during the iterative inference phase. It is not clear where or how this information could be stored in the brain, especially in low-level sensory systems. This storage must be local and ubiquitous as a naive implementation of the multi-step algorithms proposed in this work would require separate storage for every single value neuron. 
%Interestingly, one way this problem may be avoided for the multi-step algorithm is that instead of computing the prediction errors between the newly computed value and the originally stored value, we instead just overwrite the storage at every step, and compute the prediction errors between the newly computed value and the previously computed value, which can be done naturally with recurrent connections with synaptic delays and which thus requires no explicit storage. %While this algorithm diverges from the mathematical derivation of predictive coding as variational inference, we have shown empirically that it still enables robust tracking performance in practice. 

The second issue relates to the time it takes for an iterative algorithm to converge. Specifically, if we model the brain as receiving visual input as a continuous stream, then multiple iterations based upon a single stimulus would necessitate ignoring the data arriving in the intervals while the iterations are taking place. Moreover, an iterative approach would also take more time to update upon information newly received, which could be crucial for survival in some cases. 

In the modeling results presented in this thesis, these two issues are addressed using the single-step tPC algorithm (\ref{fig:recurrent_multistep}C). However, the performance of this algorithm could not perfectly match that of the iterative tPC algorithm. To achieve superior performance in temporal processing tasks, it is therefore likely that the brain may still employ some strategies to achieve a more ``complete'' inference. One potential solution to this issue is that the cortex may implement both an iterative and an amortized feedforward pass solution simultaneously \citep{tscshantz2023hybrid}, similar to what has been demonstrated in machine learning \citep{kingma2013auto}. There is also evidence for precisely this: Core visual functions can often occur within 100-200 ms \citep{thorpe1996speed,keysers2001speed,carlson2013representational,thunell2019memory} which is too short to allow multiple steps of iterative optimization. Alternatively, neurons with different inferential speed may serve different purposes: In lower levels of the cortex, neurons operate at a lower frequency to process more refined features, whereas in higher levels such as the hippocampus, neurons perform faster inference, which is sufficient to form some high-level representation of stimuli. This is reflected in our Fig.~\ref{tracking_performance_fig}C, where the fitted curve of the single-step algorithm is essentially a smoothed version of the True curve. Biologically, this hierarchically distributed inference speed may be achieved through different neural oscillation frequencies \citep{buzsaki2004neuronal,buzsaki2006rhythms}. In particular, the fast inference can be achieved through hippocampal sharp wave-ripples, which serves as a compression of sensory information \citep{buzsaki2015hippocampal}.

\subsection{Other temporal predictive coding models}
Before our proposal of tPCN, some earlier studies have also explored predictive coding for temporal predictions. For instance, early studies \citep{rao1997correlates,rao1997dynamic,rao1999optimal} utilized a Kalman filter combined with sparse image representations to make future predictions of visual stimuli. However, these works did not describe how the Kalman filter can be implemented in biological circuits, and how their Kalman filter-based models can be extended to the nonlinear case.
More recently, \citet{jiang2022dynamic} trained a temporal version of predictive coding on sequences of natural video (filmed by a person walking through a forest) and observed that neurons have spatiotemporal receptive fields resembling those in the primary visual cortex. {However, unlike our model, their model relies on an external hyper-network to perform temporal predictions.}
Temporal versions of predictive coding networks have also been extended to include multiple levels of hierarchy \citep{jiang2022dynamic, ororbia2020continual}. 
Analysis of these networks revealed that neurons on higher levels change their activity with a slower time scale than the neurons at the lower levels of the hierarchy \citep{jiang2022dynamic}. 
It has been also demonstrated that hierarchical temporal predictive coding networks can achieve performance comparable to BPTT in standard machine learning benchmarks \citep{ororbia2020continual}. However, these models require complex neural network implementations to perform these temporal tasks. It is thus an interesting future direction to see whether our tPC model, which inherits the simple neural implementation of the static predictive coding network, can present similar performance and neural responses.

A number of studies have used normative models to generate spatiotemporal receptive fields resembling those of direction-selective V1 simple cells \citep{bogacz2000emergence, olshausen2002sparse, van1998independent, pachitariu2012learning,pachitariu2017visual,berkes2005slow,hyvarinen2003bubbles, palm2012prediction,singer2018sensory,singer2023hierarchical}. Some models involved mechanisms related to predictive coding, such as sparse coding \citep{bogacz2000emergence, olshausen2002sparse} and independent component analysis \citep{van1998independent} and applied them to spatiotemporal stimuli. However, 
%these models treated time like a spatial dimension and hence lacked realistic dynamics. Also, 
while the resulting receptive fields were sometimes direction-selective, they did not have the asymmetric temporal profile seen in real STRFs. Application of slowness principles to spatiotemporal stimuli can also produce direction-selective STRFs \citep{berkes2005slow,hyvarinen2003bubbles}, but they again lack the asymmetric temporal profile.  Other models have trained single-hidden-layer neural networks to perform temporal prediction on movies of natural scenes \citep{palm2012prediction,singer2018sensory,singer2023hierarchical}. This has been shown to reproduce direction-selective STRFs with an appropriate asymmetric temporal profile \citep{singer2018sensory}, and indeed when stacked hierarchically reproduces units resembling motion-sensitive simple, complex and pattern-motion cells \citep{singer2023hierarchical}. However, these models were trained by backpropagation, hence while they are informative about whether the cortex might have temporal prediction as a normative objective, they are agnostic about the potential learning mechanism. This thesis helps bridge this gap, demonstrating how temporal prediction can be combined with a more plausible learning mechanism to produce units with spectrotemporal receptive fields resembling those of direction-selective V1 simple cells.

\subsection{Kalman filtering}
Several earlier works have tried to approach the problem of Kalman filtering in the brain.  \citet{wilson2009neural} repurpose a line attractor network and show that it recapitulates the dynamics of a Kalman filter in the regime of low prediction error. However their model only works for a single-dimensional stimulus, does not encode uncertainty, and also only works when a linearisation around zero prediction error holds. \citet{deneve2007optimal} encoded Kalman filter dynamics in a recurrent attractor network. Their approach however encodes stimuli by means of basis functions, which leads to an exponentially growing number of basis functions required to tile the space as the dimensionality of the input grows. In the predictive coding approach, neurons directly encode the mean of the estimated posterior distribution, which means that the network size scales linearly with the number of dimensions. Our gradient method also completely eschews the direct computation of the Kalman gain, which simplifies the required computations significantly. Additionally, \citet{beck2007probabilistic} show that probabilistic population coding approaches can compute posteriors for the exponential family of distributions of which the Gaussian distribution is a member. However, no explicitly worked-out application of the population coding approach to Kalman filtering exists, to our knowledge. \citet{rao1997dynamic} have used Kalman filtering as a model of dynamical processing in the brain, {and recent works have also been interested in Kalman filtering in a biologically plausible setting \citep{friedrich2021neural}}. We have found that the crucial distinction between tPCN and Kalman filtering is in the representation of the model's uncertainty, mathematically represented in the posterior variance $\bm{\Sigma}_k$. Specifically, in the tPC model, although it represents the two ``objective'' uncertainties $\bm{\Sigma}_x$ and $\bm{\Sigma}_y$ of the dynamical system it is inferring, it does not represent the uncertainty in the estimated posterior distribution (unlike the extended but inefficient model of \citet{kutschireiter2017nonlinear}, which uses multiple copies of the network to create samples from the posterior distribution). Crucially, it is the assumption of the prior at each time step that prevents the tPC model from successfully propagating the posterior uncertainty through time. 

Although tPCN does not optimally update and represent the dynamically changing posterior uncertainty of the agent, it has some computational advantages over Kalman filtering, which may render it more suitable for implementation in neural circuitry. The key advantage of the predictive coding approach is its computational simplicity compared to the Kalman filtering update rules which require complex matrix algebra and especially matrix inversions to compute the Kalman Gain matrix which are unlikely to be implementable in neural circuitry directly, while the predictive coding equations are simple and only require local and Hebbian updates and can be directly translated into relatively straightforward neural circuits. Moreover, as seen in Fig \ref{tracking_performance_fig}, the predictive coding estimate and the Kalman filter estimates of a given dynamical system end up closely converging anyway, which means that tPC networks could provide the brain an efficient and cheap way to approximate the highly effective Kalman filter using only simple circuitry.

One reason why predictive coding networks achieve performance similar to the Kalman filter is that the posterior uncertainty decays rapidly over trials (for an illustration see Fig 5A in \citet{moeller2022uncertainty}). Furthermore, for deterministic transition processes (with $\bm{\Sigma}_x=0$), the posterior uncertainty decays to $\bm{\Sigma}_k=0$ \citep{moeller2022uncertainty}. Therefore, as the learning progresses, the Kalman filter becomes more similar or even identical (for deterministic processes) to tPC.

Our work thus raises an interesting empirical question as to what kinds of uncertainties various agents -- such as humans or animals -- actually appear to represent in dynamical inference tasks. For instance, it is not clear in the literature that subjective confidence ratings are highly correlated with the true dynamical uncertainty of the decisions in a task \citep{navajas2017idiosyncratic}.
%(-/?) and appear to display many unusual properties unexpected from the normative Bayesian account (-/citations Rafal).  
Although the Kalman filter has been used to describe reinforcement learning \citep{daw2006cortical}, direct comparison with simpler reinforcement learning models did not favour the Kalman filter \citep{howard2010neural}.
It may be interesting, therefore, to compare predictions of the Kalman filter (and the extended model representing posterior uncertainty \citep{kutschireiter2017nonlinear}) and the original tPC models to experimental data directly.
For example, one could compare if the learning rate in reinforcement learning tasks is better described by the Kalman Gain or the value from tPC. Another interesting line related to Kalman filtering is to compare nonlinear tPC with extended \citep{ribeiro2004kalman} and unscented \citep{julier2004unscented} Kalman filters in nonlinear tasks. Since this work focuses on linear tPC and its properties, we leave the exploration of these directions for future investigations.


\subsection{Generalized coordinates}
% A further interesting theoretical property of the model is its potential to autonomously learn to represent the dynamics of systems in generalized coordinates of motion \citep{friston2008DEM, friston2008hierarchical} if provided with generalized coordinates as inputs.
Generalized coordinates, introduced into the predictive coding literature by \citet{friston2008DEM}, and well known in engineering practice, involve representing the n'th order derivatives of a state as additional coordinate dimensions. In effect, a point in generalized coordinates of motion to n'th order reflects the approximate trajectory of the state given by the n'th order Taylor expansion around that point. 
Original predictive coding models involving generalized coordinates typically required hard-coding the relationships between the coordinates, and the relative precisions between different dynamical orders \citep{friston2008hierarchical,friston2008variational,parr2018discrete} which results in intricate and complex hard-coded connectivity, reducing the ultimate biological plausibility of such models. On the other hand, our tPC model can learn the $\bfW$ and $\bfF$ matrices directly from data, following Hebbian rules. It is therefore an intriguing future direction to explore whether these matrices encode the derivatives through learning.


\subsection{Machine learning models}

A further avenue for future work relates to the challenge of learning long-term dependencies which span over many time steps. This has long been a central challenge with these recurrent models, and emerges essentially due to the fact that information is permuted or lost at every step of the recurrent pass, and thus tracking dependencies across many recurrent loops becomes increasingly difficult \citep{hochreiter1997long,hochreiter1998vanishing,tallec2018can}. Numerous solutions to this have been suggested in the literature, ranging from specially designed cell architectures that can explicitly store or pass along information unaltered \citep{hochreiter1997long,tallec2018can} 
to having a nested hierarchy of recurrent models which allow for the propagation of information over longer and longer timescales \citep{koutnik2014clockwork}. Recent work has also addressed the question of how biological recurrent neural networks can be trained for temporal prediction \citep{bellec2020solution}.
They rearranged learning dynamics of BPTT and proposed that synapses maintain eligibility traces encoding to what extent they contributed to neural activity over time, and when combined with error signals, such traces enable effective credit assignment. Work by \citet{lotter2016deep} adapted deep recurrent neural networks to perform a kind of predictive coding whereby the network was trained to predict future  \emph{prediction errors} of each layer. They demonstrated that the resulting network was capable of correctly predicting sequences of video frames. While substantially scaling up predictive coding architectures to challenging machine learning tasks, the networks of \citet{lotter2016deep} diverged in many ways from classical predictive coding architectures, and also utilized many non-biologically plausible components from machine learning such as convolutional and LSTM layers as well as training their network with BPTT.



% \subsection{Conclusion}

% % In this work, we have proposed neural implementations of temporal predictive coding model in recurrent networks which utilize only biologically plausible Hebbian learning rules. We demonstrate that the model can successfully perform both linear and nonlinear tracking tasks successfully, along with online system identification of both transition and observation dynamics, and can come close to the optimal tracking performance of the Kalman filter. Moreover, by analyzing closely the mathematical relationship of the model to the Kalman filter, we have identified the representation of subjective uncertainty as the key difference between the Kalman filter and the temporal predictive coding networks.

% % \begin{revised}
% In this work, we have proposed temporal predictive coding, a model of how the processing of dynamic stimuli may take place in the brain. We have demonstrated that the tPC model performs simple and localised computations that afford a biologically plausible neural network implementation while approximating the Kalman filtering model in theory. We have empirically verified this approximation, and further demonstrated that tPC can reproduce motion-sensitive receptive fields in the visual areas of the brain. Moreover, the model can be straightforwardly extended to account for nonlinear dynamics, bringing an insight on how the brain may perform dynamic sensory processing that is often nonlinear in nature.
% % \end{revised}