Conditional Finite Mixtures of Poisson Distributions for Context-Dependent Neural Correlations

TL;DR

This paper introduces a novel conditional finite mixture model of Poisson distributions that captures context-dependent neural correlations, improving understanding and prediction of neural population responses.

Contribution

The paper develops a new model that captures low-dimensional, context-dependent neural correlations and provides an efficient EM algorithm for fitting large neural datasets.

Findings

Successfully models stimulus-dependent correlations in macaque V1 neurons.

Incorporates nonlinear context-dependence for improved neural activity prediction.

Estimates correlation dimensionality via cross-validation.

Abstract

Parallel recordings of neural spike counts have revealed the existence of context-dependent noise correlations in neural populations. Theories of population coding have also shown that such correlations can impact the information encoded by neural populations about external stimuli. Although studies have shown that these correlations often have a low-dimensional structure, it has proven difficult to capture this structure in a model that is compatible with theories of rate coding in correlated populations. To address this difficulty we develop a novel model based on conditional finite mixtures of independent Poisson distributions. The model can be conditioned on context variables (e.g. stimuli or task variables), and the number of mixture components in the model can be cross-validated to estimate the dimensionality of the target correlations. We derive an expectation-maximization algorithm to efficiently fit the model to realistic amounts of data from large neural populations. We then demonstrate that the model successfully captures stimulus-dependent correlations in the responses of macaque V1 neurons to oriented gratings. Our model incorporates arbitrary nonlinear context-dependence, and can thus be applied to improve predictions of neural activity based on deep neural networks.

Equations

$${\theta }_{C,j}^{\ast }$$

$${\text{uptheta}}_{X,j}^{\ast }$$

$${\partial }_{{\text{uptheta}}_X}\log q_X(X_i)$$

$${\partial }_{{\text{uptheta}}_Y}\log q_X(X_i)$$

$${\partial }_{{\mathbf{\Theta }}_{XY}}\log q_X(X_i)$$

$$\begin{gathered} \mathcal{L}({\text{uptheta}}_X,{\text{uptheta}}_C,{\mathbf{\Theta }}_{XC},\text{uprho},{\text{upeta}}_{C,i},X_i,Z_i)={\mathbf{s}}_X(X_i)\cdot ({\text{uptheta}}_X+{\text{uptheta}}_{X\mid Z}(Z_i;\text{uprho})) \\ +{\text{uptheta}}_C\cdot {\text{upeta}}_{C,i}+{\mathbf{s}}_X(X_i)\cdot {\mathbf{\Theta }}_{XC}\cdot {\text{upeta}}_{C,i}-{\psi }_{XC}({\text{uptheta}}_X+{\text{uptheta}}_{X\mid Z}(Z_i;\text{uprho}),{\text{uptheta}}_C,{\mathbf{\Theta }}_{XC}) \end{gathered}$$

$$\underset{{\text{uptheta}}_N}{\mathrm{argmax}}\sum _{i=1}^{m_S}\mathcal{L}({\text{uptheta}}_N,{\text{uptheta}}_C,{\mathbf{\Theta }}_{NC},\text{uprho},{\text{upeta}}_{C,i},N_i,Z_i)={\text{uptheta}}_{N,0}^{†},$$

$$\underset{{\text{uptheta}}_{N,j}}{\mathrm{argmax}}\sum _{i=1}^{m_S}\mathcal{L}({\text{uptheta}}_N,{\text{uptheta}}_C,{\mathbf{\Theta }}_{NC},\text{uprho},{\text{upeta}}_{C,i},N_i,Z_i)={\text{uptheta}}_{N,j}^{†}-{\text{uptheta}}_{N,0}^{†},$$

$$\theta^{\dagger}_{N,j,k}=\log\Big{(}\frac{\sum_{i=1}^{m_{S}}\eta_{C,i,j}N_{i,k}}{\sum_{i=1}^{m_{S}}\eta_{C,i,j}e^{\theta_{N\mid Z,k}(Z_{i})}}\Big{)}.$$

$$\mathcal{L}_{i}^{\dagger}({\{{\bm{\uptheta}}^{\dagger}_{N,j}\}}_{j=0}^{m_{C}})=\sum_{j=0}^{m_{C}}\eta_{C,i,j}\big{(}N_{i}\cdot{\bm{\uptheta}}^{\dagger}_{N,j}-\psi_{N}({\bm{\uptheta}}^{\dagger}_{N,j}+{\bm{\uptheta}}_{N\mid Z}(Z_{i}))\big{)}\\ +N_{i}\cdot{\bm{\uptheta}}_{N\mid Z}(Z_{i})+{\bm{\uptheta}}^{\dagger}_{C}(Z_{i})\cdot{\bm{\upeta}}_{C,i}-\psi_{C}({\bm{\uptheta}}^{\dagger}_{C}(Z_{i})),$$

$$\underset{{\theta }_{N,j,k}^{†}}{\mathrm{argmax}}\sum _{i=1}^{m_S}{\mathcal{L}}_i^{†}({\{{\text{uptheta}}_{N,j}^{†}\}}_{j=0}^{m_C})$$

$$\displaystyle=\operatorname*{arg\,max}_{\theta^{\dagger}_{N,j,k}}\Big{\{}\theta^{\dagger}_{N,j,k}\frac{\sum_{i=1}^{m_{S}}\eta_{C,i,j}N_{i,k}}{\sum_{i=1}^{m_{S}}\eta_{C,i,j}e^{\theta_{N\mid Z,k}(Z_{i})}}-e^{\theta^{\dagger}_{N,j,k}}\Big{\}}.$$

Full text

1 Introduction

Computational neuroscientists have made significant progress in understanding how correlations amongst neural spike-counts influence information processing [1, 2]. In particular, correlated fluctuations in the responses to a fixed stimulus (termed noise correlations) depend on stimulus identity, and can both hinder and facilitate information processing in model neural circuits [3, 4, 5, 6, 7]. However, due to the challenges in assessing the strength and order of noise correlations empirically, it remains unclear how noise correlations affect computation in biological neural circuits.

Measuring correlations of every order in large neural populations is infeasible, and modelling neural correlations requires knowledge of the order and dimensionality of significant correlations. Most models measure exclusively pair-wise correlations [8, 9], which can often explain most of the variability in spike-count data [10, 11]. Moreover, dimensionality-reduction methods have shown that the complete set of pair-wise neural correlations often has a low-dimensional structure for both total [12] and noise [13, 14, 15, 16, 17] correlations. Nevertheless, higher-order correlations become more significant as both stimulus complexity and neural population size increase [8], motivating the use of higher-order models when modelling noise correlations in large populations of neurons.

Dimensionality reduction methods have provided many insights into neural correlations, but correlations must be studied through models of population codes in order to understand their effect on neural coding. Maximum entropy models of binary neural spiking drive much of the contemporary research on neural correlations in model population codes [8]. The flexibility of the maximum entropy framework allows the techniques for analyzing and fitting pair-wise models to be generalized to both sparse, higher-order models [18] and stimulus-dependent models [11]. Nevertheless, there is still a need for a rate-based model of population coding that can reliably capture the complex yet low-dimensional noise correlations found in neural data.

To address this we propose a conditional maximum entropy model, specifically a conditional finite mixture of independent Poisson distributions (CMP). When conditioned on a context variable (e.g. a stimulus) a CMP is a mixture of component collections of independent Poisson distributions [19]. A CMP with only one component reduces to a standard rate-based population code with Poisson neurons that are conditionally independent given the stimulus [20, 21], and adding components to the CMP introduces noise correlations with dimensionality controlled by the number of components.

By combining the properties of maximum entropy models with additional features of Poisson neurons we derive a hybrid expectation-maximization algorithm that can efficiently and reliably fit the CMP to data. We apply the CMP to neural population data recorded in macaque primary visual cortex (V1). We find that noise correlation structure is best captured with 3–5 components across several datasets, and that while overall dimensionality of correlations is largely stimulus-independent, the structure of noise correlations (i.e. the relative weights of the mixture components) depends on the stimulus.

2 Conditional Finite Mixtures of Independent Poisson Distributions

Many of the equations that describe how to analyze and train the CMP model can be solved by expressing the parameters of the model in an appropriate coordinate system. Two of the coordinate systems we consider are the mean and natural parameters that arise in the context of exponential families111Maximum entropy models are also known as exponential families, and we prefer that name in this paper.. We thus begin this section by reviewing exponential families, primarily for the purpose of developing notation (for a thorough development see [22, 23]). We continue by showing that the standard weighted-sum form of a finite mixture model is a third parameterization of a kind of exponential family known as an exponential family harmonium [24]. Finally, we introduce and develop conditional finite mixture models and CMPs, and we derive training algorithms for such models including the hybrid expectation-maximization algorithm for training CMPs.

2.1 Exponential Families

Consider a random variable $X\in\mathcal{X}$ with an unknown distribution $P_{X}$, and suppose we are given an independent and identically distributed sample ${\{X_{i}\}}_{i=1}^{m_{S}}$ from $P_{X}$ such that every $X_{i}\sim P_{X}$. We may model $P_{X}$ based on the sample ${\{X_{i}\}}_{i=1}^{m_{S}}$ by first defining a statistic $\mathbf{s}_{X}\colon\mathcal{X}\to{\text{H}}_{X}$, where ${\text{H}}_{X}$ is the codomain of $\mathbf{s}_{X}$ with dimension $m_{X}$. We then look for a probability distribution $Q_{X}$ that satisfies ${\mathbb{E}}_{Q}[\mathbf{s}_{X}(X)]=\frac{1}{m_{S}}\sum_{i=1}^{m_{S}}\mathbf{s}_{X}(X_{i})$, where ${\mathbb{E}}_{Q}[f(X)]=\int_{\mathcal{X}}fdQ_{X}$ is the expected value of $f(X)$ under $Q_{X}$. This is of course an under-constrained problem, but if we also assume that $Q_{X}$ must have maximum entropy, then we arrive at a family of distributions which uniquely satisfy the constraints for every value of $\frac{1}{m_{S}}\sum_{i=1}^{m_{S}}\mathbf{s}_{X}(X_{i})$.

An $m_{X}$-dimensional exponential family $\mathcal{M}_{X}$ is defined by the sufficient statistic $\mathbf{s}_{X}$, as well as a base measure $\mu_{X}$ which helps define integrals and expectations within the family. An exponential family is parameterized by a set of natural parameters $\Theta_{X}$, such that each element of the family $Q_{X}\in\mathcal{M}_{X}$ may be identified with some parameters ${\bm{\uptheta}}_{X}\in\Theta_{X}$. The density $q_{X}$ of the distribution $Q_{X}$ is given by $\log q_{X}(x)=\mathbf{s}_{X}(x)\cdot{\bm{\uptheta}}_{X}-\psi_{X}({\bm{\uptheta}}_{X})$, where $\psi_{X}({\bm{\uptheta}}_{X})=\log\int_{\mathcal{X}}e^{\mathbf{s}_{X}(x)\cdot{\bm{\uptheta}}_{X}}\mu_{X}(dx)$ is the log-partition function. Expectations of any $Q_{X}\in\mathcal{M}_{X}$ are then given by ${\mathbb{E}}_{Q}[f(X)]=\int_{\mathcal{X}}fdQ_{X}=\int_{\mathcal{X}}f\cdot q_{X}d\mu_{X}$. Because each $Q_{X}\in\mathcal{M}_{X}$ is uniquely defined by ${\mathbb{E}}_{Q}[\mathbf{s}_{X}(X)]$, the means of the sufficient statistic also parameterize $\mathcal{M}_{X}$. The space of all mean parameters is ${\text{H}}_{X}$, and we denote them by ${\bm{\upeta}}_{X}={\mathbb{E}}_{Q}[\mathbf{s}_{X}(X)]$. Finally, a sufficient statistic is minimal when its component functions are non-constant and linearly independent. If the sufficient statistic $\mathbf{s}_{X}$ of a given family $\mathcal{M}_{X}$ is minimal, then $\Theta_{X}$ and ${\text{H}}_{X}$ are isomorphic. Moreover, the transition functions between them $\bm{\uptau}_{X}\colon\Theta_{X}\to{\text{H}}_{X}$ and $\bm{\uptau}_{X}^{-1}\colon{\text{H}}_{X}\to\Theta_{X}$ are given by $\bm{\uptau}_{X}({\bm{\uptheta}}_{X})=\partial_{{\bm{\uptheta}}_{X}}\psi_{X}({\bm{\uptheta}}_{X})$, and $\bm{\uptau}_{X}^{-1}({\bm{\upeta}}_{X})=\partial_{{\bm{\upeta}}_{X}}\phi_{X}({\bm{\upeta}}_{X})$, where $\phi_{X}({\bm{\upeta}}_{X})={\mathbb{E}}_{Q}[\log q_{X}(X)]$ is the negative entropy of $Q_{X}$.

Before we conclude this subsection, we highlight the exponential families of categorical distributions $\mathcal{M}_{C}$ and of independent Poisson distributions $\mathcal{M}_{N}$. The $m_{C}$-dimensional categorical exponential family $\mathcal{M}_{C}$ contains all the distributions over integer values between 0 and $m_{C}$. The base measure of $\mathcal{M}_{C}$ is the counting measure, and the $j$th element of its sufficient statistic $\mathbf{s}_{C}(j)$ at $j$ is 1, and 0 for any other elements. The sufficient statistic $\mathbf{s}_{C}$ is thus a vector of all zeroes when $j=0$, and all zeroes except for the $j$th element when $j>0$. The $m_{N}$-dimensional family of independent Poisson distributions $\mathcal{M}_{N}$ contains distributions over $m_{N}$-dimensional vectors of natural numbers, where each element of the vector is independently Poisson distributed. The sufficient statistic of $\mathcal{M}_{N}$ is the identity function, and the base measure is given by $\mu_{N}(\{\mathbf{n}\})=\prod_{j=1}^{m_{N}}\frac{1}{n_{j}!}$.

2.2 Finite Mixture Models

In this subsection we show how finite mixture models can be expressed as latent-variable models such that the complete joint model is a particular kind of exponential family. To begin, let us consider an exponential family $\mathcal{M}_{X}$. The density of a finite mixture distribution $Q_{X}^{*}$ of $m_{C}+1$ distributions in $\mathcal{M}_{X}$ has the form $q_{X}^{*}(x)=\sum_{j=0}^{m_{C}}\eta^{*}_{C,j}q^{*}_{X,j}(x)$, where each distribution $Q^{*}_{X,j}\in\mathcal{M}_{X}$ is a mixture component with natural parameters ${\bm{\uptheta}}^{*}_{X,j}$, each $\eta^{*}_{C,j}$ is a mixture weight, and the weights satisfy the constraints $0<\eta^{*}_{C,j}<1$ and $\eta^{*}_{0}=1-\sum_{j=1}^{m_{C}}\eta^{*}_{C,j}$.

Building a model out of mixture distributions is theoretically problematic because swapping component indices has no effect on a mixture distribution, and their parameters are therefore not identifiable. Nevertheless, we may express mixture distributions as latent-variable distributions with the identifiable form $\sum_{j=0}^{m_{C}}\eta^{*}_{C,j}q^{*}_{X,j}(x)=\sum_{j=0}^{m_{C}}q^{*}_{C}(j)q^{*}_{X\mid C}(x\mid j)=\sum_{j=0}^{m_{C}}q^{*}_{XC}(x,j)$, where $\eta^{*}_{C,j}=q^{*}_{C}(j)$ and $q^{*}_{X,j}(x)=q^{*}_{X\mid C}(x\mid j)$. Moreover, $Q_{C}^{*}$ is a categorical distribution with mean parameters ${\bm{\upeta}}_{C}^{*}={(\eta^{*}_{C,j})}_{j=1}^{m_{C}}$, and $Q_{C}^{*}$ may thus be described by the natural parameters ${\bm{\uptheta}}_{C}^{*}=\bm{\uptau}_{C}^{-1}({\bm{\upeta}}_{C}^{*})$. Given $\mathcal{M}_{X}$ and $m_{C}$, we therefore define a finite mixture model $\mathcal{M}_{XC}^{*}$ as the set of all distributions $Q^{*}_{XC}$, and parameterize it by the mixture weights ${\bm{\uptheta}}_{C}^{*}$ and mixture components parameters ${\{{\bm{\uptheta}}^{*}_{X,j}\}}_{j=0}^{m_{C}}$.

Now, an exponential family harmonium is a kind of product exponential family which includes restricted Boltzmann machines as a special case [24]. We may construct an exponential family harmonium $\mathcal{M}_{XY}$ out of $\mathcal{M}_{X}$ and $\mathcal{M}_{Y}$ by defining the base measure of $\mathcal{M}_{XY}$ as the product measure $\mu_{X}\times\mu_{Y}$, and by defining the sufficient statistic of $\mathcal{M}_{XY}$ as the vector which contains the concatenation of all the elements in $\mathbf{s}_{X}$, $\mathbf{s}_{Y}$, and the outer product $\mathbf{s}_{X}\otimes\mathbf{s}_{Y}$. Intuitively, $\mathcal{M}_{XY}$ contains all the distributions with densities of the form $q_{XY}(x,y)\propto e^{\mathbf{s}_{X}(x)\cdot{\bm{\uptheta}}_{X}+\mathbf{s}_{Y}(y)\cdot{\bm{\uptheta}}_{Y}+\mathbf{s}_{X}(x)\cdot{\bm{\Theta}}_{XY}\cdot\mathbf{s}_{Y}(y)}$, where ${\bm{\uptheta}}_{X}$, ${\bm{\uptheta}}_{Y}$, and ${\bm{\Theta}}_{XY}$ are the natural parameters of $Q_{XY}\in\mathcal{M}_{XY}$. As we show, a harmonium defined partially by the categorical exponential family is equivalent to a mixture model.

Theorem 1**.**

Let $\mathcal{M}_{X}$ be an exponential family, let $\mathcal{M}_{C}$ be the categorical exponential family of dimension $m_{C}$, let $\mathcal{M}_{XC}^{*}$ be the mixture model defined by $\mathcal{M}_{X}$ and $m_{C}$, and let $\mathcal{M}_{XC}$ be the exponential family harmonium defined by $\mathcal{M}_{X}$ and $\mathcal{M}_{C}$. Then $Q_{XC}\in\mathcal{M}_{XC}$ with natural parameters ${\bm{\uptheta}}_{X}$, ${\bm{\uptheta}}_{C}$, ${\bm{\Theta}}_{XC}$ and $Q_{XC}\in\mathcal{M}_{XC}^{*}$ with mixture parameters ${\bm{\uptheta}}^{*}_{C}$ and ${\{{\bm{\uptheta}}^{*}_{X,j}\}}_{j=0}^{m_{C}}$ if and only if

[TABLE]

where $\theta_{C,j}$ and $\theta^{*}_{C,j}$ are the $j$th elements of ${\bm{\uptheta}}_{C}$ and ${\bm{\uptheta}}_{C}^{*}$, ${\bm{\uptheta}}_{X,j}$ is the $j$th row of ${\bm{\Theta}}_{XC}$, and ${\bm{\uptheta}}_{X,0}=\mathbf{0}$.

Proof.

We identify $Q_{XC}$ and $Q^{*}_{XC}$ by identifying the parameters of their conditionals and marginals. First note that the conditional distribution $Q_{X\mid C=j}$ at $j$ is an element of $\mathcal{M}_{X}$ with parameters ${\bm{\uptheta}}_{X}+{\bm{\Theta}}_{XC}\cdot\mathbf{s}_{C}(j)$ [24], and that the sufficient statistic $\mathbf{s}_{C}$ in this expression essentially selects a row from ${\bm{\Theta}}_{XC}$, or returns $\mathbf{0}$. Thus, where ${\bm{\uptheta}}_{X,j}$ is the $j$th row of ${\bm{\Theta}}_{XC}$, $Q_{X\mid C=j}=Q^{*}_{X\mid C=j}$ if and only if ${\bm{\uptheta}}^{*}_{X,0}={\bm{\uptheta}}_{X}$, and ${\bm{\uptheta}}^{*}_{X,j}={\bm{\uptheta}}_{X,j}+{\bm{\uptheta}}_{X}$ for $j>0$.

To equate $Q^{*}_{C}$ and $Q_{C}$, first note that the marginal density of a harmonium distribution satisfies $\log q_{C}(j)={\bm{\uptheta}}_{C}\cdot\mathbf{s}_{C}(j)+\psi_{X}({\bm{\uptheta}}_{X}+{\bm{\Theta}}_{XC}\cdot\mathbf{s}_{C}(j))-\psi_{XC}({\bm{\uptheta}}_{X},{\bm{\uptheta}}_{C},{\bm{\Theta}}_{XC})$ [24]. Because $\mathcal{M}_{XC}$ is partially defined by $\mathcal{M}_{C}$, we may show with a bit of algebra that $\psi_{XC}({\bm{\uptheta}}_{X},{\bm{\uptheta}}_{C},{\bm{\Theta}}_{XC})=\psi_{X}({\bm{\uptheta}}_{X})+\psi_{C}({\bm{\uptheta}}_{C}^{*})$, where ${\bm{\uptheta}}^{*}_{C}$ satisfies Equation 1, and by substituting this into the expression for the harmonium marginal, we may conclude that $Q_{C}=Q^{*}_{C}$ if and only if Equation 1 is satisfied. ∎

The advantage of formulating finite mixture models as exponential family harmoniums is that we may apply the theory of harmoniums and exponential families to analyzing and training them. For example, given a sample ${\{X_{i}\}}_{i=1}^{m_{S}}$ and a harmonium distribution $Q_{XY}\in\mathcal{M}_{XY}$ with parameters ${\bm{\uptheta}}_{X}$, ${\bm{\uptheta}}_{Y}$, and ${\bm{\Theta}}_{XY}$, an iteration of the expectation-maximization algorithm (EM) may be formulated as:

Expectation Step:

compute the unobserved means ${\bm{\upeta}}_{Y,i}=\bm{\uptau}_{Y}({\bm{\uptheta}}_{Y}+\mathbf{s}_{X}(X_{i})\cdot{\bm{\Theta}}_{XY})$ for every $i$,

Maximization Step:

evaluate $\bm{\uptau}^{-1}_{XY}(\frac{1}{m_{S}}\sum_{i=1}^{m_{S}}\mathbf{s}_{X}(X_{i}),\frac{1}{m_{S}}\sum_{i=1}^{m_{S}}{\bm{\upeta}}_{Y,i},\frac{1}{m_{S}}\sum_{i=1}^{m_{S}}\mathbf{s}_{X}(X_{i})\otimes{\bm{\upeta}}_{Y,i})$.

On the other hand, the stochastic log-likelihood gradients of the parameters of $Q_{XY}$ are

[TABLE]

where $({\bm{\upeta}}_{X},{\bm{\upeta}}_{Y},\mathbf{H}_{XY})=\bm{\uptau}_{XY}({\bm{\uptheta}}_{X},{\bm{\uptheta}}_{Y},{\bm{\Theta}}_{XY})$.

Whether or not the various transition functions and their inverses can be evaluated depends on the defining manifolds $\mathcal{M}_{X}$ and $\mathcal{M}_{Y}$. When $\mathcal{M}_{Y}=\mathcal{M}_{C}$ is the categorical family, its transition function $\bm{\uptau}_{C}$ is computable, and whether $\bm{\uptau}_{XY}$ or $\bm{\uptau}_{XY}^{-1}$ are computable is reducible to whether $\bm{\uptau}_{X}$ or $\bm{\uptau}_{Y}^{-1}$ are computable, respectively. EM is typically preferred when maximizing the likelihood of finite mixture model parameters, especially for finite mixtures of normal distributions. However, gradient descent can perform just as well as EM for certain mixture models, as we demonstrate in figure 1.

2.3 Conditional Finite Mixture Models

According to Equations 1 and 2 we may induce context dependence in all the mixture parameters of a mixture distribution $Q_{XC}\in\mathcal{M}_{XC}$ by allowing ${\bm{\uptheta}}_{X}$ to depend on additional variables. We thus define a conditional finite mixture model $\mathcal{M}_{XC\mid Z}$ as the set of all the conditional distributions $Q_{XC\mid Z}$ with densities $q_{XC\mid Z}(x,j\mid z)\propto e^{\mathbf{s}_{X}(x)\cdot({\bm{\uptheta}}_{X}+{\bm{\uptheta}}_{X\mid Z}(z))+{\bm{\uptheta}}_{C}\cdot\mathbf{s}_{C}(j)+\mathbf{s}_{X}(x)\cdot{\bm{\Theta}}_{XC}\cdot\mathbf{s}_{C}(j)}$, where ${\bm{\uptheta}}_{X\mid Z}\colon\mathcal{Z}\to\Theta_{X}$ is defined by an additional set of parameters ${\bm{\uprho}}$. In the conditional setting we aim to maximize the conditional log-likelihood $\sum_{i=1}^{m_{S}}q_{X\mid Z}(X_{i}\mid Z_{i})$ of the parameters ${\bm{\uptheta}}_{X}$, ${\bm{\uptheta}}_{C}$, ${\bm{\Theta}}_{XC}$, and ${\bm{\uprho}}$, given a sample ${\{(X_{i},Z_{i})\}}_{i=1}^{m_{S}}$ from a target conditional distribution $P_{X\mid Z}$. Stochastic gradient ascent on the conditional log-likelihood remains relatively straightforward; given a sample point $(X_{i},Z_{i})$, one must simply backpropagate the error-gradient $S_{X}(X_{i})-{\bm{\upeta}}_{X}(Z_{i})$ through the function ${\bm{\uptheta}}_{X\mid Z}$, where $({\bm{\upeta}}_{X}(Z_{i}),{\bm{\upeta}}_{C}(Z_{i}),\mathbf{H}_{XC}(Z_{i}))=\bm{\uptau}_{XC}({\bm{\uptheta}}_{X}+{\bm{\uptheta}}_{X\mid Z}(Z_{i}),{\bm{\uptheta}}_{C},{\bm{\Theta}}_{XC})$.

The expectation step of the EM algorithm also remains trivial for conditional mixtures. It is easy to show that for any $Q_{XC\mid Z}\in\mathcal{M}_{XC\mid Z}$ the conditional distribution $Q_{C\mid X,Z}=Q_{C\mid X}$, and we may therefore continue to evaluate the unobserved means as ${\bm{\upeta}}_{C,i}=\bm{\uptau}_{C}({\bm{\uptheta}}_{C}+\mathbf{s}_{X}(X_{i})\cdot{\bm{\Theta}}_{XC})$. On the other hand, the maximization step can be expressed as maximizing the sum over $i$ of the function

[TABLE]

with respect to the parameters ${\bm{\uptheta}}_{X}$, ${\bm{\uptheta}}_{C}$, ${\bm{\Theta}}_{XC}$, and ${\bm{\uprho}}$. It is more difficult to evaluate the maximization step for conditional finite mixtures due to the dependence on $Z_{i}$, although we may still find maxima of the objective by gradient ascent. Nevertheless, when $\mathcal{M}_{X}=\mathcal{M}_{N}$ is the family of $m_{N}$ independent Poisson distributions, we can fact evaluate the maximization step in closed-form with respect to the parameters ${\bm{\uptheta}}_{N}$ and ${\bm{\Theta}}_{NC}$. We refer to the conditional finite mixture model $\mathcal{M}_{NC\mid Z}$ as the family of conditional finite mixtures of independent Poisson distributions (CMPs).

Theorem 2**.**

Let $\mathcal{M}_{NC\mid Z}$ be a CMP model defined by $m_{C}$, $m_{N}$, and $\mathcal{Z}$. Then

[TABLE]

and

[TABLE]

where

[TABLE]

Proof.

Similar to Theorem 1, consider the change of variables ${\bm{\uptheta}}^{\dagger}_{N,0}={\bm{\uptheta}}_{N}$,${\bm{\uptheta}}^{\dagger}_{N,j}={\bm{\uptheta}}_{N,j}+{\bm{\uptheta}}_{N}$, and $\theta^{\dagger}_{C,j}(Z_{i})=\theta_{C,j}+\psi_{N}({\bm{\uptheta}}_{N,j}+{\bm{\uptheta}}_{N}+{\bm{\uptheta}}_{N\mid Z}(Z_{i}))-\psi_{N}({\bm{\uptheta}}_{N}+{\bm{\uptheta}}_{N\mid Z}(Z_{i}))$. We may express our optimization in this alternative parameterization as

[TABLE]

such that $\max_{{\bm{\uptheta}}_{N},{\bm{\Theta}}_{NC}}\sum_{i=1}^{m_{S}}\mathcal{L}=\max_{{\{{\bm{\uptheta}}^{\dagger}_{N,j}\}}_{j=0}^{m_{C}}}\sum_{i=1}^{m_{S}}\mathcal{L}_{i}^{\dagger}$. Moreover, the log-partition function of $\mathcal{M}_{N}$ is given by $\psi_{N}({\bm{\uptheta}}_{N})=\sum_{k=1}^{m_{N}}e^{\theta_{N,k}}$ and we may therefore conclude that

[TABLE]

This is an $\operatorname*{arg\,max}$ variant of a Legendre transform and its solution is \log\big{(}\frac{\sum_{i=1}^{m_{S}}\eta_{C,i,j}N_{i,k}}{\sum_{i=1}^{m_{S}}\eta_{C,i,j}e^{\theta_{N\mid Z,k}(Z_{i})}}\big{)}. ∎

2.4 Training Algorithms

Based on these results, we propose three algorithms for training conditional finite mixture models, the third of which is specific to CMPs, and we abbreviate them as EM, SGD, and Hybrid. The first approximates expectation-maximization (EM) by computing the expectation step in closed-form, and approximating the maximization step by gradient ascent of the objective in Definition 6. The second is stochastic gradient descent (SGD) of the negative log-likelihood based on Equations 3, 4, and 5. The third is the Hybrid algorithm based on Theorem 2. The Hybrid algorithm alternates between one of two phases every training epoch: the first phase is simply SGD with respect to all the CMP parameters, and the second phase is to update the parameters ${\bm{\uptheta}}_{N}$ and ${\bm{\Theta}}_{NC}$ based on Theorem 2.

3 Applications to Synthetic and Real Data

In this section we model synthetic and real data with a CMP. In both cases we define the context variable as an element of the half-circle such that $\mathcal{Z}=[0,\pi]$, which represents e.g. the orientation of a grating, as is common in V1 experiments. We define the natural parameter function by ${\bm{\uptheta}}_{N\mid Z}(z)={\bm{\Theta}}_{NZ}\cdot\mathbf{s}_{Z}(z)$, where $\mathbf{s}_{Z}(z)=(\cos(2z),\sin(2z))$. This choice of ${\bm{\uptheta}}_{N\mid Z}$ allows us to express the conditional firing rates of the model neurons as ${\mathbb{E}}_{P}[N_{k}\mid Z=z,C=j]=\gamma_{j,k}f_{k}(z)$, where $\gamma_{j,k}$ is the gain and $f_{k}(z)\propto e^{\rho_{k}\cos(x-\mu_{k})}$ is a von Mises density function with preferred stimulus $\mu_{k}$ and precision $\rho_{k}$. The dependence of the model on the categorical variable is expressed solely through gain components $\bm{\upgamma}_{j}={(\gamma_{j,k})}_{k=1}^{m_{N}}$, and where ${\mathbb{E}}_{P}[N_{k}\mid Z=z]$ is the tuning curve of the $i$th neuron, this form ensures that each tuning curve retains a von Mises shape regardless of correlation structure.

All simulations were run on a Lenovo P51 laptop with an Intel Xeon processor, and individual simulations of the algorithms in Subsection 2.4 took on the order of tens of seconds to execute. All algorithms were implemented in the Haskell programming language.

3.1 Synthetic Data

We first consider synthetic data generated from a CMP $P_{CN\mid Z}$ with $m_{N}=20$ neurons and $m_{C}=7$ mixture components. The preferred stimuli $\mu_{k}$ of the tuning curves are drawn from a uniform distribution over $\mathcal{Z}$, and the precisions $\rho_{k}$ and gains $\gamma_{j,k}$ are log-normally distributed with a mean of 0.8 and 2, respectively. We plot the von Mises part $f_{k}$ of each CMP tuning curve in Figure 2A. In Figure 2B we plot $p_{C\mid Z}(j\mid z)$ for every $j$, as a function of $z$. This provides a low-dimensional picture of the stimulus-dependence of the correlations; the number of active components determines the dimensionality of the correlations, and their identity determines the structure of correlations.

To synthesize a dataset consistent with typical recordings in visual cortex we synthesize 62 responses to 8 stimuli tiled evenly over the half-circle, for a total of $m_{S}=496$ sample points ${\{(N_{i},Z_{i})\}}_{i=1}^{m_{S}}$ from $P_{N\mid Z}$. In Figure 2C we plot the descent of the resulting negative log-likelihood of the CMP parameters for each of the three algorithms proposed in Subsection 2.4. Each epoch of each algorithm involves exactly one pass over the dataset, ensuring that their computation time is of a similar order. In each epoch the dataset is broken up into randomized minibatches of 50 samples, resulting in ten gradient steps per epoch for each algorithm, except for the closed-form phase of the Hybrid algorithm which processes the data as a single batch. Gradient descents are implemented with the Adam algorithm [25] with standard momentum parameters and a learning rate of 0.005, and the Adam parameters are reset every epoch. To avoid local minima we run 10 descents in parallel for each algorithm and select the best one.

In addition, we plot the true negative log-likelihood $-\sum_{i=1}^{m_{S}}\log p_{N\mid Z}(N_{i}\mid Z_{i})$ in Figure 2C to establish a training lower-bound. Since a CMP $Q^{0}_{N\mid Z}$ with $m_{C}+1=1$ component reduces to a set of conditionally independent Poisson neurons, and fitting such a CMP is a convex optimization problem, we also plot an upper-bound in the form of $-\sum_{i=1}^{m_{S}}\log q_{N\mid Z}^{0}(N_{i}\mid Z_{i})$ for the optimal $Q^{0}_{N\mid Z}$. This upper-bound tells us how much the addition of correlations to the model improves the quality of the fit. As seen in Figure 2, the Hybrid algorithm converges more quickly and to a lower value than either SGD or EM. All algorithms overfit the data relative to the lower-bound, which is unsurprising given the small sample size. We later apply cross-validation in our analysis of real data to help avoid this.

In Figure 2D we plot $q_{C\mid Z}(j\mid z)$ for the CMP model trained with the Hybrid algorithm. We highlight two curves to emphasize that the weight-curves learned by model are similar to those of the ground truth CMP. In Figure 2E we plot stimulus-dependent, paired correlation matrices of the ground truth and hybrid-trained models. As can be seen, the learned-model and ground truth correlations are nearly identical. Also note that when one mixture weight dominates, the correlations disappear. Finally, in Figure 3 we repeat the previous analysis on a CMP target and model with $m_{N}=200$ neurons. As can be seen, the results are essentially identical to those in figure 2. As such, at least for synthetic data, a small number of sample points (e.g. 496) are sufficient for modelling low-dimensional noise correlations in large populations of neurons.

3.2 Response Recordings from Macaque Primary Visual Cortex

In this subsection we repeat our analysis from the previous subsection on multi-electrode recordings of macaque primary visual cortex (V1) in response to the oriented grating stimuli (the data were originally presented in [26]). We analyze 8 datasets corresponding to 8 multi-electrode recording sessions of between 28–76 well-tuned neurons. Each dataset is composed of 80 repeated trials of 16 stimuli222Each set of 80 trials pools 4 distinct phases of the stimulus, which could inflate measured correlations. Nevertheless, recorded neurons were found to be roughly phase independent [26]., where the stimuli are spread evenly over the half circle.

In Figure 4A we depict the 10-fold cross-validation of the log-likelihood of each of the 8 datasets as a function of the number of components $m_{C}+1$ in the Hybrid-trained CMP model, and we subtract the 1-component value from the value computed for the remaining components. This figure thus quantifies the amount of information (in nats) that is gained about the true data distribution by modelling the data with the given number of components. As can be seen, between 3–5 components is optimal for the various datasets before information gain yields to overfitting. In Figure 4B we plot the results of fitting the complete selected dataset with the three proposed algorithms. The Hybrid algorithm continues to perform well and converge quickly, while it appears as though EM and SGD require significantly more computation to find good solutions.

In Figure 4C we plot $q_{C\mid Z}(j\mid z)$ for the hybrid-trained model $Q_{CN\mid Z}$, but in this case we colour the four curves to more easily distinguish them. The curves indicate that the correlations of the true data distribution are stimulus-dependent, but less so than those of our random models. Moreover, at no stimulus value does a single component dominate, and the effective dimensionality is largely stimulus-independent. Figure 4D shows that the CMP finds a subtle yet significant amount of stimulus-dependence in model correlations, which are not clearly visible in the empirical correlations.

4 Conclusion

In this paper we introduced a novel model (CMPs) of context-dependent neural correlations, derived an expectation-maximization method for training it, and demonstrated that CMPs can capture significant correlation structure amongst neurons. We also demonstrated that EM and SGD can often perform equally well in the context of finite mixture models, which stands in contrast with common practice [27]. In future work we will compare the correlations learned with a CMP to the results of alternative dimensionality-reduction techniques for context-dependent neural correlations [11, 13, 12]. At the same time, CMPs are uniquely suited for studying rate-based neural coding. For example, information-limiting correlations [7] can be directly incorporated into a CMP by defining the gain components as identical but shifted with respect to neuron index.

In our applications in Section 3, the stimulus-dependent parameters had a simple generalized linear form. Theoretically, however, these parameters could be modelled by a deep neural network, thereby allowing the output of the deep network to exhibit correlations. When the output of the network has a complex structure, such as when modelling the stimulus-dependent neural responses of mid-level brain regions to naturalistic images [28], incorporating correlations could significantly improve model predictions. Moreover, because CMPs model noise correlations in a manner that is compatible with theories of population coding, such a combined model can extend our understanding of neural computation beyond low-level sensory areas and low-dimensional variables.