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

## Sacha Sokoloski

## Ruben Coen-Cagli

Systems and Computational Biology Albert Einstein College of Medicine The Bronx, NY 10461

Systems and Computational Biology Albert Einstein College of Medicine The Bronx, NY 10461

sacha.sokoloski@einstein.yu.edu

ruben.coen-cagli@einstein.yu.edu

## Abstract

Parallel recordings of neural spike counts have revealed the existence of contextdependent 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 stimulusdependent 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.

## 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-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-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

Preprint. Under review.

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 families . We thus begin this section by reviewing exponential families, primarily for the purpose 1 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 ∈ X with an unknown distribution P X , and suppose we are given an independent and identically distributed sample { X i } m S i =1 from P X such that every X i ∼ P X . We may model P X based on the sample { X i } m S i =1 by first defining a statistic s X : X → H X , where H X is the codomain of s X with dimension m X . We then look for a probability distribution Q X that satisfies E Q [ s X ( X )] = 1 m S ∑ m S i =1 s X ( X i ) , where E Q [ f ( X )] = ∫ 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 1 m S ∑ m S i =1 s X ( X i ) .

An m X -dimensional exponential family M X is defined by the sufficient statistic s X , as well as a base measure µ X which helps define integrals and expectations within the family. An exponential family is parameterized by a set of natural parameters Θ X , such that each element of the family Q X ∈ M X may be identified with some parameters θ X ∈ Θ X . The density q X of the distribution Q X is given by log q X ( x ) = s X ( x ) · θ X -ψ X ( θ X ) , where ψ X ( θ X ) = log ∫ X e s X ( x ) · θ X µ X ( dx ) is the log-partition function. Expectations of any Q X ∈ M X are then given by E Q [ f ( X )] = ∫ X fdQ X = ∫ X f · q X dµ X . Because each Q X ∈ M X is uniquely defined by E Q [ s X ( X )] , the means of the sufficient statistic also parameterize M X . The space of all mean parameters is H X , and we denote them by η X = E Q [ s X ( X )] . Finally, a sufficient statistic is minimal when its component functions are non-constant and linearly independent. If the sufficient statistic s X of a given family M X is minimal, then Θ X and H X are isomorphic. Moreover, the transition functions between them τ X : Θ X → H X and τ -1 X : H X → Θ X are given by τ X ( θ X ) = ∂ θ X ψ X ( θ X ) , and τ -1 X ( η X ) = ∂ η X φ X ( η X ) , where φ X ( η X ) = 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 M C and of independent Poisson distributions M N . The m C -dimensional categorical exponential

1 Maximum entropy models are also known as exponential families, and we prefer that name in this paper.

family M C contains all the distributions over integer values between 0 and m C . The base measure of M C is the counting measure, and the j th element of its sufficient statistic s C ( j ) at j is 1, and 0 for any other elements. The sufficient statistic s C is thus a vector of all zeroes when j = 0 , and all zeroes except for the j th element when j &gt; 0 . The m N -dimensional family of independent Poisson distributions 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 M N is the identity function, and the base measure is given by µ N ( { n } ) = ∏ m N j =1 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 M X . The density of a finite mixture distribution Q ∗ X of m C +1 distributions in M X has the form q ∗ X ( x ) = ∑ m C j =0 η ∗ C,j q ∗ X,j ( x ) , where each distribution Q ∗ X,j ∈ M X is a mixture component with natural parameters θ ∗ X,j , each η ∗ C,j is a mixture weight , and the weights satisfy the constraints 0 &lt; η ∗ C,j &lt; 1 and η ∗ 0 = 1 -∑ m C j =1 η ∗ 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 ∑ m C j =0 η ∗ C,j q ∗ X,j ( x ) = ∑ m C j =0 q ∗ C ( j ) q ∗ X C | ( x | j ) = ∑ m C j =0 q ∗ XC ( x, j ) , where η ∗ C,j = q ∗ C ( j ) and q ∗ X,j ( x ) = q ∗ X C | ( x | j ) . Moreover, Q ∗ C is a categorical distribution with mean parameters η ∗ C = ( η ∗ C,j ) m C j =1 , and Q ∗ C may thus be described by the natural parameters θ ∗ C = τ -1 C ( η ∗ C ) . Given M X and m C , we therefore define a finite mixture model M ∗ XC as the set of all distributions Q ∗ XC , and parameterize it by the mixture weights θ ∗ C and mixture components parameters { θ ∗ X,j } m C j =0 .

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 M XY out of M X and M Y by defining the base measure of M XY as the product measure µ X × µ Y , and by defining the sufficient statistic of M XY as the vector which contains the concatenation of all the elements in s X , s Y , and the outer product s X ⊗ s Y . Intuitively, M XY contains all the distributions with densities of the form q XY ( x, y ) ∝ e s X ( x ) · θ X + s Y ( y ) · θ Y + s X ( x ) · Θ XY · s Y ( y ) , where θ X , θ Y , and Θ XY are the natural parameters of Q XY ∈ M XY . As we show, a harmonium defined partially by the categorical exponential family is equivalent to a mixture model.

Theorem 1. Let M X be an exponential family, let M C be the categorical exponential family of dimension m C , let M ∗ XC be the mixture model defined by M X and m C , and let M XC be the exponential family harmonium defined by M X and M C . Then Q XC ∈ M XC with natural parameters θ X , θ C , Θ XC and Q XC ∈ M ∗ XC with mixture parameters θ ∗ C and { θ ∗ X,j } m C j =0 if and only if

$$\theta _ { C, j } ^ { * } = \theta _ { C, j } + \psi _ { X } ( \theta _ { X, j } + \theta _ { X } ) - \psi _ { X } ( \theta _ { X } ),$$

$$\theta _ { X, j } ^ { * } & = \theta _ { X, j } + \theta _ { X },$$

where θ C,j and θ ∗ C,j are the j th elements of θ C and θ ∗ C , θ X,j is the j th row of Θ XC , and θ X, 0 = 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 C | = j at j is an element of M X with parameters θ X + Θ XC · s C ( j ) [24], and that the sufficient statistic s C in this expression essentially selects a row from Θ XC , or returns 0 . Thus, where θ X,j is the j th row of Θ XC , Q X C | = j = Q ∗ X C | = j if and only if θ ∗ X, 0 = θ X , and θ ∗ X,j = θ X,j + θ X for j &gt; 0 .

To equate Q ∗ C and Q C , first note that the marginal density of a harmonium distribution satisfies log q C ( j ) = θ C · s C ( j ) + ψ X ( θ X + Θ XC · s C ( j )) -ψ XC ( θ X , θ C , Θ XC ) [24]. Because M XC is partially defined by M C , we may show with a bit of algebra that ψ XC ( θ X , θ C , Θ XC ) = ψ X ( θ X ) + ψ C ( θ ∗ C ) , where θ ∗ 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

Figure 1: Training a finite mixture of an independent product of 2 von Mises distributions on synthetic data with EM vs batch gradient descent (GD). Ellipsoids indicate the precisions of the component densities of the true mixture (black), the EM-trained mixture (red), and the GD-trained mixture (dashed blue). Black circles are synthetic data. EM and GD find the same solutions starting from the same initial conditions and with similar computation time.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00637v2-figure-b7c76713-e9ec-4988-aded-e39c2ec52a7a.jpg)

example, given a sample { X i } m S i =1 and a harmonium distribution Q XY ∈ M XY with parameters θ X , θ Y , and Θ XY , an iteration of the expectation-maximization algorithm (EM) may be formulated as:

Expectation Step: compute the unobserved means η Y,i = τ Y ( θ Y + s X ( X i ) · Θ XY ) for every i ,

$$\text{Maximization Step: evaluate } \tau _ { X Y } ^ { - 1 } ( \frac { 1 } { m _ { S } } \sum _ { i = 1 } ^ { m _ { S } } s _ { X } ( X _ { i } ), \frac { 1 } { m _ { S } } \sum _ { i = 1 } ^ { m _ { S } } \eta _ { Y, i }, \frac { 1 } { m _ { S } } \sum _ { i = 1 } ^ { m _ { S } } s _ { X } ( X _ { i } ) \otimes \eta _ { Y, i } ).$$

On the other hand, the stochastic log-likelihood gradients of the parameters of Q XY are

$$\partial _ { \theta _ { X } } \log q _ { X } ( X _ { i } ) = \mathbf s _ { X } ( X _ { i } ) - \mathbf \eta _ { X },$$

$$\partial _ { \theta _ { Y } } \log q _ { X } ( X _ { i } ) = \tau _ { Y } ( \theta _ { Y } + s _ { X } ( X _ { i } ) \cdot \Theta _ { X Y } ) - \eta _ { Y },$$

$$\partial _ { \Theta _ { X Y } } \log q _ { X } ( X _ { i } ) = s _ { X } ( X _ { i } ) \otimes \tau _ { Y } ( \theta _ { Y } + s _ { X } ( X _ { i } ) \cdot \Theta _ { X Y } ) - H _ { X Y }.$$

where ( η X , η Y , H XY ) = τ XY ( θ X , θ Y , Θ XY ) .

Whether or not the various transition functions and their inverses can be evaluated depends on the defining manifolds M X and M Y . When M Y = M C is the categorical family, its transition function τ C is computable, and whether τ XY or τ -1 XY are computable is reducible to whether τ X or τ -1 Y 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 ∈ M XC by allowing θ X to depend on additional variables. We thus define a conditional finite mixture model M XC Z | as the set of all the conditional distributions Q XC Z | with densities q XC Z | ( x, j | z ) ∝ e s X ( x ) ( · θ X + θ X Z | ( z ))+ θ C · s C ( j )+ s X ( x ) · Θ XC · s C ( j ) , where θ X Z | : Z → Θ X is defined by an additional set of parameters ρ . In the conditional setting we aim to maximize the conditional log-likelihood ∑ m S i =1 q X Z | ( X i | Z i ) of the parameters θ X , θ C , Θ XC , and ρ , given a sample { ( X , Z i i ) } m S i =1 from a target conditional distribution P X Z | . Stochastic gradient ascent on the conditional log-likelihood remains relatively straightforward; given a sample point ( X , Z i i ) , one must simply backpropagate the error-gradient S X ( X i ) -η X ( Z i ) through the function θ X Z | , where ( η X ( Z i ) , η C ( Z i ) , H XC ( Z i )) = τ XC ( θ X + θ X Z | ( Z i ) , θ C , Θ XC ) .

The expectation step of the EM algorithm also remains trivial for conditional mixtures. It is easy to show that for any Q XC Z | ∈ M XC Z | the conditional distribution Q C X,Z | = Q C X | , and we may therefore continue to evaluate the unobserved means as η C,i = τ C ( θ C + s X ( X i ) · Θ XC ) . On the other hand, the maximization step can be expressed as maximizing the sum over i of the function

$$\mathcal { L } ( \theta _ { X }, \theta _ { C }, \Theta _ { X C }, \rho, \eta _ { C, i }, X _ { i }, Z _ { i } ) = s _ { X } ( X _ { i } ) \cdot ( \theta _ { X } + \theta _ { X | Z } ( Z _ { i } ; \rho ) ) \\ + \, \theta _ { C } \cdot \eta _ { C, i } + s _ { X } ( X _ { i } ) \cdot \Theta _ { X C } \cdot \eta _ { C, i } - \psi _ { X C } ( \theta _ { X } + \theta _ { X | Z } ( Z _ { i } ; \rho ), \theta _ { C }, \Theta _ { X C } ) \quad ( 6 )$$

with respect to the parameters θ X , θ C , Θ XC , and ρ . 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 M X = 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 θ N and Θ NC . We refer to the conditional finite mixture model M NC Z | as the family of conditional finite mixtures of independent Poisson distributions (CMPs).

Theorem 2. Let M NC Z | be a CMP model defined by m C , m N , and Z . Then

$$\underset { \theta _ { N } } { \arg \max } \sum _ { i = 1 } ^ { m _ { S } } \mathcal { L } ( \theta _ { N }, \theta _ { C }, \Theta _ { N C }, \rho, \eta _ { C, i }, N _ { i }, Z _ { i } ) = \theta _ { N, 0 } ^ { \dagger },$$

and

$$\underset { \theta _ { N, j } } { \arg \max } \underset { i = 1 } { \overset { m _ { S } } { \sum } } \mathcal { L } ( \theta _ { N }, \theta _ { C }, \Theta _ { N C }, \boldsymbol \theta, \boldsymbol \eta _ { C, i }, N _ { i }, Z _ { i } ) = \theta _ { N, j } ^ { \dagger } - \theta _ { N, 0 } ^ { \dagger },$$

where

$$\theta _ { N, j, k } ^ { \dagger } = \log \left ( \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 | Z, k } ( Z _ { i } ) } } \right ).$$

Proof. Similar to Theorem 1, consider the change of variables θ † N, 0 = θ N , θ † N,j = θ N,j + θ N , and θ † C,j ( Z i ) = θ C,j + ψ N ( θ N,j + θ N + θ N Z | ( Z i )) -ψ N ( θ N + θ N Z | ( Z i )) . We may express our optimization in this alternative parameterization as

$$\mathcal { L } _ { i } ^ { \dagger } ( \{ \theta _ { N, j } ^ { \dagger } \} _ { j = 0 } ^ { m _ { C } } ) = \sum _ { j = 0 } ^ { m _ { C } } \eta _ { C, i, j } ( N _ { i } \cdot \theta _ { N, j } ^ { \dagger } - \psi _ { N } ( \theta _ { N, j } ^ { \dagger } + \theta _ { N | Z } ( Z _ { i } ) ) )$$

$$+ N _ { i } \cdot \theta _ { N | Z } ( Z _ { i } ) + \theta _ { C } ^ { \dagger } ( Z _ { i } ) \cdot \eta _ { C, i } - \psi _ { C } ( \theta _ { C } ^ { \dagger } ( Z _ { i } ) ),$$

such that max θ N , Θ NC ∑ m S i =1 L = max { θ † N,j } mC j =0 ∑ m S i =1 L † i . Moreover, the log-partition function of M N is given by ψ N ( θ N ) = ∑ m N k =1 e θ N,k and we may therefore conclude that

$$\cdots _ { N } \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \alpha \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \lambda \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \chi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \eta \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \Xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xeta \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi \xi
 \cdots _ { N - 1 }$$

$$\nu _ { N, j, k }$$

This is an arg max variant of a Legendre transform and its solution is log ( ∑ mS i =1 η C,i,j N i,k ∑ mS i =1 η C,i,j e θ N Z,k | ( Z i ) ) .

## 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 θ N and Θ 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 Z = [0 , π ] , which represents e.g. the orientation of a grating, as is common in V1 experiments. We define the natural parameter function by θ N Z | ( z ) = Θ NZ · s Z ( z ) , where s Z ( z ) = (cos(2 z , ) sin(2 z )) . This choice of θ N Z | allows us to express the conditional firing rates of the model neurons as E P [ N k | Z = z, C = ] = j γ j,k f k ( z ) , where γ j,k is the gain and f k ( z ) ∝ e ρ k cos( x -µ k ) is a von Mises density function with preferred stimulus µ k and precision ρ k . The dependence of the model on the categorical variable is expressed solely through

Figure 2: Training a CMP model on synthetic data generated from a ground truth CMP with m N = 20 neurons and m C = 7 mixture components. ( A ) Tuning curves of the ground truth CMP. ( B ) Stimulusdependent mixture weights of the ground truth CMP, with two weight-curves highlighted in black. ( C ) Best of 10 descents of the negative log-likelihood by the EM (red) SGD (blue) and Hybrid (purple) algorithms, with lower- and upper-bounds given by the ground truth CMP and the 1-component model, respectively. ( D ) Weights of the CMP learned by the hybrid algorithm, with two weight-curves highlighted that are qualitatively similar to those in (B). ( E ) Stimulus-dependent, paired correlation matrices of the ground truth and hybrid-trained models. Upper-right and lower-left triangles are correlations of the ground truth and learned CMP, respectively. Bottom, middle, and top matrices are conditioned on z = 0 π z , = 0 33 . π , and z = 0 67 . π , respectively, matching axis ticks in (D).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00637v2-figure-e45e01c2-496b-4843-8842-4b5c13d72b5b.jpg)

gain components γ j = ( γ j,k ) m N k =1 , and where E P [ N k | 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 Z | with m N = 20 neurons and m C = 7 mixture components. The preferred stimuli µ k of the tuning curves are drawn from a uniform distribution over Z , and the precisions ρ k and gains γ 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 Z | ( j | 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 , Z i i ) } m S i =1 from P N 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

Figure 3: Training a CMP model on synthetic data generated from a ground truth CMP with m N = 200 neurons. ( A ) Negative log-likelihood descents as per figure 2C. ( B &amp; C ) Mixture weights of the ground truth CMP (B) with two weight-curves matched qualitatively with weight-curves from the hybrid-train CMP model (C).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00637v2-figure-a5a1cde1-a8d8-418f-857d-2ebc37e24ca9.jpg)

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 -∑ m S i =1 log p N Z | ( N i | Z i ) in Figure 2C to establish a training lower-bound. Since a CMP Q 0 N 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 -∑ m S i =1 log q 0 N Z | ( N i | Z i ) for the optimal Q 0 N 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 Z | ( j | 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 stimuli 2 , 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.

2 Each 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].

Figure 4: Training a CMP model on response recordings from macaque V1. ( A ) 10-fold crossvalidated log-likelihood of 8 datasets as a function of the number of mixture components, relative to the 1-component log-likelihood for that dataset. We highlight the curve of a selected dataset which underlies the rest of the plots in this figure. ( B ) Negative log-likelihood descents on the full selected dataset, as per figure 2C. ( C ) Stimulus-dependent mixture weights of the CMP learned on the selected dataset trained with the hybrid algorithm, with tick marks indicating the conditioning angle of the correlation matrices in (D). ( D ) Paired correlation matrices, where the upper-right and lower-left triangles are the empirical and learned CMP correlations, respectively. Matrices are conditioned on x = 0 1 . π x , = 0 35 . π x , = 0 65 . π , and x = 0 85 . π , respectively, matching axis labels in (C).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00637v2-figure-8789d3a9-5df0-4b9b-8d87-0c9fe39ed01d.jpg)

In Figure 4C we plot q C Z | ( j | z ) for the hybrid-trained model Q CN 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.

## References

| [1]   | Bruno B. Averbeck, Peter E. Latham, and Alexandre Pouget. Neural correlations, population coding and computation. Nature Reviews Neuroscience , 7(5):358-366, May 2006.                                                                                                                          |
|-------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [2]   | Adam Kohn, Ruben Coen-Cagli, Ingmar Kanitscheider, and Alexandre Pouget. Correlations and Neuronal Population Information. Annual Review of Neuroscience , 39(1):237-256, July 2016.                                                                                                             |
| [3]   | Larry F. Abbott and Peter Dayan. The effect of correlated variability on the accuracy of a population code. Neural computation , 11(1):91-101, 1999.                                                                                                                                             |
| [4]   | Haim Sompolinsky, Hyoungsoo Yoon, Kukjin Kang, and Maoz Shamir. Population coding in neuronal systems with correlated noise. Physical Review E , 64(5), October 2001.                                                                                                                            |
| [5]   | Maoz Shamir and Haim Sompolinsky. Implications of neuronal diversity on population coding. Neural computation , 18(8):1951-1986, 2006.                                                                                                                                                           |
| [6]   | Alexander S. Ecker, Philipp Berens, Andreas S. Tolias, and Matthias Bethge. The Effect of Noise Correlations in Populations of Diversely Tuned Neurons. Journal of Neuroscience , 31(40):14272-14283, October 2011.                                                                              |
| [7]   | Rubén Moreno-Bote, Jeffrey Beck, Ingmar Kanitscheider, Xaq Pitkow, Peter Latham, and Alexandre Pouget. Information-limiting correlations. Nature Neuroscience , 17(10):1410-1417, October 2014.                                                                                                  |
| [8]   | Elad Schneidman. Towards the design principles of neural population codes. Current Opinion in Neurobi- ology , 37:133-140, April 2016.                                                                                                                                                           |
| [9]   | Christophe Gardella, Olivier Marre, and Thierry Mora. Modeling the Correlated Activity of Neural Populations: A Review. Neural Computation , 31(2):233-269, December 2018.                                                                                                                       |
| [10]  | Elad Schneidman, Michael J. Berry, Ronen Segev, and William Bialek. Weak pairwise correlations imply strongly correlated network states in a neural population. Nature , 440(7087):1007-1012, April 2006.                                                                                        |
| [11]  | Einat Granot-Atedgi, Gašper Tkaˇik, c Ronen Segev, and Elad Schneidman. Stimulus-dependent Maximum Entropy Models of Neural Population Codes. PLOS Computational Biology , 9(3):e1002922, March 2013.                                                                                            |
| [12]  | Benjamin R. Cowley, Matthew A. Smith, Adam Kohn, and Byron M. Yu. Stimulus-Driven Population Activity Patterns in Macaque Primary Visual Cortex. PLOS Computational Biology , 12(12):e1005185, December 2016.                                                                                    |
| [13]  | Alexander S. Ecker, Philipp Berens, R. James Cotton, Manivannan Subramaniyan, George H. Denfield, Cathryn R. Cadwell, Stelios M. Smirnakis, Matthias Bethge, and Andreas S. Tolias. State Dependence of Noise Correlations in Macaque Primary Visual Cortex. Neuron , 82(1):235-248, April 2014. |
| [14]  | John P Cunningham and Byron MYu. Dimensionality reduction for large-scale neural recordings. Nature Neuroscience , 17(11):1500-1509, November 2014.                                                                                                                                              |
| [15]  | Robbe L. T. Goris, J. Anthony Movshon, and Eero P. Simoncelli. Partitioning neuronal variability. Nature Neuroscience , 17(6):858-865, June 2014.                                                                                                                                                |
| [16]  | Michael Okun, Nicholas A. Steinmetz, Lee Cossell, M. Florencia Iacaruso, Ho Ko, Péter Barthó, Tirin Moore, Sonja B. Hofer, Thomas D. Mrsic-Flogel, Matteo Carandini, and Kenneth D. Harris. Diverse coupling of neurons to populations in sensory cortex. Nature , 521(7553):511-515, May 2015.  |
| [17]  | Robert Rosenbaum, Matthew A. Smith, Adam Kohn, Jonathan E. Rubin, and Brent Doiron. The spatial structure of correlated neuronal variability. Nature Neuroscience , 20(1):107-114, January 2017.                                                                                                 |
| [18]  | E. Ganmor, R. Segev, and E. Schneidman. Sparse low-order interaction network underlies a highly correlated and learnable neural population code. Proceedings of the National Academy of Sciences , 108(23):9679-9684, June 2011.                                                                 |
| [19]  | Dimitris Karlis and Loukia Meligkotsidou. Finite mixtures of multivariate Poisson distributions with application. Journal of Statistical Planning and Inference , 137(6):1942-1960, June 2007.                                                                                                   |
| [20]  | Wei Ji Ma, Jeff Beck, Peter Latham, and Alexandre Pouget. Bayesian inference with probabilistic population codes. Nature Neuroscience , 9(11):1432-1438, October 2006.                                                                                                                           |
| [21]  | Srdjan Ostojic and Nicolas Brunel. From Spiking Neuron Models to Linear-Nonlinear Models. PLoS                                                                                                                                                                                                   |

| [22]   | Shun-ichi Amari and Hiroshi Nagaoka. Methods of information geometry , volume 191. American Mathematical Soc., 2007.                                                                                                                                                        |
|--------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [23]   | Martin J. Wainwright and Michael I. Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends R © in Machine Learning , 1(1-2):1-305, 2008.                                                                                         |
| [24]   | Max Welling, Michal Rosen-zvi, and Geoffrey E Hinton. Exponential Family Harmoniums with an Application to Information Retrieval. In L. K. Saul, Y. Weiss, and L. Bottou, editors, Advances in Neural Information Processing Systems 17 , pages 1481-1488. MIT Press, 2005. |
| [25]   | Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 , 2014.                                                                                                                                                            |
| [26]   | Ruben Coen-Cagli, Adam Kohn, and Odelia Schwartz. Flexible gating of contextual influences in natural vision. Nature Neuroscience , 18(11):1648-1655, October 2015.                                                                                                         |
| [27]   | Geoffrey J. McLachlan, Sharon X. Lee, and Suren I. Rathnayake. Finite Mixture Models. Annual Review of Statistics and Its Application , 6(1):355-378, 2019.                                                                                                                 |
| [28]   | Daniel L. K. Yamins and James J. DiCarlo. Using goal-driven deep learning models to understand sensory cortex. Nature Neuroscience , 19(3):356-365, March 2016.                                                                                                             |
Sacha Sokoloski Ruben Coen-Cagli
Conditional Finite Mixtures of Poisson Distributions for Context-Dependent Neural Correlations
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.