\section*{Appendix}

\subsection{Model Description}
\label{sec:model_details}
The GNCN model has $L+1$ layers of neurons $\mathfrak{N}^0,\mathfrak{N}^1,\dots, \mathfrak{N}^L$, where $\mathfrak{N}^0$ is the output layer. 
Each layer $\mathfrak{N}^\ell$ has $J^\ell$ neurons and each neuron has a latent state represented by a single number. The combined latent state of all neurons in  layer $\mathfrak{N}^\ell$ is represented by the  vector $\latent{z}^\ell$, with $\latent{z}^0$ being the same as the data $\mathbf{x}$. The network is interpreted as a specification of the probability $P(\latent{z}^0~|~\latent{z}^1)\cdots P(\latent{z}^{L-1}~|~\latent{z}^L)P(\latent{z}^L)$, similar to a Bayesian network. Thus layer $\mathfrak{N}^\ell$ represents the conditional probability $P(\latent{z}^{\ell}~|~\latent{z}^{\ell+1})$, with the last layer $\mathfrak{N}^L$ representing $P(\latent{z}^L)$.
For image data, the distribution $P(\latent{z}^0~|~\latent{z}^1)$ associated with the output layer is multivariate Bernoulli with mean vector $\mean{z}^0$ (which depends on $\latent{z}^1$). All of the other distributions $P(\latent{z}^{\ell}~|~\latent{z}^{\ell+1})$ are multivariate Gaussians with mean vector $\mean{z}^\ell$ (which depends on $\latent{z}^{\ell+1}$) and covariance matrix $\mathbf{\Sigma}^{\ell}$.
The mean vectors are obtained in a feed-forward manner from the latent state of the neighboring layer:
\begin{align}
    \mean{z}^\ell \gets g^\ell(\mathbf{W}^\ell \cdot  \phi^{\ell+1}(\latent{z}^{\ell+1})),\label{eqn:ff}
\end{align}
where $\mathbf{W}$ is a weight matrix, $g^\ell$ and $\phi^{\ell+1}$ are activation functions, and $\cdot$ indicates matrix multiplication. Thus each layer $\mathfrak{N}^\ell$ is specified by two functions $g^\ell$ and $\phi^\ell$, a trainable weight matrix $\mathbf{W}^\ell$ and a covariance matrix $\mathbf{\Sigma}^\ell$ (the last layer $\mathfrak{N}^L$ is specified just by the activation function $\phi^L$).
Note in Equation \ref{eqn:ff}, the mean vector $\mean{z}^\ell$ depends on the sampled realization  $\latent{z}^{\ell+1}$ from the previous layer, making this a hierarchical, continuous Gaussian mixture model.

\noindent
\textbf{Training the Model:} 
The complete data log-likelihood $\phi$ (also referred to as total discrepancy \cite{ororbia2020continual}) of the observed data $\mathbf{x}$ and the latent variables $\latent{z}^1, \dots, \latent{z}^L$ is defined formally as follows:
\begin{align}
\psi &=    \sum_j \Big(\mathbf{x}[j]\log \mean{z}^0[j] + (1-\mathbf{x}[j])\log (1-\mean{z}^0[j])\Big) + \sum_{\ell=1}^L \left(-\frac{1}{2}\log|(\Sigma^\ell)^{-1}| - \frac{1}{2}(\latent{z}^\ell - \mean{z}^\ell)^T \cdot (\Sigma^\ell)^{-1} \cdot (\latent{z}^\ell - \mean{z}^\ell)\right) 
\end{align}
Since all of the latent variables are continuous, the updates below follow the form of the exact gradient, i.e., backprop (allowing for gradient descent), or alternatives like feedback alignment \cite{lillicrap2016random} or LRA \cite{ororbia2018conducting}, to optimize the latent variables and the parameters.
The log likelihood has the following partial derivatives:
\begin{align}
    \frac{\partial \psi}{\partial (\Sigma^\ell)^{-1}} &= \frac{1}{2}\Sigma^\ell - \frac{1}{2}(\latent{z}^\ell - \mean{z}^\ell) \cdot (\latent{z}^\ell - \mean{z}^\ell)^T \label{eqn:precision_update} \\
    %
    \frac{\partial \psi}{\partial \Sigma^\ell} &= -\frac{1}{2}\Sigma^{-1} + \frac{1}{2}\Sigma^{-1} \cdot (\latent{z}^\ell - \mean{z}^\ell) \cdot (\latent{z}^\ell - \mean{z}^\ell)^T \cdot \Sigma^{-1} \label{covariance_update} \\
    %
    \frac{\partial \psi}{\partial \mathbf{W}^0} &= \left(\frac{\partial g^0}{\partial \mathbf{h}^0} \otimes(\mathbf{x} \oslash \mean{z}^0 - (\mathbf{1}-\mathbf{x})\oslash (\mathbf{1}-\mean{z}^0))\right) \cdot  \left(\phi^{1}(\latent{z}^{1})\right)^T  \label{eqn:W_0_update}\\
    &\phantom{=}\text{where }\oslash \text{ is element-wise division}, 
     \otimes \text{ is element-wise product, and } \mathbf{h}^\ell =\mathbf{W}^\ell \cdot \phi^{\ell+1}(\latent{z}^{\ell+1})\nonumber\\
    %
    \frac{\partial \psi}{\partial \mathbf{W}^\ell} &=\frac{\partial \psi}{\partial \mathbf{h}^\ell} \cdot  \left(\phi^{\ell+1}(\latent{z}^{\ell+1})\right)^T \label{eqn:W_l_update}\\
    %
    \frac{\partial \psi}{\partial \latent{z}^1} &= \frac{\partial}{\partial \latent{z}^1}\left( \sum_j \Big(\mathbf{x}[j]\log \mean{z}^0[j] + (1-\mathbf{x}[j])\log (1-\mean{z}^0[j])\Big)\right)
                    - (\Sigma^{\ell})^{-1} \cdot (\latent{z}^{1}-\mean{z}^{1}) \label{eqn:z_0_deriv} \\
   
    \frac{\partial \psi}{\partial \latent{z}^\ell} &=  
                    \left( \frac{\partial \mean{z}^{\ell-1}}
                    {\partial \latent{z}^\ell} \cdot \Big( 
                    (\Sigma^{\ell-1})^{-1} \cdot (\latent{z}^{\ell-1}-\mean{z}^{\ell-1}) \Big) \right)-
                    (\Sigma^{\ell})^{-1} \cdot (\latent{z}^{\ell}-\mean{z}^{\ell}) \label{eqn:z_l_deriv}
   
\end{align}
In this work, we incorporate two key concepts from LRA \cite{ororbia2018conducting,ororbia2019biologically}: 1) the use of error synapses to directly resolve the weight-transport problem, and 2) the omission of derivatives of activation functions which yield synapse rules that function like error-Hebbian updates. To incorporate these modifications, first note that the error neurons themselves are derived from the likelihood function:
\begin{align}
    \mathbf{e}^0 &= \frac{\partial \psi }{\partial \mean{z}^0 } = (\mathbf{x} \oslash \mean{z}^0 - (\mathbf{1}-\mathbf{x})\oslash (\mathbf{1}-\mean{z}^0))  \label{eqn:e_0} \\
    \mathbf{e}^\ell &= \frac{\partial \psi }{\partial \mean{z}^\ell } =  (\Sigma^{\ell-1})^{-1} \cdot (\latent{z}^{\ell-1}-\mean{z}^{\ell-1}) \label{eqn:e_l}
\end{align}
and are implemented as separate sets of activities (like the state neurons). To carry the activities/messages of these particular neurons in order to calculate the update to $\latent{z}^\ell$, we replace the term $\frac{\partial \latent{z}^{\ell-1}}{\partial \latent{z}^\ell}$ with a learnable error matrix $\mathbf{E}^\ell$. This substitution allows us to ultimately rewrite Equations \ref{eqn:z_0_deriv} and \ref{eqn:z_l_deriv} as follows:
\begin{align}
   \frac{\partial \psi}{\partial \latent{z}^1}  &\approx \Delta \latent{z}^1 = \mathbf{E}^1 \cdot \Big( \frac{g^0}{\mathbf{h}^0} \otimes e^0 \Big) - \mathbf{e}^1 = \mathbf{E}^1 \cdot \Big( \latent{z}^\ell - \mean{z}^\ell  \Big) - \mathbf{e}^1 \label{eqn:z_1_delta} \\
    \frac{\partial \psi}{\partial \latent{z}^\ell} &\approx \Delta \latent{z}^\ell =  \mathbf{E}^\ell \cdot \mathbf{e}^{\ell-1} - \mathbf{e}^\ell \label{eqn:z_l_delta}
\end{align}
noting that no partial derivatives of the activation function $\phi^\ell$ is required -- we are using an error-Hebbian style update rule, which has been shown to been effective empirically and crucially removes the need for state neurons to be aware of their own point-wise activation derivative. Note that this means $\phi^\ell$ could easily be a discrete function, such as the signum or Heaviside.
The synaptic matrix updates in Equations \ref{eqn:W_0_update} and \ref{eqn:W_l_update} can be written in terms of the error neurons:
\begin{align}
    \frac{\partial \psi }{\partial \mathbf{W}^0} &\propto \Delta \mathbf{W}^0 = \mathbf{e}^0 \cdot (\phi^1( \mathbf{z}^1) )^T \label{eqn:W_0_delta} \\
     \frac{\partial \psi }{\partial \mathbf{W}^\ell} &\propto \Delta \mathbf{W}^\ell = \mathbf{e}^\ell \cdot (\phi^{\ell + 1}( \mathbf{z}^{\ell+1}) )^T \label{eqn:W_l_delta}
\end{align}
and, following in line with LRA, the error synapses can be updated as:
\begin{align}
    \Delta \mathbf{E}^\ell = \lambda \Big( \phi^{\ell+1}(\latent{z}^{\ell+1}) \cdot  (\mathbf{e}^\ell)^T \Big) \label{eqn:E_l_update}
\end{align}
where $\lambda \in [0,1]$ controls the strength of the error synaptic adjustment (usually set to values $< 1$).

Putting all of the components above together, given a data point(s) $\mathbf{x}$, we train the GNCN following an online alternating maximization approach as follows:
\begin{enumerate}
    \item[] // \textbf{Initialization:}
    \item Set $\latent{z}^0 = \mathbf{x}$, and $\latent{z}^\ell = \mathbf{0}, \; \forall \ell \ge 1$
    \item Compute mean estimates $\mean{z}^\ell, \; \forall \ell \ge 0$ (Equation \ref{eqn:ff}) \& use Equation \ref{eqn:e_0} and Equation \ref{eqn:e_l} to initialize $\mathbf{e}^\ell, \; \forall \ell \ge 0$
    \item[] // \textbf{Latent Update Step:} Search for most probable value of $\latent{z}^\ell, \forall \ell$
    \item Use Equation \ref{eqn:z_1_delta} to get $\Delta \latent{z}^1$ \& Equation \ref{eqn:z_l_delta} to get $\Delta \latent{z}^\ell, \; \forall \ell > 1$. The latent states are modified as follows:
    \begin{align*}
        \latent{z}^1 &= \latent{z}^1 + \beta \Big( \Delta \latent{z}^1 - \mathbf{V}^1 \cdot \latent{z}^1 - \gamma \latent{z}^1 \Big) \\
        \latent{z}^\ell &= \latent{z}^\ell + \beta \Big( \Delta \latent{z}^\ell - \mathbf{V}^\ell \cdot \latent{z}^\ell - \gamma \latent{z}^\ell \Big), \; \forall \ell > 1
    \end{align*}
    \item Run Equation  \ref{eqn:e_0} and  \ref{eqn:e_l} to obtain updated error neurons, i.e., $\mathbf{e}^\ell, \; \forall \ell \ge 0$
    \item Repeat Steps 3 and 4 for $T$ iterations
    \item[] // \textbf{Parameter Update Step:} Update parameters given estimated latent states
    \item Update synaptic matrices using Equations \ref{eqn:W_0_delta}, \ref{eqn:W_l_delta}, \ref{eqn:precision_update}, and \ref{eqn:E_l_update} \mbox{.}
\end{enumerate}
Note that at test time, to reconstruct $\mathbf{x}$, one should only use Steps 1-5 of the recipe above. $\gamma$ controls the strength of the decay/leak applied to $\mathbf{z}^\ell$ and corresponds to placing an additional $\mathcal{N}(\mu=0,\Sigma=\lambda \mathbf{I})$ prior over the latent states.

Above, note that we also introduce a (fixed) lateral competition matrix $V^\ell$ that affects how a latent state $\latent{z}^\ell$. It is created to contain the self-excitation weights and lateral inhibition weights by using the following matrix equation: $V^\ell = \alpha_h (M^\ell) \otimes (1 - I) - \alpha_e (I) $, where $I$ is the identity matrix and the masking matrix $M^\ell \in \{0,1\}^{J_\ell \times J_\ell}$ is set by the experimenter (placing ones in the slots where it is desired for neuron pairs to laterally inhibit one another). In this study, we set $\alpha_e = 0.13 $ (the self-excitation strength) and $\alpha_h = 0.125$ (the lateral inhibition strength). Our mask matrix $M^\ell$, which emphasized a type of group or neural-column form of competition, was generated by the following process:
\begin{enumerate}
\item create $J_\ell / K$ matrices of shape $J_\ell \times K$ of zeros, i.e., $\{S_1, S_2, \cdots, S_k, \cdots, S_{J_\ell / K}\}$
\item in each matrix $S_k$  insert ones at all combinations of coordinates $c = \{ 1, \cdots, k, \cdots, K \}$ (column index) and $r = \{ 1 + K * (k-1), \cdots, k + K * (k-1), \cdots, K + K * (k-1) \}$ (row index)
\item concatenate the $J_\ell / K$ matrices along the horizontal axis, i.e., $M^\ell = <S_1, S_2, \cdots, S_C>$.
\end{enumerate}

\subsection{Baseline Model Descriptions}
\label{sec:baseline_details}

To make learning the decoder function ($\mbox{NN}$) described in the main paper tractable, it is common practice in the deep learning literature to introduce a supporting function known as the \emph{encoder} \cite{kingma2013auto}. The encoder ($\mbox{NN}_e$), parameterized by a feedforward network, takes in the input stimulus $\mathbf{x}$ and maps it to $\mathbf{z}$ or to a distribution over $\mathbf{z}$. Depending on the choice of encoder, one can recover one of the four main baselines we experimented with in this paper.

For all backprop-based baseline models in this paper, the decoder of each was regularized with an additional L2 penalty. Specifically, this meant that their data log likelihood objectives always took the form: $\psi_{reg} = \psi + \Omega(\Theta_{\mbox{NN}})$, where $\Theta_{NN} = \{\mathbf{W}^L, \cdots, \mathbf{W}^\ell, \cdots, \mathbf{W}^1 \}$ contains all of the weight matrix parameters of the decoder $\mbox{NN}$. $\Omega_{\mbox{NN}}$ is the regularization function applied to the decoder, i.e., $\Omega_{\mbox{NN}} = -\lambda \sum_{\mathbf{W}^\ell \in \Theta} ||\mathbf{W}^\ell||^2_2 $ where $||\mathbf{W}^\ell||_2$ denotes computing the Frobenius norm of $\mathbf{W}^\ell$. During training/optimization with gradient ascent, we do not constrain the column norms of any of the weight matrices for any of the baseline models (as we do for the GNCN) as we found that doing so worsened their generalization ability. 

Furthermore, the number of total layers in the decoder for any model was set to be four -- one output and one input layer with two hidden layers in between. The encoder was constrained to be the same -- one input and one output layer with layers in between (in the case of the GVAE, CV-GVAE, and GAN-AE, the encoder's output is technically split into two blocks, as described later). The sizes of the hidden layers were set such that the total number of learnable model weights were approximately equal across all baselines and GNCNs (maximum was $1,400,000$ synapses), which means that all models were forced to have the same parameter complexity to avoid any unfair advantages that might come from over-parameterization.

\noindent 
\textbf{Regularized Auto-encoder (RAE):} The encoder $\mbox{NN}_e$ is designed to be a feedforward network of $L$ layers of neurons. Each layer is a nonlinear transformation of the one before it, where $\mathbf{\hat{z}}^\ell = \phi^\ell(\mathbf{E}^\ell \cdot \mathbf{\hat{z}}^{\ell-1})$. Like in the decoder, $\phi^\ell$ is an activation function and $\mathbf{E}^\ell$ is a set of tunable weights. In this paper, we chose $\phi^\ell$ to be the linear rectifier, i.e., $\phi^\ell(v) = \mbox{max}(0, v)$. The bottom layer activation was chosen to be the logistic link, i.e., $\phi^0(\mathbf{z}) = 1/(1 + exp(-\mathbf{z})$.

Note that in the RAE, the input to the decoder is now $\mathbf{z} = \mathbf{\hat{z}}^L$, i.e., the noise sample vector is set equal to top-most layer of neural activities of the encoder.
The data log likelihood that the RAE optimizes is:
\begin{align}
    \psi = \sum_j \Big(\mathbf{x}[j]\log \mathbf{z}^0[j] + (1-\mathbf{x}[j])\log (1-\mathbf{z}^0[j])\Big) 
\end{align}
where updates to each weight matrix $\mathbf{E}^\ell$ (of the encoder) and $\mathbf{W}^\ell$ (of the decoder) are updated by computing the relevant gradients $\frac{\partial \psi}{\partial \mathbf{E}}$ and $\frac{\partial \psi}{\partial \mathbf{W}}$, respectively. The weight gradients are then used to update model parameters via gradient ascent.

\noindent 
\textbf{Gaussian Variational Auto-encoder (GVAE):} Instead of using an encoder to only produce a single value for $\mathbf{z}$, we could instead modify this network to produce the parameters of a distribution over $\mathbf{z}$ instead. If we assume that this distribution is a multivariate Gaussian with a mean $\mu_z$ and a diagonal covariance $\sigma^2_z = \mathbf{\Sigma}_z \otimes \mathbf{I}$, we can then modify the RAE's encoder function to instead be: $(\mathbf{\mu}_z, \sigma^2_z) = \mbox{NN}_e(\mathbf{x})$. Specifically, the top-most layer of $\mbox{NN}_e(\mathbf{x})$ is actually split into two separate output layers as follows: $\mu_z = \mathbf{E}^L_\mu \cdot \mathbf{z}^{L-1}$ and $\sigma^2_z = \exp( \mathbf{E}^L_\sigma \cdot \mathbf{z}^{L-1} )$ (this is also known as the variational autoencoder, or VAE \cite{kingma2013auto}).
$\mathbf{E}^L_\mu$ is the tunable weight matrix for the mean and $\mathbf{E}^L_\sigma$ is the tunable weight matrix for the variance.
The data log likelihood for the GVAE is as follows:
\begin{align}
    \psi = \sum_j \Big(\mathbf{x}[j]\log \mathbf{z}^0[j] + (1-\mathbf{x}[j])\log (1-\mathbf{z}^0[j])\Big) - \mbox{D}_{KL}\Big( q(\mathbf{z}|\mathbf{x}) || p(\mathbf{z}) \Big) \label{eqn:vae_obj}
\end{align}
where $q(\mathbf{z}|\mathbf{x}) = \mathcal{N}(\mu_z, \sigma^2_z)$ (the Gaussian parameters produced by the encoder $\mbox{NN}_e(\mathbf{x})$) and $p(\mathbf{z}) = \mathcal{N}(\mu_p,\sigma^2_p)$ where $\mu_p = 0$ and $\sigma^2_p = 1$ (an assumed unit Gaussian prior over $\mathbf{z}$). The second term in the above objective is the Kullback-Leibler divergence $D_{KL}$ between the distribution defined by the encoder and the assumed prior distribution. This term is meant to encourage the output distribution of the encoder to match a chosen prior distribution, acting as a powerful probabilistic regularizer over the model's latent space. Note that this divergence term serves as a top-down pressure on the top-most layer of the encoder while the gradients that flow from the encoder (via the chain rule of calculus) act as a sort of bottom-up pressure. Note that since $NN_e(\mathbf{x})$ is a distribution, the input to the decoder is, unlike the RAE, a sample of the encoder-controlled Gaussian, i.e., $\mathbf{z} = \mu_z + \sqrt{\sigma^2_z} \otimes \epsilon$ where $\epsilon \sim \mathcal{N}(0,1)$.

Gradients of the likelihood in Equation \ref{eqn:vae_obj} are then taken with respect to all of the encoder and decoder parameters, including the new mean and variance encoder weights $\mathbf{E}^L_\mu$ and $\mathbf{E}^L_\sigma$, which are subsequently updated using gradient ascent. All the other activation functions of the GVAE are set to be the linear rectifier, except for the output function $\phi^0$ of the decoder, which, like the RAE, is set to be the logistic sigmoid.

\noindent 
\textbf{Constant-Variance Gaussian Variational Auto-encoder (CV-GVAE):} This model \cite{ghosh2020variational} is identical to the GVAE except that the variance parameters $\sigma^2_z$ of the encoder are omitted and a fixed (non-learnable) value is chosen instead for the variance (meaning that the diagonal covariance is collapsed further to a single scalar). The exact value for this variance meta-parameter, for each benchmark, was chosen from the range $[0.025, 1.0]$ by tuning performance to a held-out set of image samples. 

\noindent
\textbf{Generative Adversarial Network Autoencoder (GAN-AE):} This model, also referred to as an adversarial autoencoder \cite{makhzani2016adversarial}, largely adheres to the architecture of the GVAE except that the second term, i.e., the Kullback-Leibler divergence term, in the data log likelihood is replaced with the adversarial objective normally used to train implicit density estimators like the generative adversarial network (GAN) \cite{goodfellow2014generative}. As a result, we integrate a third feedforward network, i.e., $p_r = \mbox{NN}_d(\mathbf{z})$, into the generative model (this module is also referred to as the discriminator). The discriminator is tasked with distinguishing whether an input vector comes from the desired prior distribution $p(\mathbf{z})$ (set to be a unit Gaussian as in the GVAE) or comes from the encoder network $\mbox{NN}_e(\mathbf{z})$ distribution. This task is posed as a binary classification problem, where a sample from the encoder $\mathbf{z}_f \sim \mathcal{N}(\mu_z, \sigma^2_z)$ is assigned the label of $c = 0$ (fake sample) and a sample drawn from the prior $\mathbf{z}_r \sim \mathcal{N}(0, 1)$ is assigned a label of $c = 1$ (real sample). These fake and real samples are fed through the discriminator which returns a scalar value for each, representing the probability $p_r = p(c=1|\mathbf{z})$. 
This leads to the modified data log likelihood objective below:
\begin{align}
    \psi = \sum_j \Big(\mathbf{x}[j]\log \mathbf{z}^0[j] + (1-\mathbf{x}[j])\log (1-\mathbf{z}^0[j])\Big) + \Big( \log(\mbox{NN}_d(\mathbf{z}_r)) + ( 1 - \log(\mbox{NN}_d(\mathbf{z}_f)) ) \Big) \mbox{.} \label{eqn:adv_obj}
\end{align}
However, to update the weights of the GAN-AE, we do not compute partial derivatives of Equation \ref{eqn:adv_obj} directly. Instead, following in line with the typical multi-step optimization of \cite{goodfellow2014generative,makhzani2016adversarial}, upon presentation of a sample or mini-batch of samples, we compute gradients with respect to $\mbox{NN}_e(\mathbf{x})$, $\mbox{NN}(\mathbf{z})$, and $\mbox{NN}_d(\mathbf{z})$ separately. Specifically, if we group all of the encoder weights under $\Theta_{NN_e}$, all of the decoder weights under $\Theta_{NN}$, and all of the discriminator weights $\Theta_{NN_d}$, then the gradients of the objective are computed in three separate but successive steps shown below:
\begin{align}
    \Delta_{auto} &= \frac{\partial \sum_j \Big(\mathbf{x}[j]\log \mathbf{z}^0[j] + (1-\mathbf{x}[j])\log (1-\mathbf{z}^0[j])\Big)}{\partial (\Theta_{NN_e} \cup \Theta_{NN})} & \mbox{(Autoencoder gradients)} \label{eqn:enc_grad} \\
    \Delta_{gen} &= \frac{\partial \Big( 1 - \log(\mbox{NN}_d(\mathbf{z}_f)) \Big)}{\partial \Theta_{NN}}  & \mbox{(Generator gradients)} \label{eqn:gen_grad} \\
    \Delta_{disc} &= \frac{\partial \Big( \log(\mbox{NN}_d(\mathbf{z}_r)) + ( 1 - \log(\mbox{NN}_d(\mathbf{z}_f)) ) \Big)}{\partial \Theta_{NN_d}}  & \mbox{(Discriminator gradients)} \label{eqn:disc_grad}
\end{align}
where the above gradient calculations are each followed by a separate gradient ascent update to their relevant target parameters, i.e., $\Delta_{auto}$ is used to update $\Theta_{NN_e} \cup \Theta_{NN}$, $\Delta_{gen}$ is used to update $\Theta_{NN_e}$, and $\Delta_{disc}$ is used to update $\Theta_{NN_e} \cup \Theta_{NN_d}$.

The number of synaptic weights associated with the discriminator were included in the model's total parameter count and had two hidden layers of linear rectifier units. Again, like the GVAE, the hidden layer functions $\phi^\ell$ of the encoder and decoder were chosen to be the linear rectifier and the output $\phi^0$ of the decoder was set to be the logistic sigmoid.




\section{Introduction}
\label{sec:intro}

One way to understand how the brain adapts to its environment is to view it as a type of generative pattern-creation model \cite{friston2006free}, one that is engaged in a never-ending process of self-correction, often without external teaching signals (or labels) \cite{ororbia2020continual}. Under this perspective, the brain is continuously making predictions about elements of its environment, a process that allows it to infer useful representations of the sensory data it receives \cite{parr2018anatomy} as well as to synthesize novel patterns, which could serve as the potential basis for long-term planning and imagination itself \cite{clark2015surfing}.
From the theoretical viewpoint of predictive processing, the brain could be likened to a hierarchical model whose levels are implemented by neurons (or clusters of neurons). If levels are likened to regions of the brain, the neurons at one level (region) attempt to predict the state of neurons at another level (region) and adjust/correct their local model synaptic parameters based on how different their predictions were from the observed signal. Furthermore, these neurons utilize various mechanisms to laterally stimulate/suppress each other \cite{liang2017interactions} to facilitate contextual processing (such as grouping/segmenting visual components of objects in a scene).
As we will demonstrate in this article, this viewpoint can be turned into a powerful framework for learning generative models.

In machine learning, one central goal is to construct agents that learn distributed representations that extract the underlying structure of data, without the use of explicit supervisory signals such as human-crafted labels, i.e., unsupervised learning. 
Generative models, or models capable of synthesizing instances of data that resemble a database of collected patterns, that are based on deep artificial neural networks (ANNs), e.g., variational autoencoders \cite{kingma2013auto} or generative adversarial networks \cite{goodfellow2014generative}, have been shown to be one way of acquiring these representations.
Once trained, an ANN model is used to ``fantasize'' patterns by injecting it with noise and propagating this noise through the system until the output nodes are reached.

However, despite the success in deploying ANNs as generative models across a variety of applications, the way that ANNs operate and learn is a far-cry from the neuro-mechanistic story we described earlier \cite{crick1989recent,zador2019critique}. Specifically, ANN generative models are trained with the popular workhorse algorithm known as back-propagation of errors (backprop) \cite{rumelhart_learning_1986}, which is an elegant mathematical solution to the credit assignment problem in deep networks -- synaptic weights are adjusted through the use of teaching signals that are created by propagating an error, which exists exclusively at the output of the ANN, backwards along a feedback pathway \cite{ororbia2019biologically} (a path created by re-using the same weights that transmitted signals forward \cite{grossberg1987competitive}) .
By virtue of this formulation, backprop imposes the constraint that the ANN take the form of a directed feedforward structure (and does not permit the use non-differentiable activation functions and other mechanisms such as lateral connectivity). While the neurons in an ANN are usually arranged hierarchically, they do not make local predictions and they do not laterally affect each other's activity. Furthermore, synaptic adjustment in backprop-based models is done non-locally, while in neurobiological networks this adjustment is often argued to be done locally \cite{hebb1949organization,magee1997synaptically,bi2001synaptic,isomura2018error} (there are far more local connections than long-range connections \cite{zhang2000universal} with the neocortex adhering a local connectivity pattern \cite{zador2019critique}) -- that is, neurons make use of the information immediately available to them (in both time and space) and do not wait on distant regions in order to adjust their synapses.

While the question as to how credit assignment is exactly implemented in the brain is an open one, it would prove useful to machine learning, (computational) neuroscience, and cognitive science to have a framework that demonstrates how a neural system can learn something as complex as a generative model without backprop, using mechanisms and rules that are brain-inspired. 
In this work, inspired by the theory of predictive processing \cite{von1855ueber,bastos2012canonical,swanson2016predictive,clark2015surfing}, we propose the neural generative coding (NGC) computational framework as a promising way to learn generative ANNs, resolving several the key backprop-centric issues described above.
We demonstrate that NGC models not only remain competitive with backprop-based generative ANNs, such as the variational autoencoder \cite{kingma2013auto}, in terms of pattern creation but that they also significantly outperform these models on tasks that they are not directly trained for, such as classification and pattern completion.
As a result, our work presents promising evidence that brain-inspired alterations to traditional deep learning techniques can be a viable source of performance gains. 


\section{Results}
\label{sec:results}

\subsection{Generative Neural Coding Learns Viable Auto-Associative Generative Models}
\label{sec:results_gen_models}

\textbf{Problem Setting:} We start with a description of the problem setting -- an agent must learn to approximate a probability distribution  from a dataset $\mathbf{X}$ of samples. For notational reasons, this dataset is presented in column-major order, so that each column $\mathbf{x}$ represents a record (also known as an example or item).  $\mathbf{X}$ has $D$ rows and and $S$ total columns. The items are assumed independent, so that $p(\mathbf{X}) = \prod_{\mathbf{x} \in \mathbf{X}} p(\mathbf{x})$ and $\log p(\mathbf{X}) = \sum_{\mathbf{x} \in \mathbf{X}} \log p(\mathbf{x})$. We are interested in  directed generative models that are capable of producing \emph{explicit} density estimates of the data distribution and we will leave the examination of most \emph{implicit} density estimators based on generative adversarial networks \cite{goodfellow2014generative} for future work. 

\begin{figure}[t]
\begin{center}
\begin{subfigure}{0.5\textwidth}
  \centering
 
  \includegraphics[width=7.8cm]{figs/backprop_path.pdf}  
 
  \label{fig:backprop_path}
\end{subfigure}%
\begin{subfigure}{0.5\textwidth}
  \centering
 
  \includegraphics[width=7.1cm]{figs/ncn_path.pdf}  
 
  \label{fig:ncg_path}
\end{subfigure}
\end{center}
\vspace{-0.2cm}
\caption{
(Left) Credit assignment in backprop requires a strict, global feedback pathway, which requires the completion of the forward pass that carries information upstream (right to left), that carries the error message $\mathbf{e}^0$ at the output layer back along (left to right) the same synapses used in the forward pass to update downstream neurons  $\mathbf{\bar{z}}^1$ and $\mathbf{\bar{z}}^2$.
(Right) Our proposed neural generative coding (NGC) model sidesteps this neurobiologically implausible requirement by learning with short, local error transmission pathways made possible through recurrent error synapses and stateful neural activities.
Credit assignment under NGC operates with local mismatch signals, $\mathbf{e}^1$ and $\mathbf{e}^2$, that readily communicate this information to their respective layers, $\mathbf{\bar{z}}^1$ and $\mathbf{\bar{z}}^2$. Black arrows indicate forward propagation while red arrows indicate backwards transmission. Solid lines indicate that a signal is transformed along a synapse(s) while dashed arrows indicate direct copying of information. $f^\prime$ shows communication of the neuron's ``first derivative'', $\Delta$ represents the computed change to the synapse (of the forward pass) that will use the (nearby) error signal, and $\otimes$ indicates multiplication of the incoming signals.
}
\label{fig:bp_v_ngc}
\vspace{-0.5cm}
\end{figure}

\textbf{The Typical Deep Learning Approach:} In modern-day deep learning practice, a feedforward ANN, also called a \emph{decoder}, would be constructed to model the desired input distribution. The decoder ($\mbox{NN}$) takes as input a noise vector or a sampled latent variable  $\mathbf{z}$ and maps it into a data sample $\mathbf{x}=\mbox{NN}(\mathbf{z})$. This artificial neural network would typically be made up of $L$ layers of neurons, where the state in layer $\ell$ is represented by a vector $\mathbf{z}^\ell$. Each layer $\ell$ is interpreted as a transformation of the layer before it. That is, setting the state $\mathbf{z}^L$ of layer $L$ to be the same as the input noise vector $\mathbf{z}$, we get $\mathbf{z}^{L-1}=\phi^{L-1}(\mathbf{W}^L \cdot \mathbf{z}^L)$, where $\cdot$ indicates matrix multiplication, $\phi^{L-1}$ is an activation function and $\mathbf{W}^L$ are tunable weights. In general,  ${\mathbf{z}^\ell= \phi^\ell(\mathbf{W}^{\ell+1} \cdot \mathbf{z}^{\ell+1}) }$.
The output $\mathbf{z}^0$ of this decoder (Figure \ref{fig:bp_v_ngc} Left) would be the parameters of a probability distribution, such as the mean of a Bernoulli distribution, or mean and covariance of a Gaussian (see Appendix for descriptions of the backprop-based networks used in this study). One can sample from this distribution to get a sample point $\mathbf{x}$ (or use the mean of the distribution directly).

To fit this model to the data, one would choose the weight parameters $\mathbf{W}^\ell$ to minimize a loss function $\psi$ such as the negative log-likelihood, typically using some variant of stochastic gradient descent. Often the backprop algorithm is used to compute the partial derivatives of  $\frac{\partial \psi}{\partial \mathbf{W}^\ell}$ needed for this optimization. 
Computing the necessary derivatives according to backprop entails first computing an error signal $\mathbf{e}^0$ at the output (downstream) layer, or $\frac{\partial \psi}{\partial \mathbf{z}^0}$. This error signal is then transmitted  to internal (upstream) neurons by carrying this signal back along the forward synapses that were originally used to  transform $\mathbf{z}^L$ (in short, this is done by multiplying the signal with the transpose of the forward weight matrices). Furthermore, knowledge of the derivative of each activation function $\phi^\ell$ is required during these computations (as shown in Figure \ref{fig:bp_v_ngc}, Left).


\textbf{Backprop-Learning versus Brain-Like Learning:} While the backprop algorithm described above has proven to be popular and effective in training ANNs (including generative models \cite{lecun2015deep}), there are several issues with this way of conducting credit assignment when considering how the brain might learn. First, according to backprop: 1) synapses that make up the forward information pathway need to directly be used in reverse to communicate teaching signals (the weight transport problem \cite{grossberg1987competitive}), 2) neurons need to be able to know about and communicate their own activation function's first derivative, 3) neurons must wait ``patiently'' for the neurons ahead of them to percolate their error signals way back so they know when and how to adjust their own synapses (the update-locking problem \cite{jaderberg2017decoupled}), 4) there is a distinct form of information propagation for error feedback, one that starts from the system's output and works its way back to the input layer (see Figure \ref{fig:bp_v_ngc} Left), which does not influence neural activity (the global feedback pathway problem \cite{ororbia2019biologically}), i.e., backprop creates signals that only affect weights but do not (at least directly) affect or improve the network's representations of the environment. However, these properties/constraints inherent to backprop do not conform to known biological feedback mechanisms underlying neural communication in the brain \cite{crick1989recent}. 

To start, the brain does not reuse its forward synapses to carry backward error information like backprop does, i.e., backprop requires that feedback and forward synaptic pathways are exactly tied together \cite{grossberg1987competitive}, which practically has been shown to create memory access pattern issues in hardware implementations \cite{crafton2019local}.
Instead, the brain is heavily recurrently connected \cite{andersen1963hippocampus,douglas1995recurrent} 
which allows for complementary circuits/pathways to form that would allow for percolation of error/mismatch information \cite{long2016hippocampal,wacongne2012neuronal}. 
In further contrast to backprop, it is more commonly accepted that neurons in the brain learn ``locally'' \cite{hebb1949organization,el2018locally},
i.e., they operate with only immediately available information (such as their own activity and that of nearby neurons that they are connected to) and error signals that are related to their operation directly, and no evidence has been presented that they also sport specialized circuitry to communicate the derivative of a loss function with respect to their activities \cite{hinton2007recognize} (in fact, this restriction imposed by backprop prevents the use of discrete-valued activation function, whereas real neurons communicate with binary spike signals). The brain, unlike backprop, does not rely on an explicit, long ``global feedback pathway'' \cite{ororbia2019biologically} to adjust neural activities and synaptic strengths, which is why the updates to synapses can happen  quickly and even likely asynchronously \cite{zeki2015massively}. Furthermore, even the form of its forward propagation does not necessarily happen in a strict sequential manner\cite{zeki2015massively}, i.e., it does not occur layer-by-layer. These two properties allow the brain to side-step the update-locking problem facilitating a great deal of parallelism and significantly reduced memory overhead for inference and synaptic weight adjustment. In addition, since the brain does not require itself to be a differentiable, acyclic graph (which is a strict requirement for backprop), neurons can be laterally connected to each other, which is a wiring pattern that has been demonstrated to facilitate contextual processing of sensory input \cite{adesnik2010lateral,liang2017interactions}

Finally, the global feedback pathway itself is one key source behind the well-known exploding and vanishing gradient problems \cite{pascanu2013difficulty}
in very deep ANNs, yielding unstable or ineffective learning unless specific heuristics are used.

\textbf{The Neural Generative Coding Approach:} In contrast to the backprop-based way of designing and training ANN generative models, our proposed framework, neural generative coding (NGC), provides one way to emulate the several neuro-biological principles and properties described above by proposing a model and its corresponding training procedure.
In this framework, the model is referred to as the generative neural coding network (GNCN). It has  a set of $L+1$ layers of neurons (also called state variables) $\mathfrak{N}^0,\mathfrak{N}^1,\dots, \mathfrak{N}^L$, where $\mathfrak{N}^0$ is the output layer. The combined latent state of all neurons in layer $\mathfrak{N}^\ell$ is represented by the  vector $\latent{z}^\ell$. The state $\latent{z}^0$ of the output layer is set to equal the data $\mathbf{x}$. We use the notation  $\mathbf{Z} = \{ \mathbf{z}^1,\cdots,\mathbf{z}^L \}$ to refer to the state of all of the \emph{intermediate} neurons (i.e., excluding the output), and we denote $z^\ell_i$ to be the $i^\text{th}$ element of  $\mathbf{z}^\ell$. 
The goal of the GNCN model is to learn a joint distribution over its $L+1$ neural states, i.e., $p(\mathbf{z}^0,\mathbf{z}^1,\cdots,\mathbf{z}^L)$, from which the marginal distribution of the data may be obtained via $p(\mathbf{x}) = \int_\mathbf{Z} p(\mathbf{x},\mathbf{z}^1,\cdots,\mathbf{z}^L) d \mathbf{Z}=\int_\mathbf{Z} p(\mathbf{z}^0,\mathbf{z}^1,\cdots,\mathbf{z}^L) d \mathbf{Z}$. The model decomposes the marginal probability of its internal state  as a product of conditional probabilities and thus models the joint distribution of data and internal states as $p(\mathbf{x},\mathbf{z}^1,\cdots,\mathbf{z}^L) = p(\mathbf{z}^0|\mathbf{z}^1) \cdots p(\mathbf{z}^{L-1}|\mathbf{z}^L)p(\mathbf{z}^L)$.
For (binary) image sensory patterns, we set $p(\mathbf{z}^0|\mathbf{z}^1)$ to be a multivariate Bernoulli distribution with mean vector denoted by the symbol $\mean{z}^0$. This mean vector is computed from the model using the state  $\mathbf{z}^1$ of the previous layer. For the non-output layers, the conditional probabilities $p(\mathbf{z}^\ell|\mathbf{z}^{\ell+1})$ are modelled as multivariate Gaussians with mean vector denoted by the symbol $\mean{z}^\ell$ (which is computed by the model from the state $\mathbf{z}^{\ell+1}$ of the previous layer). The covariance matrix $\mathbf{\Sigma}^\ell$ is a fixed, trainable parameter for each layer. 

\begin{figure}[t]
\begin{center}
\begin{subfigure}{0.5\textwidth}
  \centering
 
  \includegraphics[width=4.7cm]{figs/ncn_process.png}  
 
  \label{fig:inhibit}
\end{subfigure}%
\begin{subfigure}{0.5\textwidth}
  \centering
 
  \includegraphics[width=3.8cm]{figs/ncn_circuit.png}  
 
  \label{fig:predict}
\end{subfigure}
\end{center}
\vspace{-0.25cm}
\caption{(Left) The two key computation steps taken by an entire NGC network (a GNCN) when processing an input ($\mathbf{z}^0 = \mathbf{x}$): (1) prediction and laterally-weighted error computation, (2) error-correction of neural states. In this diagram, we depict a toy network with $3$ layers of $2$ state neurons (grey circles), i.e., $\mathbf{z}^2 = ([z^2_0, z^2_1 ])^T$, $\mathbf{z}^1 = ([z^1_0, z^1_1 ])^T$, $\mathbf{z}^0 = ([z^0_0, z^0_1 ])^T$, that are updated iteratively over $T$ time steps. Two of these layers are linked to error neurons (green diamonds), i.e., $\mathbf{e}^1 = ([e^1_0, e^1_1 ])^T$, $\mathbf{e}^0 = ([e^0_0, e^0_1 ])^T$, which compute mismatch messages that are propagated throughout the system. The bottom layer $\mathbf{z}^0$ receives sensory input, i.e., an image.
(Right) The basic neural computational unit that a GNCN is composed of, consisting of a single state neuron $z^\ell_i$ and an error neuron $e^\ell_i$ at layer $\ell$. In the circuit, observe that a state neuron not only receives messages from error neurons (carried by $E^{\ell-1}$ synapses) but also a self-excitation signal and inhibition signals from laterally connected neurons (via $V^\ell$ synapses). The error neuron receives a gain signal (via $(\Sigma^\ell)^{-1}$ synapses) from laterally connected error neurons.  Note that filled diamonds indicate inhibitory signals, non-filled circles indicate excitatory signals, and empty squares indicate multiplicative signals.}
\label{fig:ngc_process}
\vspace{-0.5cm}
\end{figure}

Note that the mean vector $\mean{z}^\ell$ depends on the sampled realization $\mathbf{z}^{\ell+1}$ of the layer state above, meaning that the GNCN could be viewed as a hierarchical continuous Gaussian mixture model. 
Furthermore, key to the learning and inference process, we introduce a vector of error neurons $\mathbf{e}^\ell$ (depicted in Figure \ref{fig:bp_v_ngc} Right) that are tasked with computing how far off the mean vector is from the relevant nearby state.
Each entry  $\mean{z}^\ell_i$ of the mean vector and  $\mathbf{e}^\ell_i$ of the error neurons is computed by the following:
\begin{align}
   \mean{z}^\ell_i =  \overbrace{ g^\ell\Big( \sum_{j \in J_{\ell+1}} \mathbf{W}^{\ell+1}_{ij} \phi^{\ell+1}( \mathbf{z}^{\ell+1}_j ) \Big) }^\text{local, top-down prediction}, \quad \mbox{ and,} \quad \mathbf{e}^\ell_i = \sum_{j \in J_\ell} \overbrace{ (\mathbf{\Sigma}^\ell)^{-1}_{ij} }^\text{lateral modulation} ( \mathbf{z}^\ell_j - \mean{z}^\ell_j ) \label{eqn:error_neuron}
\end{align}
where $J_\ell$ is the number of neurons in layer $\ell$, $\mathbf{W}^\ell$ is the weight matrix and $g^\ell$ and $\phi^{\ell+1}$ are activation functions for layer $\ell$.
Thus each layer $\mathfrak{N}^\ell$ is specified by two functions $g^\ell$ and $\phi^\ell$, a trainable weight matrix $\mathbf{W}^\ell$ and a covariance matrix $\mathbf{\Sigma}^\ell$.
In Figure \ref{fig:ngc_process} (Right), we depict the circuit for a single pair neurons at a layer $\ell$, i.e., a state neuron $\mathbf{z}^\ell_i$ and an error neuron $\mathbf{e}^\ell_i$.
Notice that in Equation \ref{eqn:error_neuron} above, that in order to compute the state of the error neuron $e^\ell_i$, the covariance matrix $\mathbf{\Sigma}^\ell$ acts as a lateral modulation matrix, which is inspired by the neuro-mechanistic concept of precision weighting in predictive processing theory \cite{moran2013free}. It allows error neuron $e^\ell_i$ to dynamically amplify/reduce the learning signal (i.e., $z^\ell_i-\mean{z}^\ell_i$) of its corresponding state neuron $z^\ell_i$, based on the learning signals of the other state neurons.
Empirically, we found these modulatory synapses to improve the crispness of model samples. Note that Equation \ref{eqn:error_neuron} would be applied to all $J_\ell$ error neurons for each layer $\ell = 0, 1, \cdots, L-1$.

Once the $i$th error neuron at layer $\ell$ has been calculated, the corresponding state neuron will proceed to update its state. By analogy to Gaussian mixture models, where the state neurons $\mathbf{z}^\ell$ correspond to latent variables, this can be viewed as an attempt to modify the states in a way that improves the complete data log-likelihood (i.e. modifying the $\mathbf{z}^\ell$ to cause  $p(\mathbf{x},\mathbf{z}^1,\cdots,\mathbf{z}^L) $ to increase; see Appendix for full definition of log likelihood $\psi$).
One possible neuroscience-inspired way to perform this update is shown in Equation \ref{eqn:naive_state} (here $\beta$ is a hyperparameter akin to machine learning concept of a learning rate). Specifically, this update is:
\begin{align}
    \mathbf{z}^\ell_i \leftarrow \mathbf{z}^\ell_i + \beta \Big( - \gamma \mathbf{z}^\ell_i  -  \mathbf{e}^\ell_i + \sum_{j \in J_{\ell-1}} ( \mathbf{E}^{\ell-1}_{ij} \mathbf{e}^{\ell-1}_j )  \Big) \mbox{.} \label{eqn:naive_state}
\end{align}
Here $z^\ell_i$ is modified through three terms. The first is a decaying pressure caused by the leak term $-z^\ell_i$, controlled by the strength factor $\gamma$. The second term, $-e^\ell_i$, can be interpreted as top-down pressure where $e^\ell_i$ is a measure of how much the neuron's state differs from the predicted state $\mean{z}^\ell$ that is computed by the layer above. Finally, the third term 
adds in the error message from each error neuron $e^{\ell-1}_j$ in the layer below, communicated by special error synapses $E^{\ell-1}$, acting as a form of bottom-up pressure.

However, instead of directly using this update rule, we further incorporate another property of the brain -- activation sparsity. Sparsity in real neuronal and artificial systems is often argued to be useful in not only learning compact representations (since most activitie value will be suppressed or will be zero), which allow for efficient storage and vastly improved energy efficiency, but also allows the simultaneous representation of distinct patterns with little interference \cite{ahmad2016neurons} while still maintaining a large representational capacity.
To emulate this type of sparsity, we integrate a mechanism to force neurons to compete for activation, where we take inspiration from the known occurrence of lateral synapses in cortical regions of the brain (which are often thought to facilitate contextual processing and determining the properties of neural receptive fields \cite{spratling2001dendritic}).
To do this, we introduce two terms to Equation \ref{eqn:naive_state} that use excitatory/inhibitory synapses stored in a matrix $\mathbf{V}^\ell$. This means that state neurons are updated as follows:
\begin{align}
    \mathbf{z}^\ell_i &\leftarrow \mathbf{z}^\ell_i + \beta \Big( \overbrace{-\gamma \mathbf{z}^\ell_i}^\text{leak} + \overbrace{ \Big( \sum_{j \in J_{\ell-1}} \mathbf{E}^{\ell-1}_{ij} \mathbf{e}^{\ell-1}_j \Big) - \mathbf{e}^\ell_i }^\text{bottom-up + top-down pressures} - \overbrace{ \Big( \sum_{j \in J_\ell, j \neq i} \mathbf{V}^\ell_{ij} \phi^\ell(\mathbf{z}^\ell_j) \Big) }^\text{lateral inhibition} + \overbrace{\mathbf{V}^\ell_{ii} \phi^\ell(\mathbf{z}^\ell_i)}^\text{self-excitation} \Big)  \label{eqn:state_unit}
\end{align}
Depending on the values set in $\mathbf{V}^\ell$, different types of sparsity patterns emerge, creating a flexible means for testing the benefits/drawbacks of different kinds of lateral competition patterns in an interpretable manner.
Figure \ref{fig:lateral_design} provides a graphical example of the type of interaction pattern we found worked well for the GNCN in this study, i.e., we forced $J_\ell/K$ groups (or columns) of $K$ neurons to compete with each other (see Appendix). 

\begin{figure}[!t]
\begin{center}
\includegraphics[width=6.25cm]{figs/lateral_design.png}  
\end{center}
\vspace{-0.25cm}
\caption{Depending on how the lateral connectivity matrices, $\mathbf{V}^1$ \& $\mathbf{V}^2$, are designed, different competition patterns emerge among neurons in network. In the system with lateral matrices as shown above, each neuron is driven by its own self-excitation (red-colored blocks) and laterally inhibits (blue-colored blocks) other neurons that inhabit the same group (of $2$ or $3$ units). In the figure, one possible outcome (or state of the agent) of such a competition is shown (after $T$ steps). Self-excitation and lateral inhibition strengths are controlled by coefficients set a priori. Dashed orange edges show actively used synapses, filled diamonds indicate inhibitory signals, and non-filled circles indicate excitatory signals.} 
\label{fig:lateral_design}
\vspace{-0.35cm}
\end{figure}

Putting the ideas above together, we can characterize the neural dynamics of an NGC system (graphically depicted in Figure \ref{fig:ngc_process} Left). When presented with sensory input $\mathbf{x}$, the GNCN sets $\mathbf{z}^0=\mathbf{x}$ and applies Equations \ref{eqn:error_neuron} and \ref{eqn:state_unit} to each layer $L,\dots, 1$ in the system (cycling through the levels $T$ times). 
Once the activity values for $\mathbf{z}^\ell$ and $\mathbf{e}^\ell_i$ have been updated at the end of a $T$-step episode, the synaptic weights $\mathbf{W}^\ell$ are adjusted by computing products between the relevant activities and error neurons (see Appendix), following an error Hebbian-like update. For example, the update to any single synapse $\mathbf{W}^\ell_{ji}$ is simply $\Delta \mathbf{W}^\ell_{ij} = (\mathbf{e}^\ell_i) \phi(\mathbf{z}^{\ell+1}_j)$. Desirably, in contrast to backprop, no activation function derivatives are needed.
Note that the weight and the neural state updates come from the objective that the GNCN optimizes online, i.e., total discrepancy (see Appendix). 
In essence, the steps described above and shown in Figure \ref{fig:ngc_process} (Left) illustrate that the NGC framework embodies the idea that neural state and synaptic weight adjustment are the result of a process of generate-then-correct, or \emph{continual error correction}, in response to samples of the agent's environment.


\textbf{Generation Results:} The model framework that we have described so far would be immediately useful for creating an auto-associative memory of sensory input, i.e., upon receiving a particular sensory input, the model would be able to recall seeing it by accurately reconstructing it. Furthermore, since the model is learning an estimator of the input distribution $p(\mathbf{x})$, we may sample from it (see Methods) to synthesize or ``hallucinate'' data patterns, as we will demonstrate later in this section

\begin{table}[!t]
\begin{center}
\caption{Generative modeling results across datasets (averaged over $10$ trials). The binary cross entropy (BCE) of the model reconstructions and the marginal log likelihood $\log p(\mathbf{x})$ (in nats) on the test set are reported (lower BCE is better and $\log p(\mathbf{x})$ closer to $0.0$ is better). We \textbf{bold} the best two scores of the models with respect to BCE (first column) as well as with respect to $\log p(\mathbf{x})$ (second column).
}
\label{results:gen_results}
 \begin{tabular}{l | l c c | l c c} 
 \hline
 
 \textbf{Model} & & \multicolumn{1}{c}{\textbf{BCE}} & \multicolumn{1}{c}{$\log p(\mathbf{x})$}  & & \multicolumn{1}{c}{\textbf{BCE}} & \multicolumn{1}{c}{$\log p(\mathbf{x})$} \\
 \hline\hline
  GMM & \parbox[t]{0.15mm}{\multirow{7}{*}{\rotatebox[origin=c]{90}{MNIST}}} & -- & $-185.37 \pm 0.47$ &  \parbox[t]{0.15mm}{\multirow{7}{*}{\rotatebox[origin=c]{90}{KMNIST}}} & -- & $-426.73 \pm 1.03$ \\
 RAE & & $\mathbf{52.55 \pm 0.11}$ & $-119.38 \pm 1.16$ &  & $\mathbf{165.20 \pm 0.40}$ & $-240.40 \pm 0.85$ \\
 GVAE-CV & & $71.19 \pm 0.40$ & $-104.93 \pm 0.38$ & & $182.05 \pm 1.06$ & $-229.67 \pm 0.83$ \\
 GVAE & & $75.95 \pm 0.96$ & $-98.87 \pm 0.75$ & & $198.28 \pm 1.89$ & $\mathbf{-218.3 \pm 0.88}$ \\
 GAN-AE & & $99.30 \pm 0.14$ & $\mathbf{-95.50 \pm 0.17}$ & & $230.64 \pm 1.57$ & $\mathbf{-219.66 \pm 0.58}$ \\
 GNCN & & $\mathbf{58.05 \pm 0.11}$ & $\mathbf{-98.29 \pm 0.12}$ & & $\mathbf{142.22 \pm 0.23}$ & $-222.76 \pm 0.30$ \\
 \hline
 \hline
 GMM & \parbox[t]{0.15mm}{\multirow{7}{*}{\rotatebox[origin=c]{90}{FMNIST}}} & -- & $-302.96 \pm 0.47$ & \parbox[t]{0.15mm}{\multirow{7}{*}{\rotatebox[origin=c]{90}{CalTech}}} & -- & $-88.63 \pm 0.43$  \\
 RAE & & $\mathbf{99.94 \pm 0.35}$ & $-137.26 \pm 0.50$  & & $\mathbf{43.01 \pm 0.51}$ & $-45.09 \pm 0.22$ \\
 GVAE-CV & & $114.87 \pm 1.01$ &  $-131.4 \pm 0.43$ &  & $52.69 \pm 0.91$ & $\mathbf{-41.62 \pm 0.26}$ \\
 GVAE & & $121.98 \pm 1.41$ & $\mathbf{-127.67 \pm 0.39}$  & & $53.91 \pm 1.40$ & $-42.91 \pm 0.54$ \\
 GAN-AE & & $138.36 \pm 1.31$ & $-136.14 \pm 0.48$ & & $59.39 \pm 2.67$ & $-43.79 \pm 1.54$ \\
 GNCN & & $\mathbf{99.79 \pm 0.24}$ & $\mathbf{-131.32 \pm 0.29}$ & & $\mathbf{36.84 \pm 0.09}$ &  $\mathbf{-40.95 \pm 0.15}$ \\ %
 \hline
\end{tabular
\end{center}
\vspace{-0.5cm}
\end{table}

To evaluate our framework, six approaches were compared across four image datasets (Table \ref{results:gen_results}). One is a Gaussian mixture model (GMM) and five are neural models -- four of these are backprop-based and one is NGC-based. All five neural models were constrained to have their top-most layer to contain $20$ processing elements (in the case of the GNCN, $20$ neural columns) and all had approximately the same total number of synapses.
The regularized auto-encoder (RAE), the auto-associative network least equipped to serve as a data synthesizer, reaches lower reconstruction error, in terms of binary cross entropy (BCE), compared to the other backprop-based generative models, i.e., the Gaussian variational autoencoder (GVAE), the constant-variance GVAE (CV-GVAE), and the adversarial autoencoder (GAN-AE). However, while the CV-GVAE, GVAE, and GAN-AE models yield worse reconstruction than the RAE, they obtain much better log likelihoods, especially compared to the GMM baseline, indicating that they strong data samplers. It makes sense that these models obtain better likelihood at the expense of BCE given that their optimization objective imposes a strong pressure to craft a proper (Gaussian) distribution over latent variables in addition to reconstructing data samples (in the case of the GAN-AE, the pressure comes from forcing the discriminator to distinguish between fake and real latent variables).
Interestingly enough, we see that \textbf{the GNCN obtains competitive log likelihood with the CV-GVAE, GVAE, and GAN-AE while yielding the best reconstruction for all four datasets}, even better than the RAE itself. On MNIST, we note that the GNCN outperforms some other related prior models \cite{bengio2013generalized} that used our same experimental setup, e.g., a restricted Boltzmann machine with $\log p(\mathbf{x}) = -112$ nats
, a denoising autoencoder with $\log p(\mathbf{x}) = -142$ nats
(trained via backprop) and $\log p(\mathbf{x}) = -116$ nats (using the walk-back algorithm). As indicated by our results, the GNCN's sampling and reconstruction ability is very promising

As shown in Figure \ref{fig:convergence}, where the final model test BCE was measured throughout training, the GNCN consistently generalizes sooner than all of the backprop-based models (similar loss curve patterns were observed on the training set). Even after only passing through the dataset once, the GNCN reaches much lower BCE than the others. For example, on MNIST, the GNCN reached a $\mbox{BCE} = 57$ nats by the end of the first pass while the GVAE reached $\mbox{BCE} = 194$ nats and the RAE reached $\mbox{BCE} = 168$ nats.
Furthermore, in Figure \ref{fig:gncn_samples}, we present samples obtained from either: 1) randomly sampling the original dataset, 2) ancestrally sampling the GVAE, or 3) ancestrally sampling (see Methods) the GNCN.
Observe that both the GNCN and GVAE yield reasonably good-looking sample images for all four datasets and the differences in perceptual quality between the models' sets of samples is marginal. This is encouraging, given that our goal was to demonstrate that the GNCN could be competitive with backprop-based generative models.

\begin{figure}[!t]
\begin{center}
\begin{subfigure}{0.25\textwidth}
  \centering
  \includegraphics[width=4.25cm]{figs/mnist_converge.png}  
  \caption{MNIST.}
  \label{fig:mnist_converge}
\end{subfigure}%
\begin{subfigure}{0.25\textwidth}
  \centering
  \includegraphics[width=4.25cm]{figs/kmnist_converge.png}  
  \caption{KMNIST.}
  \label{fig:kmnist_converge}
\end{subfigure}%
\begin{subfigure}{0.25\textwidth}
  \centering
  \includegraphics[width=4.25cm]{figs/fmnist_converge.png}  
  \caption{FMNIST.}
  \label{fig:fmnist_converge}
\end{subfigure}%
\begin{subfigure}{0.25\textwidth}
  \centering
  \includegraphics[width=4.25cm]{figs/ctech_silh_converge.png}
  \caption{CalTech.}
  \label{fig:ctech_converge}
\end{subfigure}
\end{center}
\caption{BCE loss curves (lower is better) for the pattern reconstruction task on each dataset (averaged over 10 trials).} 
\label{fig:convergence}
\vspace{-0.35cm}
\end{figure}

\begin{figure}[!t]
\centering
\minipage{0.165\textwidth}%
  \includegraphics[width=\linewidth]{figs/mnist_data_patterns.jpg}\\
  \vspace{-0.26cm}
  \includegraphics[width=\linewidth]{figs/gvae_mnist_samples.jpg}\\
  \vspace{-0.26cm}
  \includegraphics[width=\linewidth]{figs/gncn_mnist_samples.jpg}
 
\endminipage\hspace{0.185cm}
\minipage{0.165\textwidth}%
  \includegraphics[width=\linewidth]{figs/kmnist_data_patterns.jpg}\\
  \vspace{-0.26cm}
  \includegraphics[width=\linewidth]{figs/gvae_kmnist_samples.jpg}\\
  \vspace{-0.26cm}
  \includegraphics[width=\linewidth]{figs/gncn_kmnist_samples.jpg}
 
\endminipage\hspace{0.185cm}
\minipage{0.165\textwidth}%
  \includegraphics[width=\linewidth]{figs/fmnist_data_patterns.jpg}\\
  \vspace{-0.26cm}
  \includegraphics[width=\linewidth]{figs/gvae_fmnist_samples.jpg}\\
  \vspace{-0.26cm}
  \includegraphics[width=\linewidth]{figs/gncn_fmnist_samples.jpg}
 
\endminipage\hspace{0.185cm}
\minipage{0.165\textwidth}%
  \includegraphics[width=\linewidth]{figs/ctech_silh_data_patterns.jpg}\\
  \vspace{-0.26cm}
  \includegraphics[width=\linewidth]{figs/gvae_ctech_samples.jpg}\\
  \vspace{-0.26cm}
  \includegraphics[width=\linewidth]{figs/gncn_ctech_samples.jpg}
 
\endminipage
\caption{Visualization of samples from models with best BCE-$\log p(\mathbf{x})$ balance. Top row shows dataset samples, middle row shows GVAE samples, and bottom row shows GNCN samples. Columns are arranged, left to right, as: a) MNIST, b) KMNIST, c) FMNIST, and d) CalTech.}
\label{fig:gncn_samples}
\vspace{-0.55cm}
\end{figure}

\subsection{Neural Generative Coding Yields Strong Downstream Pattern Classifiers}
\label{sec:results_classify}

All of the generative models we have experimented with in this paper are unsupervised in nature, meaning that by attempting to learn a density estimator of the data's underlying distribution, the representations acquired by each might prove useful for downstream applications, such as image categorization. To evaluate each how useful each model's latent representations might be when attempting to discriminate between samples, we evaluate the performance of a simple log-linear classifier, i.e., maximum entropy, that is fit to each model's topmost latent variable using the labels accompanying each dataset. 
For all models, we measure the classification error (\textbf{Err}), as a percentage, on the test set of each benchmark in Table \ref{results:diagonistics}, where the closer a model is to $0\%$, the better.
In addition, we provide the results for simple, purely discriminative baseline for context (the DSRN), which is simply a backprop-trained sparse recitfier network \cite{glorot2011deep} that is constrained to have the same number of synapses as the generative models. As we see in Table \ref{results:diagonistics}, in terms of test error, \textbf{the GNCN is competitive with the purely discriminative DSRN and outperforms all of the other generative models} (though it does not outperform the DSRN in two out of four cases).

In Figure \ref{fig:tsne_latents}, we provide qualitative evidence that the latent representations of the GNCN appear to yield a stronger, natural separation of the test data points into seemingly class-respective clusters as compared to the GVAE (one of the best performing backprop-based models with respect to both log likelihood and reconstruction error). Again, we emphasize that the GNCN acquired these representations without labeled information, meaning that the class-based relationships have emerged as a result of its very sparse neural activities.  This offers some promising evidence that the GNCN's representations offer benefits beyond the original density estimation task, allowing reuse of the same system for downstream tasks like image categorization. 

\begin{table}[t]
\begin{center}
\caption{(Left) Downstream diagnostics across datasets. The classification error (Err) of a log-linear model fit to each model's latent codes and the masked mean squared error (M-MSE) of the model's pattern completion ability on the test set are reported (lower is better -- two best scores are in \textbf{bold} with respect to each metric/column). (Right) Masked patterns in first row and pattern completions from the GNCN in the second row. }
\label{results:diagonistics}
 \begin{tabular}{l | l c c | l c c } 
 \hline
 
 \textbf{Model} & & \multicolumn{1}{c}{\textbf{Err} (\%)} & \multicolumn{1}{c}{\textbf{M-MSE}} &  & \multicolumn{1}{c}{\textbf{Err} (\%)} & \multicolumn{1}{c}{\textbf{M-MSE}}  \\
 \hline\hline

 DSRN & \parbox[t]{0.15mm}{\multirow{7}{*}{\rotatebox[origin=c]{90}{MNIST}}} & $\mathbf{1.93 \pm 0.04}$ & -- & \parbox[t]{0.15mm}{\multirow{7}{*}{\rotatebox[origin=c]{90}{KMNIST}}} & $\mathbf{10.02 \pm 0.08}$ & -- \\

 RAE & & $10.55 \pm 0.31$ & $25.17 \pm 0.61$ & & $37.36 \pm 1.16$ & $32.47 \pm 0.32$ \\
 GVAE-CV & & $13.84 \pm 0.75$ & $20.20 \pm 0.46$ & & $36.64 \pm 0.89$ & $30.20 \pm 0.47$ \\
 GVAE & & $7.47 \pm 0.59$ & $19.81 \pm 0.31$ & & $43.93 \pm 1.60$ & $28.91 \pm 0.37$ \\
 GAN-AE & & $9.68 \pm 0.50$ & $\mathbf{16.22 \pm 0.15}$ & & $35.12 \pm 0.80$ & $\mathbf{27.15 \pm 0.33}$ \\
 GNCN & & $\mathbf{1.32 \pm 0.23}$ & $\mathbf{14.19 \pm 0.11}$ & & $\mathbf{12.67 \pm 0.27}$ & $\mathbf{26.90 \pm 0.11}$ \\

 \hline
 \multicolumn{5}{c}{ } \\

 \hline

 DSRN & \parbox[t]{0.15mm}{\multirow{7}{*}{\rotatebox[origin=c]{90}{FMNIST}}} & $\mathbf{9.75 \pm 0.09}$ & -- & \parbox[t]{0.15mm}{\multirow{7}{*}{\rotatebox[origin=c]{90}{CalTech}}} & $\mathbf{31.56 \pm 0.01}$ & -- \\

 RAE & & $21.09 \pm 0.52$ & $54.20 \pm 2.92$ & & $37.42 \pm 0.58$ & $23.56 \pm 0.91$ \\
 GVAE-CV & & $19.95 \pm 0.26$ & $53.05 \pm 2.84$ & & $39.57 \pm 0.42$ & $27.22 \pm 1.43$ \\
 GVAE & &  $26.44 \pm 0.86$ & $50.72 \pm 3.28$ & & $42.6 \pm 0.63$ & $24.08 \pm 2.31$ \\
 GAN-AE & & $20.45 \pm 0.79$ & $\mathbf{38.79 \pm 0.90}$ & & $39.28 \pm 0.92$ & $\mathbf{19.47 \pm 0.88}$ \\
 GNCN & & $\mathbf{14.84 \pm 0.44}$ & $\mathbf{25.33 \pm 1.48}$ & & $\mathbf{29.88 \pm 0.62}$ & $\mathbf{9.60 \pm 0.48}$ \\

 \hline
\end{tabular
\begin{tabular}{cccc}
\hline
\multicolumn{2}{c}{\hspace{0.25cm}\textbf{MNIST}} & \multicolumn{2}{c}{\textbf{KMNIST}} \\
 \includegraphics[width=0.055\linewidth]{mnist_mask/samp1_o.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{mnist_mask/samp2_o.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{kmnist_mask/samp1_o.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{kmnist_mask/samp2_o.jpg} \\
 \includegraphics[width=0.055\linewidth]{mnist_mask/samp1_r.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{mnist_mask/samp2_r.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{kmnist_mask/samp1_r.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{kmnist_mask/samp2_r.jpg} \\
 \hline
 \multicolumn{4}{c}{} \\
 \hline 
 \multicolumn{2}{c}{\hspace{0.25cm}\textbf{FMNIST}} & \multicolumn{2}{c}{\textbf{CalTech}} \\
 \includegraphics[width=0.055\linewidth]{fmnist_mask/samp1_o.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{fmnist_mask/samp2_o.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{ctech_mask/samp1_o.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{ctech_mask/samp2_o.jpg} \\
 \includegraphics[width=0.055\linewidth]{fmnist_mask/samp1_r.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{fmnist_mask/samp2_r.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{ctech_mask/samp1_r.jpg}\hspace{-0.35cm} & \includegraphics[width=0.055\linewidth]{ctech_mask/samp2_r.jpg} \\
 \hline 
\end{tabular}
\end{center}
\vspace{-0.55cm}
\end{table}

\begin{figure}[!tbp]
\centering
\begin{tabular}{c@{\hskip -0.8cm}c@{\hskip -0.4cm}}
\includegraphics[width=5.0cm]{figs/gvae_mnist_latents.jpg} &
\includegraphics[width=5.0cm]{figs/gncn_mnist_latents.jpg}
\end{tabular}%
\begin{tabular}{ c|c|c|c } 
 \textit{GNCN} &  $\rho_a^1$ &  $\rho_a^2$ &  $\rho_a^3$ \\
 \hline
 \textbf{MNIST} & $0.16$ & $0.13$ & $0.11$\\ 
 \textbf{KMNIST} & $0.21$ & $0.17$ & $0.15$\\ 
 \textbf{FMNIST} & $0.18$ & $0.16$ & $0.15$ \\
 \textbf{CalTech} & $0.18$ & $0.14$ & $0.13$ \\   
 \hline
\end{tabular}
\vspace{-0.45cm}
\caption{(Left) t-SNE visualization of latent codes of the GVAE \& GNCN. 
(Right) Sparsity levels in the GNCN, $\rho^\ell_a$ indicates sparsity for layer $\ell$.}
\label{fig:tsne_latents}
\vspace{-0.685cm}
\end{figure}

We hypothesize that the GNCN's latent codes make it easier for a linear classifier to separate out patterns by category as a result of the fact that they are sparse. Crucially, the lateral structure of our model creates columns of neural processing units that fight for the right to represent and explain the input, meaning that only a few (possibly even just one) within a group would be active when processing particular patterns.
To quantity this sparsity in the GNCN, we measure and report in the table of Figure \ref{fig:tsne_latents} (Right) the (mean) proportion $\rho^\ell_{a}$ of neurons at layer $\ell$ that were active (we counted whether each unit for a given pattern satisfied $z^\ell_i > \epsilon$, where $\epsilon = 1e^{-6}$)) in the model for each dataset. 
Note that the activities of the GNCN are quite sparse (the greatest value found was on KMNIST in layer $1$ with a sparsity of $21$\%) and, furthermore, that the number of active neurons is lower for deeper layers, with $11$-$15$\% sparsity for all datasets for the top layer $\mathbf{z}^3$.
We believe that the GNCN's structured form of sparsity aids it in modeling the input, i.e., yielding strong log likelihood and reconstruction error, while facilitating better separation of highly nonlinear data. Though models like the GVAE and CV-GVAE reach good log likelihoods too, it is possible that since their learning process to focus on shaping their latent spaces to dense, multivariate Gaussian distributions, there is little chance for useful sparsity to emerge as a by-product, reducing the possibility that the latent space might result in beneficial side-effects.

From a neurobiological perspective, it is well-known that lateral synaptic connections often facilitate contextual processing \cite{cui2017sparseHTM} and offer a natural form of ``activity sharpening'' \cite{adesnik2010lateral}.
We believe that this is an important structural prior that biases the GNCN towards acquiring more economical representations \cite{barlow1972single}, where progressively fewer neurons at layers progressively farther away from the input stimulus work to encode information (or are non-zero).
Much like the more recent incarnations of spike-and-slab coding \cite{goodfellow2013spikeslab}, our form of structurally-enforced sparsity has a much stronger regularization effect on the model and ensures that the generative distribution of the neural activities is truly sparse. This is unlike more traditional approaches where sparsity is weakly enforced through the use of factorial kurtotic priors applied during inference  \cite{olshausen1997sparse,rao1999predictive}.
Our NGC frameowrk not only yields naturally sparse codes but also offers flexibility in exploring other types of lateral connectivity patterns beyond the choice made in this study

\subsection{Neural Generative Coding Can Conduct Pattern Completion}
\label{sec:results_completion}
Another interesting ability attributed to auto-associative memory models is their ability to complete partially-corrupted or incomplete patterns. In the real world, this type of scenario would often occur in the form of object occlusion, where the view of an object might be partially obstructed from the agent's view. Being able to ``imagine'' the rest of the object might prove useful when planning to grasp it or manipulate it in some fashion.
To test each model's ability to complete patterns, we conducted an experiment where the right half of each image in each dataset was masked and each model was tasked with predicting the deleted portions. In Table \ref{results:diagonistics}, we report the masked mean squared error (M-MSE, see Methods) of each model on each dataset's test set.

Interestingly enough, we see that \textbf{the GNCN significantly outperforms all of the other baselines in terms of pattern completion}, with respect to measured M-MSE. We furthermore provide some examples of original data patterns that were masked in Table \ref{results:diagonistics} (in the first and third rows) contrasted with the GNCN's completed patterns  (in the second and fourth rows). 
Much like the analysis on downstream classification, we hypothesize that the GNCN is able to conduct better pattern completion because it is not as constrained to learn an actual distribution over its top-most hidden layer and thus able to extract sparse latent codes given its lateral synapses.


\section{Discussion}
\label{sec:discussion}
Generative models based on artificial neural networks (ANNs) have yielded promising tools for estimating and sampling from complicated probability data distributions \cite{kingma2013auto,goodfellow2014generative}. By re-considering how these models might operate and learn, drawing inspiration from one promising neuro-mechanistic account of how the brain interacts and adapts to its environment, i.e., predictive processing \cite{clark2015surfing},
we have shown that learning a viable generative model is possible. Specifically, we propose the neural generative coding (NGC) computational framework for learning neural probabilistic models of data, implementing one concrete instantiation of NGC which we called the generative neural coding network (GNCN). 
Surprisingly, the GNCN is not only competitive with several powerful, modern-day backprop-based models on the task of estimating the marginal distribution of the data but it can generalize beyond the task it was trained to do. Specifically, we investigated the performance of the GNCN on downstream tasks such as pattern completion and pattern categorization and discovered that our unsupervised model outperformed all of the examined backprop-based baselines and was even competitive with an ANN that was directly trained to specialize for classification.
As a result, for systems as complex as probabilistic generative models, we have demonstrated that crafting a more fundamentally brain-inspired approach to information processing and credit assignment can yield artificial neural systems that extract rich representations of input data in an unsupervised fashion.
Our results demonstrate, on four datasets, that even though extra computation is needed to process each input for a fixed (yet small) stimulus presentation time, the GNCN converges far sooner than comparable backprop-based ones, generalizing well earlier.

In recent research, developing learning procedures that enable backprop-level learning while embodying elements of actual neuronal function has seen increasing interest in the machine learning community. However, while insights provided by each development have proven valuable, increasing the evidence that shows how a backprop-free form of adaptation can be consistent with some aspects  real networks of neurons, many of these ideas only address one or a few of the issues described earlier.
Random feedback alignment algorithms \cite{baldi2016learning,lillicrap2016random} address the weight transport problem, and to varying degrees, the update locking problem \cite{nokland2016direct,mostafa2018deep,frenkel2019learning}, but are fundamentally emulating backprop's differentiable global feedback pathway to create teaching signals (and pay a reduction in generalization ability the farther away they deviate from backprop \cite{bartunov2018assessing,frenkel2019learning}). In addition, experimentally, the success of these approaches depends on how well the feedback weights are chosen a prior (instead of learning them). Other procedures, like local representation alignment \cite{ororbia2019biologically} and target propagation \cite{hinton1988learning,lee2015difference}, which also resolve the weight transport problem and eshew the need for differentiable activations, fail to address the update locking problem, since the various incarnations of these require a full forward pass to initiate inference. 
Procedures have been proposed to address the update locking problem, such as the method of synthetic gradients \cite{jaderberg2017decoupled}, but still require backprop to compute local gradients. It remains to be seen how these algorithms could be adapted to learn generative models.
Other systems, such as those related to contrastive Hebbian learning (CHL)  \cite{movellan1991contrastive,o1996biologically,scellier2018generalization},
are much more biologically-plausible but often require symmetry between the forward and backward synaptic pathways, i.e., failing to address weight transport. More importantly, CHL requires long settling phases in order to compute activities and teaching signals, resulting in long computational simulation times.
However, while some of these algorithms have shown some success in ANN training \cite{lee2015difference,lillicrap2016random,ororbia2019biologically}, they focus on classification, which is purely supervised and arguably a simpler problem than generative modeling. 

Boltzmann machines \cite{ackley1985learning}, which are generalizations of Hopfield networks \cite{little1974existence,hopfield1982neural} to incorporate latent variables, are a type of generative network that also contains lateral connection between neurons much as the models in our NGC framework do. While training for the original model was slow, a simplification was later made to omit lateral synapses, yielding a bipartite graphical model referred to as a harmonium, trained by contrastive divergence \cite{hinton2002training}, a local contrastive Hebbian learning rule.
However, while powerful, the harmonium could only synthesize reasonable-looking samples with many iterations of block Gibbs sampling and the training algorithm suffered from mixing problems (leading low sample diversity among other issues) \cite{dumoulin2013challenges}.
Another kind of generative model \cite{hinton2005kind} can be trained with the wake-sleep algorithm
(or the up-down algorithm in the case of deep belief networks  \cite{hinton2005kind}), where an inference (upward) network and a generative (downward) network are jointly trained to invert each other. Unfortunately, wake-sleep suffers from instability, struggling to produce good samples of data most of the time, due to the difficulty both networks have in inverting each other due to layer-wise distributional shift.
Motivated by the deficiencies in models learned by contrastive divergence/wake-sleep, algorithms have been created for auto-encoder-based models \cite{bengio2013generalized} but most efforts today rely on backprop.

How the brain conducts credit assignment is a fundamental and open question in both (computational) neuroscience and cognitive science. There are many theories that posit how this might happen and our proposed NGC framework represents a scalable, computational instantiation of only one of them, the theory of predictive processing, suggesting that cortical regions in the brain communicate prediction errors and/or predictions across different regions through a hierarchical message passing scheme \cite{wacongne2011evidence}. Observe that the direction that this study takes is one that starts from cognitive neuroscientific concepts and ends in the development of a statistical learning algorithm. Specifically, we have shown that one way of emulating some neuro-biological principles yield agents that learn more general-purpose representations, as our results on downstream classification and pattern completion provide evidence for. As for implications for computational neuroscience itself, while the NGC/GNCN does embody concepts such as lateral competition-driven sparsity, hierarchical message passing, and local Hebbian-like synaptic adjustment, our framework lacks many important details that might allow it to serve as a means to make falsifiable claims that can be proven or disproven by neurobiological experiments. However, if the NGC framework is modified significantly to more faithfully model neurobiological details, e.g., synapses are constrained to take on only non-negative values and neurons communicate via spike trains (for example, the NGC could work with  leaky integrate-and-fire neurons, similar to that of \cite{ororbia2019spiking}), it could serve as the means to facilitate refinement of predictive brain theories such as predictive processing \cite{clark2015surfing} and principles such as free energy \cite{friston2006free}. Doing so could facilitate  stronger synergy between neural computational modeling and the design of agents that solve complex tasks examined in statistical learning. 

Given the difficulty in imagining how backprop, in the form it is implemented when training deep ANNs today, might occur in the brain \cite{grossberg1987competitive,crick1989recent}, there is value in not only developing approximations of it that might be more brain-like (moving machine learning a bit closer to computational neuroscience) but also in exploring alternatives that start from neuro-cognitive principles, theories, and mechanisms (creating new algorithms that embody particular ideas in cognitive neuroscience at the outset). Taking the second of the last two directions, which is in the spirit of this work and several others \cite{guerguiev2017towards,sacramento2018dendritic}, might allow us to more easily shed the constraints imposed by backprop in the effort to construct general-purpose learning agents capable of emulating more complex human cognitive function. Doing so might also allow the machine learning community to make further progress on problems even harder than generative modeling, such the problem of active inference \cite{friston2011action} and continual temporal prediction \cite{ororbia2020continual}

\section{Methods}
\label{sec:methods}


\paragraph{Training Setup:} 
The parameters of all models, whether they were updated via backprop or by the NGC learning process, were all optimized using stochastic gradient descent using mini-batches (or subsets, randomly sampled without replacement) of $200$ samples for $50$ passes through the data (epochs).
For the backprop-based models, we re-scaled the gradients \cite{pascanu2013difficulty} by re-projecting them to a Gaussian ball with radius of $5.0$, which we found ensured stable performance across trials.
For the GNCN models, after a weight update was made, the weight matrices were normalized such that the Euclidean norms of their columns were $1.0$
For each model, upon completion of training, we fit a Gaussian mixture model (GMM) to the topmost neural activity layer (to serve as the model prior). This density estimator was trained with expectation-maximization and contained $K = 65$ components, where each component $k$ defines a  multivariate Gaussian with mean $\mu_k$ and covariance $\Sigma_k$ parameters as well as a mixing coefficient $\pi_k$.

\paragraph{The Density Modeling Task:}
Given a dataset $\mathbf{X} \in \{0,1\}^{D \times S}$, where $S$ is the number of vector pattern samples and $D$ is the dimensionality of any given pattern, the goal is to learn a density model of $p(\mathbf{X})$, or the probability of sensory input data, where subset of $B$ vectors is denoted as $\mathbf{x} \in \{0,1\}^{D \times B}$. We parameterize the probability distribution $p(\mathbf{X})$ via $p_\Theta(\mathbf{X})$ by introducing learnable parameters $\Theta$. Since computing the marginal log likelihood $\log p_\Theta(\mathbf{X})$ directly is intractable for all the models in this paper, we estimate it by calculating a Monte Carlo estimate using $5000$ samples according to the recipe: 1) sample the GMM prior, 2) ancestrally sample directed neural generative model (whether a baseline or a GNCN) given the samples of the prior.

The reconstruction metric, binary cross-entropy (BCE), also known as the negative Bernoulli log likelihood, was computed as follows: $BCE(\mathbf{X}, \mathbf{\hat{X}}) = -\frac{1}{S} \sum^S_{s=1} \sum^D_{d=1} \Big( \mathbf{X} \log(\mathbf{\hat{X}}) + (1 - \mathbf{X}) \log(1 - \mathbf{\hat{X}}) \Big)[d,s] $ (measured in nats). $\mathbf{\hat{X}} \in \{0,1\}^{D \times S}$ is the predicted sensory input (matrix) produced from a model under evaluation.

\paragraph{The Pattern Completion Task:}
For this task, we test each model's ability to complete patterns where the images of each dataset were partially masked and each model was tasked with completing the masked images. Specifically, half of each image $\mathbf{x}$ (containing $\sqrt{D}$ columns of $\sqrt{D}$ pixels) was masked according to a binary mask $\mathbf{m}$ where $\sqrt{D}/2$ columns were set to $1$ and the rest $0$, i.e., $\mathbf{x}_m = \mathbf{x} \otimes \mathbf{m}$. We report the masked mean squared error (M-MSE) of each model on each dataset's test set, which is computed per image as follows: $\mbox{M-MSE}(\mathbf{X}, \mathbf{\hat{X}}, \mathbf{M}) = ((\mathbf{\hat{X}} - \mathbf{X}) \otimes (1 - \mathbf{M}))^T \cdot ((\mathbf{\hat{X}} - \mathbf{X}) \otimes (1 - \mathbf{M})) $ (measured in nats), where $\mathbf{M} \in \{0,1\}^{D \times S}$ is a masking matrix $\mathbf{m}$ where each column contains one mask vector per image vector $\mathbf{x}$.

\paragraph{The Classification Task:}
In this task, each model's prediction error is measured on the test sample's label set. This is possible given that each dataset $\mathbf{X}$ also comes with a set of target annotations that be encoded into a $C$-dimesnional space yielding the label matrix $\mathbf{Y} \in \mathcal{R}^{C \times D}$ (one label vector per sample), where $C$ is the number of unique categories labeled a priori.
We measure the classification error (\textbf{Err}) (as a percentage \%) on the test set, as follows:
\begin{align*}
    \mbox{Err}(\mathbf{Y},\mathbf{\hat{Y}}) = \bigg(1 - \frac{1}{S} \sum^S_{s=1} 
    \begin{array}{cc}
      \bigg\{ & 
        \begin{array}{cc}
          1 & \hat{y}_s = y_s \\
          0 & \hat{y}_s \neq y_s
        \end{array}
    \end{array}
    \bigg)
\end{align*}
where $\hat{y}_s = \arg \max_d (\mathbf{\hat{Y}}[:,s] )$, the class index chosen by the model, and $y_s = \arg \max_d (\mathbf{Y}[:,s] )$, the index of the actual class label. Note that $\mathbf{\hat{Y}} \in \mathcal{R}^{C \times D}$ is the collected set of predictions from the linear classifier fit to a model's latent space.

\subsection{Datasets}
\label{sec:datasets}
All of the datasets used in this paper, except for CalTech 101 which already contained binary images, contained gray-scale pixel feature values in the range of $[0,255]$. The images in these databases were first pre-processed by normalizing the pixels to the range of $[0,1]$ by dividing them by $255$ and finally converted to binary values by thresholding at $0.5$ as in \cite{bengio2013generalized}.

The MNIST dataset
contains $28\times28$ images containing handwritten digits across $10$ categories.
Fashion MNIST (FMNIST) \cite{xiao2017fashion}, which was proposed as a challenging drop-in replacement for MNIST, contains $28x28$ grey-scale images depicting clothing items (out of $10$ item classes).
Kuzushiji-MNIST (KMNIST) is another $28\times28$ image dataset containing  hand-drawn Japanese Kanji characters \cite{clanuwat2018deep}.
The Caltech 101 Silhouettes database
contains $16 \times 16$ binary pixel images across $100$ different categories.
Each training sample had $60000$ samples and the testing subset had $10000$ (the standard test split of each dataset was used in this paper). A validation subset of $2000$ samples was drawn from each training set to be used for tuning model meta-parameters.

\subsection{Baseline Model Descriptions}
\label{sec:baselines}
The baseline backprop-based models implemented for this article included a regularized auto-associative (autoencoding) network (RAE), a Gaussian variational autoencoder with fixed (spherical) variance (GVAE-CV) \cite{ghosh2020variational}, a Gaussian variational autoencoder (GVAE) \cite{kingma2013auto}, and an adversarial autoencoder (GAN-AE) \cite{makhzani2016adversarial}. All models were constrained to have four layers (see Appendix) like the GNCN. For all backprop-based models, the sizes of the layers in between the latent variable and the input layer $\widetilde{\mathbf{z}}^0$ were chosen such that the total synaptic weight count of the model was approximately equal to the GNCN, the linear rectifier was used for the internal activation function, i.e., $\phi^\ell(v) = max(0, v)$, and weights were initialized from a centered Gaussian distribution with a standard deviation that was tuned on held-out validation data. To further improve generalization ability, the decoder weights of all autoencoders were regularized and some autoencoder models had specific meta-parameters that were tuned using validation data. Notably, the GAN-AE was the only model that required a specialized gradient descent rule, i.e., Adam \cite{kingma2014adam}, in order to obtain good log likelihood and to stabilize training (stability is a known issue related to GANs \cite{goodfellow2014generative}).
Finally, we implemented an optimized Gaussian mixture model (GMM)
that was fit to the training data via expectation-maximization with the number of mixture components chosen based on preliminary experiments that yielded best performance. The only GMM implementation detail worthy of note was that we clipped the images sampled from the mixture model to lie in the range $[0,1]$ (improving likelihood slightly).

\subsection{Neural Generative Coding Model Description}
\label{sec:model}
In this section, we detail the sampling and image completion procedures for the GNCN. See the appendix for a detailed, formal specification of the GNCN as well as its learning and inference processes.

\noindent
\textbf{Sampling from the Model:} 
Since the optimization procedure, whether via gradient descent or another method is not guaranteed to find globally optimal parameter settings (since the objective function is not convex), the distribution of the latent state $\latent{z}^L$ of the final layer will not be Gaussian. To estimate it, after training on the data, we obtain the corresponding $\latent{z}_i^L$ value for each data point $\mathbf{x}_i$. We fit a Gaussian mixture model to this collection of values $\latent{z}_1^L,\dots, \latent{z}^L_N$ (where $N$ is the number of training points). Then, in order to reduce variance during sampling, we take the following approach. First, we sample $\latent{z}^L$ from the Gaussian mixture model, recursively set  $\mean{z}^{\ell-1} \gets g^{\ell-1}(\mathbf{W}^{\ell-1} \phi^{\ell}(\latent{z}^{\ell}))$ and output $\latent{z}^0$. This is similar to how variational auto-encoders are sampled in practice, where the input is a Gaussian and the output is the mean vector.

\noindent
\textbf{Image Completion:} 
In the event that incomplete sensory input $\mathbf{x}$ is provided to the GNCN, i.e., portions of $\mathbf{x}$ are masked out by the variable $\mathbf{m} \in \{0,1\}^{J_0 \times 1}$, we may infer the remaining portions of $\mathbf{x}$ by utilizing the output error neurons of the GNCN and treating the bottom sensory layer $\mathbf{z}^0$ as a partial latent state. Specifically, we update the missing portions, i.e., $1 - \mathbf{m}$, of $\mathbf{z}^0$ as follows:
\begin{align}
    \latent{z}^0 = \mathbf{x} \otimes \mathbf{m} + \Big( \latent{z}^0 + \beta( -\frac{\partial \psi }{\partial \mean{z}^0} ) \Big) \otimes (1 - \mathbf{m}) = \mathbf{x} \otimes \mathbf{m} + \Big( \latent{z}^0 - \beta \mathbf{e}^0 \Big) \otimes (1 - \mathbf{m}) \mbox{.}
\end{align}


\subsection{Programming}
The experiments, the baseline models, and the proposed computational framework were written in Python, using the Numpy and TensorFlow 2 libraries. Experiments were sped up using a GPU card.

\subsection{Data Availability}
The datasets used in this study are publicly available.

\bibliographystyle{acm}
