\chapter{Supplementary: Temporal Predictive Coding}\label{app:chap4}

\section{Derivation of recursive Bayesian estimation}
% \begin{comment}
Here we provide detailed derivations of Equation \ref{eq:bayesian-filtering}, showing that the original Bayesian filtering problem can be simplified as a recursive estimation process that finally leads to Kalman filtering and our own temporal predictive coding:
\begin{eqnarray}\label{eq:bayesian-filtering-derive}
\begin{aligned}
    p (\mathbf{x}_k | \mathbf{y}_{1:k}) &= \int p(\mathbf{x}_{1:k} | \mathbf{y}_{1:k}) d\mathbf{x}_1,\dots,d\mathbf{x}_{k-1} \\
    &\stackrel{(a)}{\propto} 
    %\int p(\mathbf{x}_{1:k}, \mathbf{y}_{1:k}) d\mathbf{x}_1,\dots,d\mathbf{x}_{k-1} \\
    %&= 
    \int p(\mathbf{y}_{1:k} | \mathbf{x}_{1:k})p(\mathbf{x}_{1:k}) d\mathbf{x}_1,\dots,d\mathbf{x}_{k-1} \\
    &\stackrel{(b)}{=} \int \prod_{i=1}^k p(\mathbf{y}_i | \mathbf{x}_i)p(\mathbf{x}_i|\mathbf{x}_{i-1}) d\mathbf{x}_1,\dots,d\mathbf{x}_{k-1} \\
    &= p(\mathbf{y}_k |\mathbf{x}_k)\int p(\mathbf{x}_k|\mathbf{x}_{k-1}) \prod_{i=1}^{k-1} p(\mathbf{y}_i | \mathbf{x}_i)p(\mathbf{x}_i|\mathbf{x}_{i-1})d\mathbf{x}_1,\dots,d\mathbf{x}_{k-1} \\
    &= p(\mathbf{y}_k |\mathbf{x}_k)\int p(\mathbf{x}_k|\mathbf{x}_{k-1}) \left [ \int \prod_{i=1}^{k-1} p(\mathbf{y}_i | \mathbf{x}_i)p(\mathbf{x}_i|\mathbf{x}_{i-1})d\mathbf{x}_1,\dots,d\mathbf{x}_{k-2}\right ] d\mathbf{x}_{k-1} \\
    &\stackrel{(c)}{\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{aligned}
\end{eqnarray}
where step (a) is achieved by Bayes' rule, step (b) is achieved by the Markov assumption we made on the hidden states, and step (c) is achieved by observing that the integral in the square bracket in the penultimate line is proportional to the Bayesian filtering posterior at the previous step $k-1$. This can be observed by comparing it with the expression on the third line of the above equations. Importantly, this expression provides a recursive solution to the problem i.e., knowing the posterior at step $k-1$, we could estimate the posterior at step $k$. 

\section{Derivation of the update rules for the models}\label{appendix:update}
% \begin{comment}
Here we derive the update rule for $\hat{\mathbf{x}}_k$ that underlies both predictive coding and Kalman filtering (Equations \ref{eq:KF-correction} and \ref{eq:PC-update}). Notice that the objective functions for these two models follow a unified form:
\begin{eqnarray}\label{eq:F_general}
    \mathcal{F}_k = \frac{1}{2}(\mathbf{y}_k - \mathbf{F}\mathbf{x}_k)^{\top}\mathbf{\Sigma}_y^{-1}(\mathbf{y}_k - \mathbf{F}\mathbf{x}_k) + \frac{1}{2}(\mathbf{x}_k - \hat{\mathbf{x}}_k^-)^{\top} \mathbf{M}^{-1}(\mathbf{x}_k - \hat{\mathbf{x}}_k^-)
\end{eqnarray}
where $\mathbf{M} = \mathbf{\Sigma}_x$ for predictive coding and $\mathbf{M} = \mathbf{W}\mathbf{\Sigma}_{k-1}\mathbf{W}^{\top} + \mathbf{\Sigma}_x$ for Kalman filtering. To obtain $\hat{\mathbf{x}}_k$ that minimizes $\mathcal{F}_k$, we first take the derivative of $\mathcal{F}_k$ with respect to $\mathbf{x}_k$:
\begin{eqnarray}
    \begin{aligned}
        \frac{\partial \mathcal{F}_k}{\partial \mathbf{x}_k} = (\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1}\mathbf{F} + \mathbf{M}^{-1})\mathbf{x}_k - (\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1}\mathbf{y}_k + \mathbf{M}^{-1}\hat{\mathbf{x}}_k^-).
    \end{aligned}
\end{eqnarray}

Then, by setting the derivative to $0$ we have the optimal $\hat{\mathbf{x}}_k$:
\begin{eqnarray}
    \begin{aligned}
        \hat{\mathbf{x}}_k &= (\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1}\mathbf{F} + \mathbf{M}^{-1})^{-1}(\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1}\mathbf{y}_k + \mathbf{M}^{-1}\hat{\mathbf{x}}_k^-) \\ 
        &\stackrel{(a)}{=} [\mathbf{M} - \mathbf{M}\mathbf{F}^{\top}(\mathbf{\Sigma}_y + \mathbf{F}\mathbf{M}\mathbf{F}^{\top})^{-1}\mathbf{F}\mathbf{M}](\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1}\mathbf{y}_k + \mathbf{M}^{-1}\hat{\mathbf{x}}_k^-) \\
        &\stackrel{(b)}{=} [\mathbf{M} - \mathbf{K}\mathbf{F}\mathbf{M}](\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1}\mathbf{y}_k + \mathbf{M}^{-1}\hat{\mathbf{x}}_k^-) \\
        &= \hat{\mathbf{x}}_k^- - \mathbf{K}\mathbf{F}\hat{\mathbf{x}}_k^- + [\mathbf{M}\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1} - \mathbf{K}\mathbf{F}\mathbf{M}\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1}]\mathbf{y}_k \\
        &= \hat{\mathbf{x}}_k^- - \mathbf{K}\mathbf{F}\hat{\mathbf{x}}_k^- + [\mathbf{K}\mathbf{K}^{-1}\mathbf{M}\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1} - \mathbf{K}\mathbf{F}\mathbf{M}\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1}]\mathbf{y}_k \\
        &= \hat{\mathbf{x}}_k^- - \mathbf{K}\mathbf{F}\hat{\mathbf{x}}_k^- + \mathbf{K}[(\mathbf{\Sigma}_y + \mathbf{F}\mathbf{M}\mathbf{F}^{\top})\mathbf{F}^{-T}\mathbf{M}^{-1}\mathbf{M}\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1} - \mathbf{F}\mathbf{M}\mathbf{F}^{\top}\mathbf{\Sigma}_y^{-1}]\mathbf{y}_k \\
        &= \hat{\mathbf{x}}_k^- - \mathbf{K}\mathbf{F}\hat{\mathbf{x}}_k^- + \mathbf{K}\mathbf{y}_k \\
        &= \hat{\mathbf{x}}_k^- + \mathbf{M}\mathbf{F}^{\top}(\mathbf{\Sigma}_y + \mathbf{F}\mathbf{M}\mathbf{F}^{\top})^{-1}(\mathbf{y}_k - \mathbf{F} \hat{\mathbf{x}}_k^-).
    \end{aligned}
\end{eqnarray}

In step (a) we use the Woodbury matrix inversion identity and in step (b) we replace $\mathbf{M}\mathbf{F}^{\top}(\mathbf{\Sigma}_y + \mathbf{F}\mathbf{M}\mathbf{F}^{\top})^{-1}$ with $\mathbf{K}$. Substituting $\mathbf{M} = \mathbf{\Sigma}_x$ for predictive coding and $\mathbf{M} = \mathbf{W}\mathbf{\Sigma}_{k-1}\mathbf{W}^{\top} + \mathbf{\Sigma}_x$ for Kalman filtering we get Equations \ref{eq:KF-correction} and \ref{eq:PC-update} respectively.

\newpage
\section{Estimated position and velocity in the tracking task}
% \begin{comment}
Fig.~\ref{fig:vel_pos} shows the impact of the number of inference steps on the velocity and position dimensions in the tracking experiment. All values are with arbitrary units (a.u.).
\begin{figure}[!h]
    \centering
    \includegraphics[width=\linewidth]{figs/chap4/other2dims.pdf}
    \caption{{\bf Estimated position and velocity (a.u.) in the tracking experiments with tPC and Kalman filtering.} } 
    \label{fig:vel_pos}
\end{figure}

\newpage
\section{Tracking performance on the other two dimensions}\label{append:twodim-learning}
The following figure shows the performance on the other two dimensions in both hidden and observation levels in the tracking experiment, where $\bfW$ and $\bfF$ are learned. All values are with arbitrary units (a.u.).
% \begin{comment}
\begin{figure}[!h]
    \centering
    \includegraphics[width=\linewidth]{figs/chap4/twodims_learning.pdf}
    \caption{{\bf Estimated position and velocity (a.u.) in the tracking experiments with tPC and Kalman filtering, where $\bfW$ and $\bfF$ are learned.} } 
    \label{fig:vel_pos_learned}
\end{figure}

% \section{Experiments on a simple nonlinear model}
% % \begin{comment}
% In this Appendix, we focus on a simple simulation in which observations are generated from a probabilistic model of the same form as that of a nonlinear predictive coding network.
% In such a setting we expected the nonlinear temporal predictive coding network to learn to predict these data well and we examined whether it is able to do it better than a linear model.
% %: 1) by simulating a series of randomly generated synthetic data and then inferring the immediate future observation; 2) by simulating the nonlinear dynamics of a pendulum swinging and testing the model's learning capability of the generalized coordinates.

% %\subsubsection{Synthetic Stimuli}
% We generated a time series of 100, two-dimensional stimuli according to Eqs. $1$ and $2$ in the main text
% % by defining the parameters of $C$ and $A$ matrices to test whether the model can learn and update those parameters. The synthetic stimuli were simulated according to the nonlinear equations as,
% % \begin{align}
% %       \label{nonlinear_data_2}
% %        x_k &= A g(x_{k-1}) + \omega_x \\
% % %       \numberthis
% %        \label{nonlinear_data_1}
% %        y_k &= C f(x_k) + \omega_y
% % \end{align}
% where the nonlinear function $f$ is chosen to be the hyperbolic tangent function, and $\bm{\omega}_y$ and $\bm{\omega}_x$ correspond to standard Gaussian noise 
% %A similar neural architecture as in Figure \ref{fig:recurrent_multistep} was used to implement the recurrent predictive coding algorithm. For each simulation noise values were 
% sampled from an i.i.d. Gaussian distribution with a mean of 0 and a variance of 0.01. In this 2D example, parameters of $\bfF$ and $\bfW$ were set to a rotation matrix and the identity matrix multiplied by 3, respectively, which are as follows,
% \begin{eqnarray}
%     \begin{aligned}
%     \bfW &= \begin{bmatrix}
%     \frac{-1}{2}\Delta k & 1 \\
%     -1 & \frac{-1}{2}\Delta k
%     \end{bmatrix} \times 3 & \\
%     \bfF &= \begin{bmatrix}
%     1 & 0 \\
%     0 & 1
%     \end{bmatrix} \times 3.
% \end{aligned}
% \end{eqnarray}



% The multiplication by 3 in both matrices above was performed to produce a system that operates beyond the linear range of the hyperbolic tangent function (we used $\Delta k = 0.5$). 
% %As a result, 
% The observed stimuli are periodic signals which are simply noisy versions of the hidden variables. 
% %The variables and stimuli were generated with a sampling frequency of 2Hz to result in a relatively large $\Delta t$ (\textit{$\Delta t$ = 1 / sampling frequency}), by solving Equations \ref{nonlinear_data_1} and \ref{nonlinear_data_2} using Euler method with integration step of 0.5. Using a large $\Delta t$ makes the learning problem substantially more challenging. 

% Both linear and nonlinear models were trained to predict these stimuli. The models were simulated with a single state of inference between samples ($\Delta k= \Delta t$). For the learning process, the parameters of the $\bfF$ matrix were initialized as an identity matrix while the weights between hidden units were initialized to zeros. Additionally, each neuron $\hat{x}$ was initialized to 0, and the model was trained according to Equations 11 in the main text.


% Fig~S3A and B show zoomed-in time series with predictions of nonlinear and linear models (for the sake of visualization, only one dimension is shown).
% Fig~S3A demonstrates that in the first 50 time steps of the simulation, both models significantly move their predictions towards true time-series. Fig~S3B shows the final 50 time steps, where 
% %The training data is a noisy time-series and the neurons are trained to represent the current time input stimulus internally by continuously learning the hidden variables i.e., the $C$ and $A$ matrices. 
% the nonlinear model performed relatively better than the linear model. 
% %Looking at the left window, the nonlinear model keeps under-predicting the true signal and also straggles slightly from the ground-truth simulation but, as the learning proceeds, the nonlinear model begins to temporally adjust itself with the the stimulus and move its prediction to the current time step. Additionally, the nonlinear model no longer under-predicts the ground-truth signal and it also predicts the shape and the pattern of the signal more accurately. 

% \begin{figure}[!t]
%     \centering
%     \includegraphics[width=\linewidth]{figs/chap4/Fig-nonlin.pdf}
%     \caption{{\bf S3 Figure. Results with the simple nonlinear model.} A, B: Predictions in the first and final $50$ steps of the models. C: MSE of both models as a function of time. D: Fitted values with the $\sin$ function. E, F: MSE between the learnt $\bfW$s and the true $\bfW$ and visualizations of the learnt and true $\bfW$s.} 
%     \label{fig:nonlin}
% \end{figure}

% The observation prediction error of each model is depicted in Fig~S3C. 
% %The plots show the average errors from the 100 simulations, with solid lines representing the mean values and shaded areas representing the standard errors. 
% Within the first 50 to 100 time steps, the linear model has a much higher rate of reduction in its prediction error compared with the nonlinear model; however, this rapid drop in the error is followed by a stationary line and even a slow divergence from the optimum toward the end of the simulation. On the other hand, the nonlinear model shows a slower, but relatively steady, rate of error reduction which ultimately leads to convergence. 
% %The mean of the average errors are 0.04 and 0.2 for the nonlinear and linear models, respectively. The standard errors for the nonlinear model is 0.003 and for the linear model is 0.01. 

% The lower loss achieved by the nonlinear model could be explained by the fact that the ground-truth simulation data has a different oscillatory pattern as a result of the nonlinearity in the Eqs \ref{eq:general_transition} and \ref{eq:general_emission} that were used to generate such data. Therefore, the shape of the signal is different than one would expect to observe from a simple linear system. To show this, we simulated a data generation process based on the following equation:
% \begin{eqnarray}
%        \label{sin_data}
%        \hat{y}_k = a \sin(b k + c),
% \end{eqnarray}
% where $a$, $b$, and $c$ were approximated by a least squares optimisation so they match the ground-truth data. The comparison is depicted in Figure S3D, illustrating the difference in the shape of the time-series.

% To better understand the behaviour of each model, we investigated how well the transition matrix $\bfW$ is learnt in each case. In Fig~S3F, we plotted the mean of the learnt matrix by colour-coding the elements and normalising it to the minimum and maximum values of the true matrix. As expected, the parameters of the $\bfW$ matrix are closer to the true values when compared against their corresponding values in the linear model. 
% To quantify it, we also plotted the average errors of the learnt parameters for both models, showing that the linear model has a much higher error (almost by three folds) versus the nonlinear model - see Fig~S3E.