Estimating Noise Correlations Across Continuous Conditions With Wishart Processes
Authors: Amin Nejatbakhsh
Isabel Garon
Alex H Williams
Venue: NeurIPS 2023 poster (2023)
Keywords: Noise Correlations
Wishart Process
Variational Inference
TL;DR: We develop a Bayesian modeling approach based on Wishart processes to accurately estimate noise correlations with very few trials by exploiting the smoothness of covariance change across neighboring conditions.

Abstract: The signaling capacity of a neural population depends on the scale and orientation of its covariance across trials. Estimating this "noise" covariance is challenging and is thought to require a large number of stereotyped trials. New approaches are therefore needed to interrogate the structure of neural noise across rich, naturalistic behaviors and sensory experiences, with few trials per condition. Here, we exploit the fact that conditions are smoothly parameterized in many experiments and leverage Wishart process models to pool statistical power from trials in neighboring conditions. We demonstrate that these models perform favorably on experimental data from the mouse visual cortex and monkey motor cortex relative to standard covariance estimators. Moreover, they produce smooth estimates of covariance as a function of stimulus parameters, enabling estimates of noise correlations in entirely unseen conditions as well as continuous estimates of Fisher information&mdash;a commonly used measure of signal fidelity. Together, our results suggest that Wishart processes are broadly applicable tools for quantification and uncertainty estimation of noise correlations in trial-limited regimes, paving the way toward understanding the role of noise in complex neural computations and behavior.

Estimating Noise Correlations Across Continuous Conditions With Wishart Processes
Amin Nejatbakhsh Isabel Garon Alex H Williams Center for Neural Science, New York University, New York, NY Center for Computational Neuroscience, Flatiron Institute, New York, NY jtbakhshh,, igaron, awill
org
Abstract
The signaling capacity of a neural population depends on the scale and orientation its covariance across trials. Estimating this "noise" covariance challenging and is thought to require large number of stereotyped trials. New approaches are therefore needed to interrogate the structure neural noise across rich, naturalistic behaviors and sensory experiences, with few trials per condition. Here, we exploit the fact that conditions are smoothly parameterized in many experiments and leverage Wishart process models to pool statistical power from trials in neighboring conditions. demonstrate that these models perform favorably on experimental data from the mouse visual cortex and monkey motor cortex relative standard covariance estimators. Moreover, they produce smooth estimates covariance function of stimulus parameters, enabling estimates of noise correlations in entirely unseen conditions as well as continuous estimates of Fisher information commonly used measure of signal fidelity. Together, results suggest that Wishart processes are broadly applicable tools for quantification and uncertainty estimation of noise correlations in trial-limited regimes, paving the way toward understanding the role of noise in complex neural computations and behavior.
Introduction
Nominally identical repetitions of sensory stimulus or behavioral action often elicit variable responses neurons. Intuitively, this "noise" degrades the information capacity of neural representa- tions. However, the precise nature of this effect is complex since neural responses are not independent random events. Indeed, pairs of co-recorded neurons exhibit significant "noise correlations" from trial-to-trial, with well-studied and intriguing consequences [1-5]. In particular, introducing noise correlations into a circuit can theoretically degrade or enhance its signaling capacity depending on the geometric configuration of neural representations in firing rate space (Fig. IA; [6. 7]).
Estimating noise covariance is therefore an important problem that arises frequently in neural data analysis.' However, this problem is widely regarded as challenging [8-10] because the number of estimated parameters (i.e., number of unique neuron pairs) grows quadratically with the number of co-recorded neurons, while the number of measurements per trial (i.e., number of neurons with a response) grows only linearly (for more discussion see, e.g., [11]). As result, experimental investigations into noise correlation have been limited to extremely simple tasks where a very large number trials can be collected over a small number of conditions. For example, Rumyantsev et al. [10] recorded neural responses over just two stimulus conditions (oriented visual gratings), repeated across hundreds of trials. We therefore have rudimentary characterization of how noise correlations
'Throughout, we focus on the problem of estimating noise covariance instead of noise correlation. By accurately estimating the former. can also estimate the latter since the correlation between two neurons is given by their covariance normalized each neuron's marginal variance.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Information Limiling Noise Covariance
Isotropic
E
Noise Covarinace
Prior Sumples
U.ub Pc.P)
Maias Respanes
Med
Incrensting Smoothmen
F
ULPo Po.P1
Tedaperal Premjency
Radiisx
Max
A/a
Angle
Orientation
Incrensing Mean Smoothness
Figure I: (A) Illustration of information limiting noise correlations. (B) Experimental datasets with smoothly parametrized conditions (see [7]). Left, a nonhuman primate makes point-to-point reaches radial targets. Right, rodent views drifting visual grating stimuli of different orientations and temporal frequencies (i.e. speed). (C) Parameterized condition spaces for the datasets in B. (D) Samples from Gaussian and Wishart process prior distributions for with C 6 conditions over 1D periodic variable (e.g. angle) and N 2 neurons at varying kernel smoothness parameters (A in eq. 6). Dots and ellipses represent the mean and covariance of neural responses. Colors appear in the increasing order of condition value (e.g. angle) and dashed lines connect neighboring conditions. Increasing A (horizontal axis) encourages smoothness the axis) encourages the ellipses to change smoothly in scale and orientation. (E-F) Graphical model of the Wishart model with Normal (E) and Poisson (F) observations (see Supplement B.5 for details of the Poisson model).
while increasing A (vertical
change across rich or naturalistic experimental conditions, despite interest among the theoretical community [5, 12, 13), and experimental evidence that noise statistics context-dependent [14- 16]. New statistical approaches that are capable of inferring noise covariance with fewer trials per condition are needed make progress this subject.
To address this need, we develop probabilistic model based on Wishart processes [17]. The core insight exploit is that similar experimental conditions- e.g. visual gratings at similar orientations or cued arm reaches to similar locations-ord ought to exhibit similar noise statistics. Thus, by building this insight into the prior distribution of probabilistic model, can pool statistical power across neighboring conditions and accurately infer noise correlation structure with fewer trials per condition. In fact, this principle also enables by interpolating between (and potentially extrapolating from) observed conditions. We demonstrate the scientific potential and generality of this approach through two experimental applications spanning different modalities (vision and motor tasks) and species (mouse and nonhuman primate).
predict covariance in entirely unseen conditions
2 Methods
2.1 Problem Setup
Consider an experiment with N neurons recorded over C' distinct stimulus conditions, each with K repeated trials. Experimental advances continually enable to record larger neural populations (large N) and neuroscientists are inereasingly interested performing rich and naturalistic experi- ments (large C). Both of these trends put us in regime where K < N, meaning there are few trials per condition relative to the number of neurons [9].
"To simplify exposition, we assume that there are equal number of trials in each condition. However, the model easily applies to the more general with unbalanced trial allocations.

---

Empirical Estimator: Our goal is to estimate the N N covariance matrix, Σ., describing trial-by-trial variability in each condition E (1....,Cy C}. If yeu E R the measured response on trial k of condition then straightforward estimator is the empirical covariance:
E = (1/K) Σ cyck [yCk - -.)" where
V. = (1/K) Σ yok
(1)
This estimator behaves poorly in our regime of interest when K < N. Indeed, while the true covariance must be full rank, the empirical covariance is singular with N K zero eigenvalues.
Regularized Estimators: Statisticians often refer to problems where K < N as the high- dimensional regime. Estimating covariance in this regime is actively studied and understood to be challenging [18]. The Ledoit-Wolf estimator [19] is a popular approach, which regularizes ("shrinks"" the empirical covariance towards a diagonal estimate. Other estimators shrink the individual co- variance estimates towards the grand covariance (see eq. 2 below) instead of a diagonal target [20]. Graphical LASSO [21] is another popular approach, which incorporates on the inverse covariance matrix. All of these methods can be understood as forms of regularization, which reduce the variance of the covariance estimate at the expense of introducing bias.
parrsit-promooinn penalty
Grand Empirical Estimator: The methods summarized above (or variants thereof) are common within the neuroscience literature [8, 22, 23]. It is also common to simply estimate the covariance as being constant across conditions [10, 24]. in which case an estimate can be found by pooling trials across all conditions to achieve an empirical estimate of the grand covariance:
Σ プー KC
Σ (yet R,C
(y - -..
This estimator can be thought of as having high bias (since covariance often thought to be stimulus- dependent [14-16]), but low variance since many more trials can be averaged over after pooling across conditions. We compare our model these baselines in section 3.
2.2 Wishart Process Model of Noise Correlations
Our insight is to incorporate alternative form of regularization whenever the experimental condi- tions are smoothly parameterized. For instance, a classical stimulus set used by visual neuroscientists is grating pattern varied over C' orientations. In fact, such stimuli are commonly used for noise correlation estimation [8, 10, 25, 26). Intuitively, two gratings with nearly identical orientations should evoke similar probabilistic responses. Thus, as stimulus orientation incremented, the mean and covariance of the response should smoothly change along a manifold (Fig. IA).
This idea readily extends to other experimental settings. For example, pure acoustic tones with smoothly varied frequencies are often used in auditory neuroscience [27]. In studies of olfactory coding. isolated odorants or mixtures of odorants can be presented at smoothly varied concentra- tions [28]. In motor neuroscience, voluntary reach movements can be cued to continuously varied target locations [29]. Finally, in the context of navigation, noise correlations have also been studied in neural ensembles coding for continuous variables like position [30] and head direction [31]. In principle, the modeling approach
develop is applicable any of these experimental settings. This is in contrast to experimental settings that employ a discrete set of stimuli with no underlying continuous parameterization (such as olfactory responses to distinct classes of odors).
To formalize these concepts, let c E {l..... C} index the set of conditions, and let x, E X denote the stimulus or animal behavior in condition c. For example, a. could be scalar denoting the orientation of a voluntary reach or a visual grating (Fig. IB). However, the model we will describe is more general than this: x. could be vector describing a multi-dimensional stimulus space. For example, we will later study motor preparatory activity to reach targets at different radii and angles (sec. 3.2). In this setting, æ. is 2-dimensional vector specifying these two parameters of the reach target (Fig. 1C, left). We will also study neural responses to drifting visual gratings (sec. 3.3). In this setting, H. is vector specifying the speed and orientation of the grating stimulus (Fig. 1C, right).
Recall that Yok E N denotes the simultaneously measured response of N neurons to condition CE {L...., C} on trial & E {L,...,K}. Our goal is to model these measured responses as samples
3

---

from some distribution F with mean and covariance Σ:
Yck ~ F( ((te)) independently for c, k.k E (T, KJ..
(3)
RR Importantly, the functions μ X R. and : R should be smooth so that the response distribution changes gradually as the experimental condition altered. To accomplish this, we leverage the machinery of Gaussian processes, which can be used to specify a prior distribution over smooth functions [32]. While Gaussian process models have been extensively used within statistical neuroscience [33-37], they have typically been used to model smoothness in temporal dynamics rather than smoothness across experimental conditions (but see [38]).
In this paper, we assume that F(U, Σ) is multivariate Gaussian distribution, which enables tractable and convenient probabilisticc inference (see sec. 2.4). Since many noise correlation analyses are applied to spike count data, it is of great interest to extend our model so that F is better suited to this case. Multivariate extensions of the Poisson distribution are an attractive option [39]. While the focus of this paper is multivariate Gaussian distributions, we develop inference tools for the Poisson model and we leave further investigation this model for future work (see IF for the corresponding graphical model and Supplement B.5 for further details). We note that it is common to pre-process spike count data and calcium fluorescence traces (e.g. through square root or log transforms) so that they are better approximated by Gaussian models [33, 40, 41].
To model px)),, we stack independent Gaussian processes into an -evector. To model E()), we use Wishart process [42]. In essence, the Wishart process (denoted by WP) stacks Gaussian processes into an N x P matrix and uses this as an overcomplete basis for the covariance. Formally,
1 ~ GpN (O kw), U() ~ NP ks),
E(x) = Ux))
where P 2 N is a hyperparameter, and GPN (m, k)) denotes an N N-dimensional Gaussian process with (constant) mean m IRN and positive definite kernel function k: : X' x ~ R, described below. For any fixed value of æ E X. the product U (x))U()) standard probabilistic model for covariance and inverse covariance (see sec. 3.6 of [43]). To our knowledge, Wishart processes have not previously been leveraged for neural data analysis.
follows Wishart distribution, which
The kernel functions kµ and ky play important roles by respectively specifying the degree of smoothness in the mean and covariance of the neural response across conditions (Fig. ID). The kernel functions must satisfy certain technical conditions, as reviewed in [44, 45]. We satisfy these conditions by defining a kernel for each coordinate of x E X' and evaluating the product, k(x, x²) II, (x,, x This is guaranteed produce valid kernel ([44], Proposition 13.2). For non-periodic variables, we use a standard squared exponential kernel
k(x, x') = 700z' + Bexp [-(x - ---xx'''' /AA
(5)
where (,,3,A] }} are hyperparameters and d = when x = x' and otherwise evaluates to zero. For periodic conditions (e.g. orientation of visual gratings or radial reach directions) use [46]:
k(x,x') = 762m' + 8 exp sin (I 1 =///T//J
(6)
where specifies the period of the kernel. For both kernels, the most critical hyperparameter is X (often called the "bandwidth"), which specifies the degree of smoothness across conditions. The remaining hyperparameters, Y and 3, respectively ensure that the Gaussian process is well- conditioned and scaled. We tune all of these hyperparameters to maximize heldout log likelihood by cross-validation. Throughout, we use 177,,,m )n) and {Tx, fx, >s:} to respectively represent the kernel hyperparameters for kn and ky.
2.3 Low-Rank Plus Diagonal Covariance and Other Model Variants
The vanilla Wishart process model in equations (3) and (4) can be extended in several ways, which we found to be useful in practical applications. First, we augment the covariance model by adding a diagonal term, resulting in a model that is reminiscent of factor analysis:
{xx) U())U(()) + A(#),
where (ax)
N condition-dependent nonnegative diagonal matrix where the diago- nal terms are independent samples from a GPN that transformed using a softplus function,

---

= + +e), to ensure positivity. We term this extension of WP as Wpird (for low-rank plus diagonal), Intuitively, this model extension has fewer parameters than a vamilla WP and hence discourages overfitting. Additionally, the diagonal term, (ax), helps ensure that the covariance matrices remain well-conditioned during inference. Indeed, we can choose values of P < N (or even P = 0), which enables to specify prior distributions over the covariance that are "low-dimensional" in the sense that a large fraction of variance can be explained by a small number of dimensions. This comports with numerous experimental observations [10, 47]. As can be seen in Supplement Fig. 5, the dimensionality of matrices sampled from a vanilla Wishart distribution scales linearly with N, while matrices sampled according to equation (7) saturate at lower values.
Next, we additionally incorporate lower- teier-trianglar matrix R RNN into the covariance model to capture aspects of covariance that are condition condition-indepenee independent (i.e. do not depend on æ):
E(x) L(U (x) + +A(x))L",
(8)
In the literature LL" is referred as the scale matrix. Some intuition can be gained by assuming (æx = O and P > N, in which case U (x) x I in expectation under the marginal prior distribution. Thus, a good choice is to set to be the Cholesky factor of the grand empirical covariance as defined in eq. (2). In our experiments, we use this as an initialization for L and we optimize to maximize the evidence lower bound (ELBO) as described in the Supplement B.3.
It is also possible, and potentially beneficial, to instead model smoothness in the inverse covariance matrix (also known as the precision matrix). This results in an inverse Wishart process. We refer the reader to Supplement B.2 for additional details of the inverse Wishart model.
2.4 Latent Variable and Parameter Inference
Each of the probabilistic models above specifies posterior distribution of the noise covariance con- ditioned on observed neural activity, { [E(=:)) (Vok ) c=1,,=1)... More precisely, each model
specifies a posterior on the latent covariance factor matrices, p (UU())) 1 which defines the desired posterior over covariance through measurable transformation (eq. 8). This distribution over latent covariance factors has NPC dimensions, precluding us from repre- senting the posterior exactly (doing so would require us to tile R.PC with a lattice of points). We instead use mean field variational inference [48] to approximate the desired posterior. Since this approach has already been described for similar Wishart process models [49]. we relegate description to Supplement B.3. All codes are implemented in numpyro [50] and available at https:// gigtt.u..... www.coom//eene http://eoossaalla/iisseeeeeeeeee wisaraabb/iiiaatt-pressss. -proces...
2.5 Connections to Population Decoding
As we discussed in section 1, the estimation of noise covariance matrices is an important step for a variety of downstream neural analyses. Here, we briefly highlight two applications of Wishart process models in the context of neural decoding and perceptual discrimination [3].
Quadratic Discriminant Analysis (QDA) A simple decoding task is infer which condition e E [L....,C} gave rise to an observed neural response, y E R on a single trial. A natural approach is to use a maximum likelihood classifier, class(y) = arg max.c(L....C) log Ex). Doing so requires that we construct estimates for the mean responses, P1...., Pc. and covariances. E,,..., Ec. Since covariance estimation is challenging, a common practice (see e.g. [10]) is to use the grand empirical covariance (eq. 2) and therefore set Σ, = Σc. Under this model, known linear discriminant analysis (LDA; see e.g. [51] Ch. 4), it can be shown that the decision boundaries of the maximum likelihood classifier become linear.
Intuitively, if we can accurately estimate distinct covariances for each condition, we may unlock more powerful decoders by relaxing the assumption of equal covariance across conditions. Under this model, known as quadratic discriminant analysis (QDA), the decision boundaries are described by hyperbolas (A schematic is shown in Fig. 2D). Furthermore, Pagan, Simoncelli, and Rust [52] described a circuit implementation of QDA using cascade of two iineir-nonlinear models, suggesting that such classifiers are simple and biologically plausible neural computation. We demonstrate the ability of Wishart process models to enable accurate decoding via QDA in Fig. 2E.
5

---

Grand-FFpiiic Empiriaal
Empiriaa
Tex
-10
10
PC1
PC1
PC1
PC1
PC1
50
Log Probobaltt
-500
Erpetual
Expetaal
-200
-80
-750
-300
-8
92
40
40-
00
24
Openiur
30
32
45
L6
20
24-
30
8
16
15
20
40
80
40
60
80
1.5
3.0
4.5
Dimension
Conditions
Smicefheess
Quadratic
1.00
Pisher Information
Fivbur Information
1.6
Diseriminant
Aceuray
0.75
Maste PoLice act
Analysis
0.50
1.2
0.25
0.8
0.0
1.5
3.0
6.0
150
300
Smoothheess
Stimalus Orientation
Figure 2: Results on synthetic data. See Supplement B. for all simulation details. (A) Simulated (AI) and inferred (A2-A5) response distributions using different methods. Mcans ((x,) and covariances E(x,) are shown in dashed lines and colored ellipses in PC space. Train and test trials you shown in black and red dots. (B1-B4) Mean log- likeliidl of held-out trials as a function of the number of dimensions, conditions, trials, and covariance smoothness. Error bars represent median + SEM across 5 independent runs. (C1-C4) Same as in B, the operator norm between the true and inferred covariances. (D) Schematie of QDA classifier. (E) Using different estimates of covariance, we use QDA to infer the condition label of held-out trials. Error bars represent the standard deviation over independent runs. (F) Schematic of Fisher Information, FI. (G) Estimation FI using grand empirical covariances and Wishart model.
Fisher Information (FI) A quantity that plays an important role in our understanding of how neural noise properties affect behavioral discriminability and representation is FI defined as E[ă² log p(z)((())'') [3]. Unfortunately, estimating
challenging because it involves mea- suring infinitesimal changes in the neural response distribution with respect to the stimulus condition (Fig. 2F). Wishart processes enable
exploit the continuous parametrization of experimental conditions to derive natural estimator for FI. For constant covariances across conditions, the FI for multivariate normal distributions is given by p p²(x)" -<>=)>>>>''" (1'()) where μu' (xx is the gradient of the mean with respect to If covariances vary across conditions, this quantity is replaced with μí "=>) (x) + 1 str(E" the conditions and tr is the trace operator. Wishart process models enable the continuous estimation
E')))()>)))) where E' E'xx) is the gradient of the covariance w.r.t.
of FI as function of z. This is achieved by sampling from the gradient of the inferred posterior distribution (see Supplement B.6).
3 Results
We qualitatively and quantitatively compared Wishart processes to common covariance estimation methods used in neural data analysis: Ledoit-Wool [19], empirical covariance (eq. 1), grand empirical covariance (eq. 2), and graphical LASSO [21]. Our primary measure of performance is the log-

---

likelihood that the models assign to held-out data. Because the baseline methods are not Bayesian, we handicap ourselves by using a single sample from the posterior distribution to calculate held-out log- likelihoods. We also visualize and qualitatively evaluate contours of the predictive probability density ("covariance elllises"") the top two principal component subspace ("PC space"). For synthetic data, we report the difference between inferred and true covariances measured in the operator norm (Fig. 2C1-C4). We observed that generally log-likelihood and operator norm results are aligned. However, since they measure the performance in two different ways it is possible for a model to outperform another in log probability but vice versa on the covariance operator norm.
3.1 Synthetic Data
We first validate our approach on simulated data over a periodic I-dimensional space, with conditions sampled at C' equispaced angles on the interval [0, 2T). The mean and covariance of neural responses were sampled according to the Wishart process generative model (eq. 4) and trials for each condition were conditionally sampled according to eq. (3) from multivariate Gaussian distribution. An example of generated data is shown in Fig. 2A1. Due to our choice for the condition space, we use periodic kernels both for the Gaussian and Wishart process priors. For simplicity, we set A (xx) to the identity matrix both in the data generative and inference models. To match experimental observations [47], we simulated low-dimensional noise covariance structure by setting scale matrix to LL = Udiagdiagssaa. 11 )U" where U is a uniformly sampled orthogonal matrix in N dimensions and 81 82 > ... > SN > 0 are logarithmically spaced numbers between I and 10-s. Except for the number of dimensions, conditions, trials, and Wishart smoothness (Az), which we vary over a range to assess the performance of different models, all other hyperparameters
kept constant (P = 2, = 0.001, = 1. and A = 1). Supplement B. provides further details for these simulations, covering all specific hyperparameter settings for each panel of Fig. 2.
We fit Wishart process models to these simulated datasets using the procedures described in sec- tion 2.4. These models qualitatively capture the ground truth covariance structure (Fig. 2A2-A5) and quantitatively outperform alternative methods across a range of synthetic datasets. Wishart process inference performs particularly well in nigh-dimensional regimes relevant to neural data analysis with large N, small K, and intermediate levels of smoothness across conditions (Fig. 2B-C). We also investigated the ability of Wishart process models to improve decoding accuracy on this dataset via QDA, as summarized in section 2.5. The Wishart process model outperformed the baseline models across range of data smoothness suggesting that it captures the boundaries between classes more accurately (Fig. 2E). Furthermore, in Fig. 2G we show that on synthetic data the Wishart process model accurately estimates the Fl continuously across conditions.
3.2 Primate Reaching
To study how the brain plans upcoming actions, animals are often trained to prepare voluntary movements to cued locations during "delay period" before executing those movements after "go cue" delivered [53-55]. During this delay period, neural firing rates in the primary motor cortex (M1) and dorsal premotor cortex (PMd) correlate with features of the upcoming movement (c.g. direction and distance) in nonhuman primates [56, 57J. Trial-to-trial variability in the delay period has been suggested to correlate with fluctuations in reaction time [58], but quantitative account of this variability with dense target configurations is challenging to achieve [29].
We analyzed 200ms of spiking activity prior to the onset of the "go cue" across N = 192 multiunit channels from both PMd and MI in a nonhuman primate executing point-to-point hand reaches. The dataset contains C 48 conditions, where targets were organized in three concentric circles radii of 4cm, Scm, and 12cm, with 16 equidistant targets around each ring (Fig. 3A). We used a squared exponential kernel (eq. 5) to model smoothness across reaches of different lengths and a prrpeoio kernel (eq. 6) to model smoothness across reach angles. Models were fitted on 60% of the available trials, radial and angular smoothness parameters for both GP and WP kernels were selected validation set of 25% of the data, and log-likelihooe scores were generated using the remaining 15% of trials. This procedure was repeated over 8 randomized data partitions.
The Wishart process model outperforms competing estimators (Fig. 3B). Since N > K, the empirical covariance (eq. 1) is singular and therefore produces negative infinity log- liglliiklssoods. Instead,
Data was previously published and made available by [29].
7

---

C
Example Neuren
-175
0.60
Log Probability
0.50
0.40
-180
Viriacee
0.30
Condition Space
100
100
Example Neuron 2
0.70
185
0.60
Kadrus
0.50
Sapirral
Wegdin And
LANSY
Lokaa
Angle
-100
Reaching Angle
Wishart
Empirical LW
Mean Response Projected onto lins PCS
2.5
O
-2.
-150
-100
-50
100
150
PC1
PCI
Reaching Angle
1.0 100
50 100 150
Renching Angle
Figure 3: Results on primate reaching task. (A) Task schematic (left) and set of experimental conditions (right). (B) Log-lite ligellikooood held-out trials across 8 randomized folds. Horizontal lines indicate the maximum, minimum, and mean of the lig-liieldhoo distribution. (C) Marginal variance of two example neurons across reaching angles, averaged over reach radii (green) and Wishart process model estimate (black). Shaded green denotes mean + SEM. (D) Wishart (leff) and Ledoit-Won Wolf covariance, projected onto the first two principal components. To
(right) estimates of mean demonstrate the ability of the Wishart Process model to predict covariance in unseen conditions, the middle ring of targets was held out in training (red ellipses). (E) Predicted mean and (F) covariance in held-out conditions shown in panel D. Covariance ellipses, grand-empirical estimates (dashed lines), and samples (black dots) are visualized in the top-2 principal component subspace.
weighted averages of the grand empirical and empirical covariance calculated according to y >n (x.) (x.) + (1 a) (WA, E, and GE denote weighted average. empirical, and grand-empirical estimators respectively). This baseline is simplified variant of [20] shrinking the per-condition covariances towards the grand-empirical estimate. Performance results produced by this baseline are shown in Fig. 3B at different values of a.
These results suggest that WP is able capture the smooth change in the covariance across conditions. Since smooth changes the covariance imply smooth changes in the marginal variance (diagonal terms in the covariance), we plot the marginal variance of individual neurons as function of the reaching angle to visualize the smooth change. We see that the variance indeed changes across conditions and that the Wishart process model provides a smooth fit to these noisy variance estimates (Fig. 3C). This smoothness is absent in the predictions of the Ledoit- Wolf and graphical LASSO estimators, while the grand covariance estimate is constant and, in essence, over-smoothed (Supp. Fig. 2). The -at-ol- log probability results suggest that the smoothness captured by WP is indeed property that is present in the data and has physiological relevance (Fig. 3B).
We also investigated the Wishart model's ability to interpolate covariance estimates for unseen conditions. Here, the model was retrained on only trials from the inner and outermost rings (4em and 12cm). The neural response means and covariances around the center ring were then estimated by sampling from the model posterior evaluated at these hold-out condition locations (8cm, all angles).

---

Condition Space
-54
Teriporal Freqdency
Log Probability
-60
600
66-
46
Brpprezd
72
72
Gral Expltial
Witter
-78
0.0
1.5
3.0
4.5
6.0
3
Orientation
Smoothness
Number of Components
Training Conditions
Unseen Test Conditions
Smprto Empirical LW
Wishart
Empirical + LW
Wishart
1.5
1.5
0,0
0.0
1.5
Contioon
Coudicc
Condition
Condition 25
Condition
Condition
2.5
0
0
-2.5
PC1
PC1
PC
PC1
Figure 4: Results on drifted gratings dataset. (A) Task schematic (left) and set of experimental conditions (right), (B) Cross-validated log probability as a function of Wishart smoothness (As) and the number of components (P) achieving their optimal values at Äg 1.5 and /' 0 and outperforming alternatives. (C) Top: mean responses and noise covariances of the full set of conditions plotted in the PC space for the Ledoit-Woo (left) vs. Wishart (right) estimators (covariances are scaled down by factor of 100 for visualization purposes). Colors are in the increasing order of condition number matching the colors (A); Bottom: Train and test trials (black and red dots) in two example conditions plotted in the top two PC subspace. Covariance ellipses are color-coded by the model as panel B. (D) Generalizing to unseen conditions: Top: 80% of the conditions are used for training (shown in colors) and the remaining conditions (shown in red) are withheld. Bottom: Inferred covariances by Wishart Process model for two held-out conditions (as in panel C, bottom). Notice that panels C, and D are from two separate analyses. In panel C, we train all conditions and test on held out trials. In panel D. we train on subset conditions and test on held out conditions (interpolation stimulus space). Condition was held out only in the second analysis.
Wishart model estimations perform comparably to empirical mean and grand-empirical covariance estimations calculated on these hold-out conditions (Fig. 3E-F).
3.3 Allen Brain Observatory (Neuropixels)
Perceptual discrimination is longstanding area of research in psychology and neuroscience, dating back to work in the early 1800s by Weber and Fechner [59]. A classic paradigm requires subjects to discriminate between pairs of oriented visual gratings, which humans and nonhuman primates accomplish with high precision [60-62]. Discrimination thresholds on this task are thought to be related to the scale and orientation of noise neural representations [7] (Fig. IA), and recent works have examined this hypothesis by simultancously recording large populations of neurons in rodent models [10, 25, 26]. The resulting analyses are delicate and require large number of trials.
To investigate, we analyzed simultaneously recorded responses in mouse primary visual cortex to drifting visual gratings in the Allen Brain Observatory (Visual Coding: Neuropixels dataset). 4 We analyzed a session (session-id 756029989) with 81 neurons. Spike count responses in a 200 ms
'nttps://portal wratoertta..aaanneeeee brain- www.aa......aa.e. wwwprr/paaeeeeatttttt/////tuueea neuroiiees

---

window from the stimulus onset were recorded across 4 grating orientations, 2 directions, and 5 temporal frequencies as illustrated in Fig. 4A. We used periodic kernels for the orientation and direction axes (eq. 6), and a squared exponential kernel (eq. 5) for the temporal frequency axis of the condition space. For each of the resulting 40 conditions, 15 trials are collected, which we split into 13 train vs. 2 test trials. For choosing the number of components (P) and the smoothness of the Wishart kernel (Ag), we quantified the cross-validated log probability of test trials over 10 randomized folds. As depicted in Fig. 4B the optimal performance of the Wishart model is achieved for Xy R 1.5 and P=0.
This selection of the model parameters allows us to capture the structure that is shared among all conditions (note the similarity between Wishart and grand empirical covariance in Fig. 4C), while representing the condition anddition-specffic specific changes in the covariance structure (note the smooth change of matrix orientation across conditions in the top-right plot of Fig. 4C). Finally, we present covariance interpolation results across unseen conditions where we fit the model according to 32 randomly selected train conditions (colored ellipses in Fig 4D top-left) and leave out the complete set of trials for the remaining 8 conditions (red ellipses in Fig 4D top-right). The interpolated means and covariances for two example test conditions present a strong agreement with their corresponding held-out trials (Fig 4D bottom row).
4 Conclusion
We proposed a method to estimate covariance across large populations of co-recorded neurons, which is widely regarded as a challenging problem [8, 9] with important implications for theories of neural coding [3. 5, 7]. Specifically, we formulated probabilistic Gaussian and Wishart process models [42] to model smoothly changing covariance structure across continuously parameterized experimental conditions. We found this approach outperforms off-the-shelf covariance estimators that do not exploit this smoothness (e.g. Ledoit-Wolf) as well the grand covariance estimator (eq. 2) which is infinitely smooth (i.e. does not allow for condition-dependent covariance). We investigated two task modalities (vision and motor) and model organisms (mouse and nonhuman primate). However, we expect these modeling principles will be even more broadly applicable, including to other sensory modalities with smoothy parameterized stimuli (c.g. acoustic tones or odorant mixtures) and in neural circuits supporting navigation [30, 31].
We believe our work also provides an important conceptual advance. Given the limited duration of experiments, it may seem that there is a strict tradeoff between accurately estimating neural responses across many experimental conditions and estimating trial-by-trial noise statisties in any fixed condition. Our results show that this tradeoff is not so simple. Indeed, if neural responses are sufficiently smooth across conditions, the Wishart process model could predict neural covariance accurately by densely sampling conditions with only single trial per condition (e.g. visual gratings random orientations, as done in [25]). Intuitively, the model works by pooling information across nearby conditions, and so can even predict covariance in entirely unseen conditions (as demonstrated in Figs. 3D and 4D). Thus, the limiting factor in estimating noise covariance not the number of trials per condition, but the smoothness of the neural responses over the space of experimental conditions.
We conclude by discussing the current limitations of our model and opportunities for future work. First, the model requires us to explicitly specify kernel functions (ky,, ks:) over the space of experi- mental conditions. We studied simple paradigms where reasonable parametric form of this kernel could be specified and the hyperparameters could be fine-tuned through cross-validation. Extending the model to more complex stimulus sets (e.g. to naturalistic image stimuli) would be of great interest, and likely require more sophisticated approaches for learning kernels and feature representations of these stimulus sets. Second, we assumed that trial-to-trial noise followed Gaussian distribution. A deeper investigation of noise models, including multivariate extensions of Poisson distributions [39, 47, 63], is warranted. Finally, we utilized a basic variational inference procedure to approximate the posterior distribution over neural covariance. While fully representing the posterior is intractable (as is typically the case), more sophisticated optimization techniques may further benefit performance [64]. Nonetheless, our results show that Wishart process models are performant even in the absence of the aforementioned extensions. These models outperform common practices in neuroscience and open the door to new analyses, such as the prediction of noise covariance in held-out conditions.
10

---

Acknowledgements
We thank Saurabh Vyas (Columbia) and Eric Trautmann (Columbia) for helpful discussions and for providing data from the nonhuman primate motor cortex. We also thank Sam Zheng (NYU) and Sarah Harvey (NYU) for their feedback on the manuscript.
References
[1] E Zohary, M Shadlen, and W T Newsome, "Correlated neuronal discharge rate and its implications for psychophysical performance". In: Nature 370.6485 (July 1994), pp. 140-143. [2] M N. Shadlen, K H Britten, W ' Newsome, and J A Movshon. "A computational analysis of the relationship between neuronal and behavioral responses to visual motion". In: Neurosci.
16.4 (Feb. 1996), pp. 1486-1510.
[3] Bruno B Averbeck, Peter E Latham, and Alexandre Pouget. "Neural correlations, population coding and computation". In: Nature reviews neuroscience 7.5 (2006), pp. 358-366.
[4] Marlene R Cohen and Adam Kohn. "Measuring and interpreting neuronal correlations". In:
Nat. Neurosci. 14.7 (June 2011), pp. 81 1-9-
[5] Stefano Panzeri, Monica Moroni, Houman Safaai, and Christopher D Harvey. "The structures and functions of correlations in neural population codes". : In: Nat. Rev. Neurosci. 23.9 (Sept. 2022). pp. 551-567.
[6] LF Abbott and P Dayan. "The effect of correlated variability on the accuracy of a population code". In: Neural Comput. 11.1 (Jan. 1999), pp. 91-101.
[7] Rubén Moreno-Bote, Jeffrey Beck, Ingmar Kanitscheider, Xaq Pitkow, Peter Latham, and Alexandre Pouget. "Information-limiting correlations". In: Nat. Neurosci. 17.10 (Oct. 2014), pp. 1410-1417.
[8] Dimitri Yatsenko, Krešimir Josié, Alexander S Ecker, Emmanouil Froudarakis, R James Cotton, and Andreas Tolias. "Improved estimation and interpretation of correlations neural circuits". In: PLOS Comput. Biol. 11.3 (Mar. 2015), e1004083.
[9] Alex H. Williams and Scott W. Linderman. "Statistical neuroscience in the single trial limit". In: Current Opinion in Neurobiology (2021). Computational Neuroscience, pp. 193-205. ISSN: 0959-4388.
[10] Oleg Rumyantsev, Jéróme Lecoq. Oscar Hernandez, Yanping Zhang, Joan Savall, Radostaw Chrapkiewicz, Jane Li, Hongkui Zeng, Surya Ganguli, and Mark Schnitzer. "Fundamental bounds the fidelity of sensory cortical coding". In: Nature 580.7801 (Apr. 2020), pp. 100- 105.
[11] lain M. Johnstone and Debashis Paul. "PCA in High Dimensions: An Orientation". In: Pro- ceedings of the IEEE 106.8 (2018), pp. 1277-1292.
[12] Joel Dapello, Jenelle Feather, Hang Le, Tiago Marques, David Cox, Josh MeDermott, James DiCarlo, and Sueyeon Chung. "Neural Population Geometry Reveals the Role of Stochasticity in Robust Perception". " In: Advances in Neural Information Processing Systems. Ed. by M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan. Vol. 34. Curran Associates, Inc., 2021, pp. 15595-15607.
[13] Lyndon R. Duong. Jingyang Zhou, Josue Nassar, Jules Berman, Jeroen Olieslagers, and Alex H. Williams. "Representational dissimilarity metric spaces for stochastic neural networks". In: International Conference on Learning Representationss 2023.
[14] Adam Kohn and Matthew A Smith. "Stimulus dependence of neuronal correlation in primary visual cortex the macaque". In: J. Neurosci. 25. 14 (Apr. 2005), pp. 3661-3673. [15] Marlene R. Cohen and William T. Newsome. "Context-Dependent Changes in Functional Circuitry in Visual Area MT In: Neuron 60.1 (2008), pp. 162-173. ISSN: 0896-6273.
[16] Adrian AanceeAlvarez, Alexander Thiele, Thomas D. Albright, Gene R. Stoner, and Gustavo Deco. "Stimulus-dependee" variability and noise correlations in cortical MT neurons". In: Proceedings of the National Academy of Sciences 1 (2013), pp. 13162-13167.
[17] C Gourieroux, Jasiak, and R Sufana. "The Wishart Autoregressive process of multivariate stochastic volatility". In: J. Econom. 150.2 (June 2009), pp. 167-181.
[18] Mohsen Pourahmadi. diig--dimensionl covariance estimation: with high-dimensional data.
Vol. 882. John Wiley & Sons, 2013.

---

[19] Olivier Ledoit and Michael Wolf. "A well-conditiones estimator for large- lrnge-dieensionll co- variance matrices". In: Journal of Multivariate Analysis 88.2 (2004). pp. 365-411, ISSN: 0047-259X.
[20] Mehdi Rahim, Bertrand Thirion, and Gaël Varoquaux. "Population shrinkage of covariance (PoSCE) for better individual brain -unceiinal-comeectiiityy estimation". In: Medical image analysis 54 (2019), pp. 138-148.
[21] Jerome Friedman, Trevor Hastic, and Robert Tibshirani. "Sparse inverse covariance estimation with the graphical lasso". In: Biostatistics 9.3 (Dec. 2007), pp. 432-441. ISSN: 1465-4644.
[22] Ming Bo Cai, Nicolas W Schuck, Jonathan W Pillow, and Yael Niv. "Representational structure task structure? Bias in neural representational similarity analysis and Bayesian method for reducing bias". In: PLoS Compul. Biol. 15.5 (May 2019), e1006299.
[23] Hsin-Hung Li, Thomas C. Sprague, Aspen H. Yoo, Wei Ji Ma, and Clayton E. Curtis. "Joint representation of working memory and uncertainty in human cortex". In: Neuron 109.22 (2021), 3699-3712.c6. ISSN: 0896-6273.
[24] Alexander Walther, Hamed Nili, Naveed Ejaz, Arjen Alink, Nikolaus Kriegeskorte, and Jöm Diedrichsen. "Reliability of dissimilarity measures for multi-voxel pattern analysis". In: Neuroimage 137 (Aug. 2016). pp. 188-200.
[25] Carsen Stringer, Michalis Michaelos, Dmitri Tsyboulski, Sarah E. Lindo, and Marius Pachitariu. "High-preeisee coding in visual cortex". In: Cell 184.10 (2021), 2767-2778.e11. ISSN: 0092- 8674.
[26] MohammadMehdi Kafashan, Anna W. Jaffe, Selmaan N. Chettih, Ramon Nogueira, Iñigo Romeeo, Christopher D. Harvey, Rubén Moreno-Bote, and Jan Drugowitsch. "Scaling of sensory information in large neural populations shows signatures of information-limiitng correlations". In: Nature Communications 12.1 (Jan. 2021), p. 473, ISSN: 2041-1723.
[27] M M Merzenich, P L Knight, and G L Roth. "Representation of cochlea within primary auditory cortex in the cat".. In: J. Neurophysiol. 38.2 (Mar. 1975), pp. 231-249.
[28) Elena Kudryavitskaya, Eran Marom, Haran Shani-Narkiss, David Pash, and Adi Mizrahi. "Flexible categorization in the mouse olfactory bulb". In: Curr. Biol. 31.8 (Apr. 2021), 1616- 1631.e4. .4..
[29] Nir Even-Chen Blue Sheffer, Saurabh Vyas, Stephen Ryu, and Krishna V Shenoy. "Structure and variability of delay activity in premotor cortex". In: PLoS Comput. Biol. 15.2 (Feb. 2019), e1006808.
[30] Omer Hazon, Victor H Minces, David P Tomas, Surya Ganguli, Mark Schnitzer, and Pablo Jercog. "Noise correlations in neural ensemble activity limit the accuracy of hippocampal spatial representations". In: Nat. Commun. 13.1 (July 2022), p. 4276.
[31] Adrien Peyrache, Marie M Lacroix, Peter C Petersen, and György Buzsáki, "Internally orga- nized mechanisms of the head direction sense". In: Nat. Neurosci. 18.4 (Apr. 2015), pp. 569- 575.
[32] Christopher KI Williams and Carl Edward Rasmussen, Gaussian processes for machine learning. Vol. 2. 3. MIT press Cambridge. MA, 2006.
[33] Byron M Yu, John P Cunningham, Gopal Santhanam, Stephen Ryu, Krishna V Shenoy, and Maneesh Sahani. "Gaussian-procee" factor analysis for low-dimensional single-trial analysis of neural population activity". In: Advances in Neural Information Processing Systems. Ed. by D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou. Vol. 21. Curran Associates, Inc., 2008
[34] Lea Duncker, Gergo Bohner, Julien Boussard, and Maneesh Sahani. "Learning interpretable continuous-time models of latent stochastic dynamical systems". In: Proceedings of the 36th International Conference on Machine Learning. Ed. by Kamalika Chaudhuri and Ruslan Salakhutdinov. Vol. 97. Proceedings of Machine Learning Research. Long Beach, California, USA: PMLR, 2019, pp. 1726-1734.
[35] Stephen Keeley, David Zoltowski, Yiyi Yu. Spencer Smith. and Jonathan Pillow. "Efficient non-conjugate Gaussian process factor models for spike count data using polynomial ap- proximations". In: International Conference on Machine Learning. PMLR. 2020, pp. 5177- 5186.
[36] Virginia Rutten, Alberto Bernacchia, Maneesh Sahani, and Guillaume Hennequin. "Non- reversible Gaussian processes for identifying latent dynamical structure in neural data". In: Advances in Neural Information Processing Systems. Ed. by H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin. Vol. 33. Curran Associates, Inc., 2020, pp. 9622-9632.
12

---

[37] Stephen Keeley, Mikio Aoi, Yiyi Yu, Spencer Smith, and Jonathan W Pillow. "Identifying signal and noise structure in neural population activity with Gaussian process factor models". In: Advances in Neural Information Processing Systems. Ed. by H. Larochelle, M. Ranzato, Hadsell, M.F. Balcan, and H. Lin. Vol. 33. Curran Associates, Inc., 2020, pp. 13795-13805.
[38] Anqi Wu, Nicholas A. Roy. Stephen Keeley, and Jonathan W Pillow. "Gaussian process based nonlinear latent structure discovery in multivariate spike train data". In: Advances in Neural Information Processing Systems. Ed. by I. Guyon, U. Von Luxburg. S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett. Vol. 30. Curran Associates, Inc., 2017.
[39] David Inouye, Eunho Yang, Genevera Allen, and Pradeep Ravikumar. "A Review of Multivari- ate Distributions for Count Data Derived from the Poisson Distribution". In: Wiley Interdiscip. Rev. Comput. Stat. 9.3 (Mar. 2017).
[40] György Buzsáki and Kenji Mizuseki. "The log-dynamic brain: how skewed distributions affect network operations". In: Nat. Rev. Neurosci. 15.4 (Apr. 2014), pp. 264-278.
[41] Feng Zhu, Harrison A Grier, Raghay Tandon, Changjia Cai, Anjali Agarwal, Andrea Gio- vannucci, Matthew T Kaufman, and Chethan Pandarinath. "A deep learning framework for inference of single-trial neural population dynamics from calcium imaging with subframe temporal resolution". In: Nat. Neurosci. 25.12 (Dec. 2022), pp. 1724-1734.
[42] Andrew Gordon Wilson and Zoubin Ghahramani. "Generalised Wishart processes". In: Pro- ceedings of the 27th Conference on Uncertainty in Artificial Intelligence, UAI 2011. AUAI Press. 2011, pp. 736-744.
[43] Andrew Gelman, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin. Bayesian data analysis (Third Edition). CRC press, 2013.
[44] Bernhard Scholkopf and Alexander J Smola. Learning with kernels. Adaptive Computation and Machine Learning series. London, England: MIT Press, Dec. 2001.
[45] Thomas Hofmann, Bernhard Schölkopf. and Alexander Smola, "Kernel methods in machine
learning". In: The Annals of Statisttics 36.3 (2008), pp. 1171-1220.
[46] David IC MacKay et al. "Introduction to Gaussian processes" ". In: NATO ASI series F computer and systems sciences 168 (1998), pp. 133-166.
[47] Neil C Rabinowitz, Robbe L. Goris, Marlene Cohen, and Eero P Simoncelli. "Attention stabilizes the shared gain of V4 populations". In: eLife 4 (Nov, 2015). Ed. by Matteo Carandini, e08998. ISSN: 2050-084X.
[48] David M Blei, Alp Kucukelbir, and Jon D McAuliffe. "Variational Inference: A Review for Statisticians". In: J. Am. Stat. Assoc. 112.518 (Apr. 2017), pp. 859-877.
[49] Creighton Heaukulani and Mark van der Wilk. "Scalable Bayesian dynamic covariance model- ing with variational Wishart and inverse Wishart processes". In: Advances in Neural Informa- tion Processing Systems. Ed. by H. Wallach, H. Larochelle, A. Beygelzimer, F. d' A'Alcéé-Bur -Bce E. Fox, and R. Garnett. Vol. 32. Curran Associates, Inc., 2019.
[50] Du Phan, Neeraj Pradhan, and Martin Jankowiak. "Composable Effects for Flexible and Accelerated Probabilistic Programming in NumPyro". In: arXiv preprint arXiv: 1554 (2019).
[51] Kevin Murphy. Machine learning: a probabilistic perspective. MIT press, 2012.
[52] Marino Pagan, Eero P Simoncelli, and Nicole C Rust. "Neural quadratic discriminant anal- ysis: Nonlinear decoding with V1-like computation" In: Neural computation 28.11 (2016), pp. 2291-2319.
[53] David A Rosenbaum. "Human movement initiation: Specification of arm, direction, and extent". In: J. Exp. Psychol. Gen. 109.4 (Dec. 1980), pp. 444-474.
[54] S P Wise. "The primate premotor cortex: past, present, and preparatory". In: Annu. Rev. Neurosci. 8 (1985), pp. 1-19.
[55] Gregorio Luis Galiñanes, Claudia Bonardi, and Daniel Huber. "Directional Reaching for Water as a Sortex-Dependent Behavioral Framework for Mice". In: Cell Rep. 22.10 (Mar. 2018), pp. 2767-2783
[56] Tanji and E V Evarts. "Anticipatory activity of motor cortex neurons in relation to direction of an intended movement". In: J. Newrophysiol. 39.5 (Sept. 1976), pp. 1062-1068.
[57] Messier and J F Kalaska. "Covariation of primate dorsal premotor cell activity with direction and amplitude during memorized-delay reaching task". In: J. Neurophysiol. 84.1 (July 2000),
pp. 152-165.
13

---

[58] Afsheen Afshar, Gopal Santhanam, Byron M. Yu, Stephen I. Ryu, Maneesh Sahani, and Krishna V. Shenoy. "Single- Trial Neural Correlates of Arm Movement Preparation". In: ISSN: 0896-6273.
Neuron 71.3 (2011), pp. 555-564.
[59]
Gustav Theodor Fechner. Elemente der psychophysik. Vol. 2. Breitkopf Härtel, 1860.
[60]
Michael A. Webster, Karen K. De Valois, and Fugene Switkes. "Orientation and spatial- frequency discrimination for luminance and chromatic gratings". In: J. Opt. Soc. Am. A 7.6 (June 1990), pp. 1034-1049.
[61] D Burr and $ A Wijesundra. "Orientation discrimination depends on spatial frequency". In: Vision Res. 31.7-8 (1991), pp. 1449-1452.
[62] R Vogels and G A Orban. "How well do response changes of striate neurons signal differences in orientation: study in the discriminating monkey". In: Neurosci. 10.11 1 (Nov. 1990), pp. 3543-3558.
[63] Olivier Hénaff, Yoon Bai, Julie Charlton, Ian Nauhaus, Eero Simoncelli, and Robbe Goris. "Primary visual cortex straightens natural video trajectories". In: Nat. Commun. 12.1 (Oct. 2021), p. 55982.
[64]
Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov. "Importance weighted autoencoders".
In: 2015.
14


Reviews:
[{"review_id": "kjGcyjqYBk", "reviewer": "Reviewer_sHno", "summary": "- Exhibits good performance, especially for recordings of large neural populations with very few trials per condition.\n- Enables dense sampling of conditions with only a single trial per condition, unlike standard estimators that require a large number of trials per condition (although based on empirical assumptions).", "questions": "- Figure 2: What does the \"Wishart oracle\" mentioned in the legend represent? Although it is mentioned, it does not appear to be shown in any panel. Additionally, why is \"Empirical\" only plotted in panels D3 and D4? It is assumed that its log probability in panel C diverges, but what about D1 and D2? Furthermore, why does the mean operator norm not align with the log probability for C3 and D3? In D4, one would expect the gap between \"Grand Empirical,\" assuming infinite smoothness, and Wishart to consistently widen with lower smoothness.\n- Figure 4: How is it possible for condition 1 to be present in both the training conditions (panel C) and the test conditions (panel D)?\n- What is the computational runtime compared to standard covariance estimators, particularly in comparison to the \"Grand Empirical\" estimator that performs well on the real experimental data shown in Figures 3 and 4?", "limitations": "The authors have discussed limitations of their work.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Exhibits good performance, especially for recordings of large neural populations with very few trials per condition.\n- Enables dense sampling of conditions with only a single trial per condition, unlike standard estimators that require a large number of trials per condition (albeit grand empirical estimator).", "weaknesses": "- Computationally more expensive than standard covariance estimators", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> **Computationally more expensive than standard covariance estimators**\n\nWe have clarified that the computational burden is reasonable: it takes 80 seconds to fit a dataset with 100 neurons, 80 conditions, 32 trials per condition. Notice that we run 10,000 iterations of our optimization algorithm, therefore the run time for such a dataset is about 8 milliseconds per iteration. The runtime for other algorithms for a similarly sized dataset are the following (in seconds): PoSCE: 45, Graphical Lasso: 24, Ledoit-Wolf: 0.3, and Empirical: 0.3.\n\nExperimental data in neuroscience is extremely expensive and time-consuming to collect. Equipment and facilities alone cost tens of thousands of dollars, and animals may need to be trained for weeks (if not months) to perform certain tasks. It is also ethically imperative to use as few animal subjects as possible (particularly for non-human primates). In short, the experimental datasets we target are extremely valuable: neuroscientists are highly motivated to perform high-quality analyses, even if they require marginally more computational resources.\n\nWe also wish to remind the reviewer that our model has fundamentally new capabilities that are not present in simple baselines. In particular, we can infer a continuous manifold for the mean and covariance of neural responses (see Fisher Information analysis in Rebuttal Fig 1D), generalize to entirely unseen conditions (see Fig. 3D), and quantify our uncertainty in a Bayesian framework. A fair comparison between inference times should take this difference into account.\n\n> **Figure 2: What does the \"Wishart oracle\" mentioned in the legend represent? Although it is mentioned, it does not appear to be shown in any panel.**\n\nThank you for noticing. This corresponds to an experiment that we removed before submission but forgot to remove from the legend. Briefly, this corresponded to a Wishart process but with parameters initialized to the ground truth. Our intention was for this to be a sanity check – an “upper bound” on our model’s performance. However, we intend to remove it in the final version for a streamlined presentation.\n\n> **Additionally, why is \"Empirical\" only plotted in panels D3 and D4? It is assumed that its log probability in panel C diverges, but what about D1 and D2?**\n\nThis has been an oversight, we apologize for this. We have included the revised figures in Rebuttal Fig 1A,B. However, empirical covariance matrices were ill-conditioned and did not produce proper log-likelihood estimates to be included in Fig. 2 C1-4.\n\n> **Furthermore, why does the mean operator norm not align with the log probability for C3 and D3?**\n\nThey are measuring performance in two different ways, so it is possible for a model to outperform another in log probability but vice versa on the covariance operator norm. Generally, these two measures of performance are qualitatively aligned (except for the portions of C3 and D3) that the reviewer notes. We will add a short explanation to section 3.1 in the revision.\n\n> **In D4, one would expect the gap between \"Grand Empirical,\" assuming infinite smoothness, and Wishart to consistently widen with lower smoothness.**\n\nWe agree. This is generally what we observe in D4. The reviewer may be noting there is an exception at very low degrees of smoothness. Here, all models perform badly since there is no good way to share information across conditions and there are too few trials per condition. It is curious that Grand Empirical “beats” the WP here, which may be due to the variational inference getting caught in sub-optimal local minima. Ultimately, we think this edge case is not a regime of primary interest, since no model performs well (i.e. covariance estimation is not tractable).\n\n> **Figure 4: How is it possible for condition 1 to be present in both the training conditions (panel C) and the test conditions (panel D)?**\n\nThere are two separate analyses. In panel C, we train on all conditions and test on a subset of trials per condition. In panel D, we train on a subset of conditions and test on held-out conditions. We show condition 1 in each panel since it was held out in the second analysis. We will update the legend titles and figure legend to make this more clear.\n\n> **What is the computational runtime compared to standard covariance estimators, particularly in comparison to the \"Grand Empirical\" estimator that performs well on the real experimental data shown in Figures 3 and 4?**\n\nPlease see our response above. In short, we believe the computational burden is negligible. Our revision will include a table showing the computational scaling of our algorithm as a function of the number of neurons, conditions, and trials per condition. Roughly speaking, we expect computation times to be on the order of minutes for large-scale neural data.\n\n"}, {"review_id": "AihukkIJlY", "reviewer": "Reviewer_mr49", "summary": "The estimation of noise covariance in neuroscience is limited by the experimental difficulty of obtaining large numbers of trials for the same neurons, but also the desirability of large numbers of neurons. Here the authors begin by recognizing this need for lower-variance estimators of the noise covariance, and especially those that leverage correct assumptions specific to this domain. This submission centers around leveraging the fact that noise covariance is often estimated not only by repeating a single stimulus many times, but also by showing very similar stimuli parameterized by a continuous number (such as grating orientation), which should evoke similar responses. To leverage this assumption, the authors build a probabilistic model of neural responses in which assumes that the response statistics vary smoothly with a known parameter x. Specifically, transformed neural responses are modeled as a Gaussian distribution whose mean and factorized covariance ($U$) are Gaussian Processes of x (with a squared exponential kernel). The parameters of this probabilistic model, including the covariance as a function of x, are inferred via mean-field variational inference. The estimation method is validated on simulated data and neural recordings in mouse and macaque.", "questions": "I have a request for some analyses that would show how well this model captures neural responses.\n - First, for the presented neural data it would be nice to see well this model compares in log likelihood to a probabilistic baseline that is more expressive than a GP mean. The ELBO of a VAE comes to mind, but there are many options in a good-faith effort to build powerful (if uninterpretable) probabilistic models.\n - It would also be nice to see an analysis of synthetic data where some of the assumptions are broken. For example, what happens on data in which the smoothness parameter is not uniform but is itself a function of the inputs? This seems like a near-certain outcome in real neural responses. (For example, in efficient codes of natural images, certain small changes in inputs result in larger changes in neural activity. As per the oblique effect, one might expect that neural responses vary more quickly w/r/t orientation for the more-frequent cardinal orientations as compared to less-frequent diagonal orientations). A discussion of these figures could help ensure this estimator is used where appropriate.\n\nFor clarity, could the graphical model of this generative model and variational inference procedure be presented as a figure?\n\nNot much was said about the intricacies of learning and inference. Figures in this direction would likely help adopters troubleshoot the problems usually found in ELBO-trained methods. E.g. what is the tradeoff of accuracy vs speed in the number of MC samples used to estimate of the gradient of the loss? \n - Related: In Fig 3B, is this the max, min, and mean (across seeds) of the median LL across 30 folds? If so, that would narrow these distributions considerably relative to just the distribution of LL across folds on a single run. This would somewhat obscure the variance of the estimator. A thorough analysis of the variance of this estimator would be interesting, help instill confidence in its use, and align with the paper's main narrative (as it is likely smaller than naïve covariance estimation).\n\nMight the kernel be learned via a more expressive function? It strikes me that contrastive learning methods learn such kernels (and are often talked about as learned distance metrics) and could be readily slotted into this framework without (in my opinion) much loss in interpretability. A transfer-learning/ foundation model approach would mitigate the large sample size requirement of a large model like this. There was a sentence about this in the Discussion, so count this as a vote of excitement in that future direction. \n", "limitations": "Some limitations were discussed. Restating some things above, it’d be nice to see more acknowledgements of the limitations relating to the implied GP encoding model and the potential shortcomings of the nonconvex, ELBO-based inference method.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "This well-written paper tackles a practical problem already faced in many experimental designs in neuroscience, and will likely find wide use if the code is well-documented. The solution is elegant, offers a nice balance of simplicity and power, and benefits from its specific tailoring of known statistical methods to a particular statistical problem in the domain of neuroscience.", "weaknesses": "Although I see the benefits of model simplicity, one wonders what this model might miss about the neural response. It is worth noting that the estimated covariance is, strictly speaking, the estimated covariance of the model's residual error of its mean estimate. If this mean is poorly estimated, the noise covariance will be larger than actuality, if not wholly inaccurate. The accuracy of the mean response, in turn, depends on the assumptions built into the model. How much is it accurate to say that neural responses, in general, have a mean firing rate that is a GP of the stimulus parameter with a squared-exponential kernel? The utility of analyzing the resulting $\\Sigma(x)$ depends on how much one trusts that this implied encoding model is a good one.\n\nWith this said, it is at least more transparent here that analyzing summary statistics like noise correlations always implies, under the hood, a statistical model of neural responses. It is nice here that this is explicit, and I think a step in the right direction.\n\nFinally, it would help the significance somewhat if this estimation method could help discover (or refine) a scientific conclusion in, say, the Allen dataset. What more can we learn from existing data with a lower-variance covariance estimate?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> **Although I see the benefits of model simplicity [...].**\n\nWe agree that accurately inferring the mean is an important aspect of accurately inferring the covariance. Using a GP to enforce some smoothness in the mean response across conditions is indeed reasonable. See, for example: Wu et al. (\"Gaussian process based nonlinear latent structure discovery in multivariate spike train data.\", NeurIPS 2017).\n\nThe details of the model will likely depend on the particular neural dataset under study. It is certainly possible that other kernels are better suited to certain situations. Our code package will of course have the standard kernels implemented so that users can try different forms of the model.\n\nWhile we think this is a reasonable first pass, it is easy to swap in new kernel functions, methods for kernel learning, or even methods for non-stationary Gaussian process inference as discussed below. We are excited to try these new directions in future work.\n\n> **[...] it is at least more transparent here [...].**\n\nThank you. We agree with these sentiments.\n\n> **[...] discover (or refine) a scientific conclusion [...].** \n\nThank you for this comment, which motivated us to flesh out how WPs can provide estimates of Fisher Information that are (a) more accurate, and (b) come with Bayesian uncertainty estimates. Fisher Information is a fundamental quantity in theoretical neuroscience that is used to characterize thresholds of perception (see Moreno-Bote et al., 2014, Nature Neuroscience). Please see the general rebuttal and Rebuttal Fig D for more information.\n\n> **[...] probabilistic baselines.**\n\nThank you for the interesting suggestion. In our initial discussions, we were not able to find a straightforward way of performing this experiment. Popular VAEs consider a factorization in the latent space which is independent across different dimensions, estimating a diagonal noise covariance in the latent space as opposed to the high-dimensional firing rate space.\n\nDeep VAEs are data-hungry (require many trials) while the applications we are interested in are very trial-limited. We also feel that doing a comprehensive comparison to VAEs is outside the scope of this paper since VAEs represent a very broad class of models with different objective functions (e.g. beta-VAE), architectures, and hyperparameters. Please let us know if you have a concrete model in mind that could be implemented in the time frame of the review process. We are happy to try.\n\n> **[...] some of the assumptions are broken [...].**\n\nWe thank the reviewer for the insightful suggestion. Indeed model misspecification is likely to occur for any model of neural data. We would like to re-emphasize that a GP or WP prior does not necessarily mean that the neighboring conditions are equi-distanced in the mean or covariance space. The prior encourages the nearby conditions to have close means and covariances, once the model is conditioned on the observed data, it allows for the distances to be influenced by the empirical values. In technical terms, GPs with stationary kernels are asymptotically consistent (see the introduction of Koepernik et al. \"Consistency of Gaussian process regression in metric spaces.\" for a review). \n\nRegarding the suggestion, again our initial discussion was not conclusive about how to perform this experiment. GPs with stationary kernels do not allow for a varying smoothness parameter. One way to incorporate the reviewer’s suggestion would be to use non-stationary kernels, but choosing a non-stationary kernel that’s consistent with neuroscientific observations is not straightforward. We believe our work opens up the possibility of adapting non-stationary Gaussian process methods to this problem, and we now discuss this as a future direction in our revision. See Paun et al. \"Stochastic variational inference for scalable non-stationary Gaussian process regression.\" for more information.\n\n> **For clarity, could the graphical model [...]?**\n\nPlease see Rebuttal Fig 1G,H. We will update the manuscript with this information.\n\n> **[...] intricacies of learning and inference [...].**\n\nIn our code package, all the parameters for every run of the algorithm are stored in configuration (yaml) files allowing for inspecting and modifying the parameters. An example config file is included in the uploaded code package attached to the supplementary.\n\nFortunately, we have not observed sensitivities to the optimization parameters. We’ve used Adam optimizer with a learning rate of 0.001 with a single MC sample and 10000 iterations for all experiments. In our experiments, using multiple MC samples did not lead to qualitatively different outcomes. As expected, the run time is moderately slower with multiple MC samples; overall run time is still on the order of a few minutes, so this is not a major concern. The number of MC samples is a tunable parameter in our code package.\n\n> **Related: In Fig 3B, is this the max, min [...].**\n\nWe are not 100% certain what the reviewer is requesting. We reported the median log-likelihood of individual observations in the test set, which is similar to the mean. We plan to substitute the mean log-likelihood in our revision since it produced essentially identical results (see Rebuttal Fig E,F).\n\nThe variance of the LL across individual samples is not related to the variance of the estimator. The variance of an estimator is a frequentist notion of uncertainty that does not neatly map onto the Bayesian approach that we’ve adopted here. We don’t think that the variance of LL across individual samples is usually of interest to report, but perhaps we are missing something or misinterpreting the reviewer’s request. We are happy to continue the discussion.\n\n> **[...] more expressive function [...].**\n\nWe thank the reviewer for sharing our excitement. Indeed this is a direction that we are considering and excited about for future work. "}, {"review_id": "hqTi0Z0a2a", "reviewer": "Reviewer_yodS", "summary": "This work proposes to use Wishart Processes, originally proposed by Wilson and Ghahramani, to estimate the covariance structure of neural activity in experiments where there are parametric variations linking stimulus conditions. By pooling estimates appropriately across task conditions, the limited number of experimental trials can be used to provide more reliable estimates of variability. Inference is performed via a mean field variational approach, and the model is applied to synthetic data and two neural data sets.\n\nThis seems like a reasonable approach to data of the type the authors are modeling, and the experiments show that it outperforms existing methods in the literature. Conversely, this is a very straightforward application of an existing method, which may be novel to neuroscientists, but does not represent much of new conceptual contribution. Moreover, the results presented, while likely better estimates of covariance, do not produce any neuroscientific insights reported here.", "questions": "- ll. 26-29: The number of parameters grows quadratically, but if all pairs of neurons are simultaneously observed, it doesn't follow that the process of estimating correlations does not scale well. Rather, the argument should be from the limited number of trials in each condition, right?\n- What is the difference between a Wishart process and a GP with a kernel that includes both trial and neuron indices? Some discussion of this might be helpful for those unfamiliar with the method.\n- Details of the variational approximation in the supplement are somewhat sparse. Do I understand correctly that the authors are performing full GP inference on $\\mathbf{U}$ and $\\boldsymbol{\\mu}$? What exactly is the mean field form of $q_\\phi$? Does it factorize over $(\\mathbf{U}, \\boldsymbol{\\mu})$ or something else? Finally, the reference on line 7 of Supplement B is not included there, and I don't think it's the same as reference 3 in the main text.", "limitations": "- This work is focused on neuroscience tasks in which trials are sampled from a set of parametrically related experimental conditions. This is true of many kinds of experiments but will not apply to many other cases that do not share that structure.\n- The authors discuss (in Supplement B) several methods for mitigating the poor scaling behavior of Gaussian Processes with $N$, $P$, and $C$ but do not implement any of the well-known scalable GP methods that would mitigate this, since their examples are for only moderate $N$.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- Application of Bayesian nonparametric techniques for estimating covariance to a new domain.\n- Pools strength across parametrically related stimulus conditions, making good use of a limited number of trials.\n- Clear presentation.", "weaknesses": "- This is a relatively straightforward application of an existing method. \n- Inferential details of the model are somewhat sparsely described. (See questions below.)\n- Results for the experiments are evaluated on various performance metrics, but it is unclear what new insights this method provides. It is perfectly reasonable to employ a better estimation method for this model, but the experiments in Figures 3 and 4 do not produce any additional findings beyond evidence that the new method works better. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> ***This is a relatively straightforward application of an existing method***\n\nThis is one of the most important points to discuss so we have laid out a detailed response in our “general rebuttal.\" To summarize our arguments:\n* “Application” papers of all stripes are requested in the “NeurIPS Call for Papers”\n* Our submission fits into a specific tradition of method papers in the “Neuroscience and Cognitive Science” track.\n* We extended the model to incorporate low-rank plus diagonal structure and we have now further extended the model to handle Poisson-distributed noise.\n* Even “vanilla” Wishart Processes are relatively exotic models that (to our knowledge) have not been widely applied outside of quantitative finance.\n\n> ***The number of parameters grows quadratically, but if all pairs of neurons are simultaneously observed, it doesn't follow that the process of estimating correlations does not scale well***\n\nOur original statement&mdash;\"number of estimated parameters grows quadratically... while the number of measurements grows only linearly\"&mdash;is correct.\n\nE.g. Imagine a case where `N=2` neurons and `K=150` trials. Thus, we have `N * K = 300` measurements of firing rate. We want to estimate a `2 x 2` symmetric covariance, which has 3 parameters. Thus, the ratio of measurements to parameters is 100 to 1.\n\nNow imagine if we have `N=100` neurons, and we want to estimate `N * (N + 1) / 2` parameters in the covariance matrix. Intuitively, if we want the ratio of measurements to parameters to stay the same as before (i.e. 100 to 1), we would need to chose `K` to satisfy `N * K = 100 * N * (N + 1) / 2`. That is, we would need `K = 100 * 201 / 2 = 10050` trials.\n\nIn summary, the difficulty of the problem depends *both* on `N` and `K`. We will edit our paper to convey this more clearly. Note that the intuition we provided above can be made fully rigorous. E.g. Vershynin (2011) “How Close is the Sample Covariance Matrix to the Actual Covariance Matrix?”\n\n> ***What is the difference between a Wishart process and a GP with a kernel that includes both trial and neuron indices?***\n\nAs a refresher, our method uses a GP to estimate the mean response and a Wishart process (WP) to estimate the covariance. Both the GP and WP have a kernel that measures similarity *across different experimental conditions*. The hyperparameters of this kernel are related to how smooth the neural response is as a function of changing the condition. In our model there is no kernel that measures similarity across neurons or trials. We consider these extensions in turn below:\n\n*GP model with kernel over neurons.* This would incorporate a smoothness prior over the neuron indices (e.g. neuron 1’s response would be highly correlated with neuron 2’s response, and de-correlated with neuron 53’s response). In most experiments, neurons are labeled arbitrarily so it typically does not make sense to add this structure.\n\n*GP model with kernel over trial.* This would incorporate a smoothness prior over trial index, which could be helpful to model slow drift or non-stationarity in the neural response. Since multiple conditions are often randomly interleaved across trials, it may be better to use the absolute time between two trials (rather than their integer index) to quantify this non-stationarity.\n\nNote that adding a time component to the GP kernel  *would only model non-stationarity in the mean* neural response.* It would *not* capture correlation structure across neurons (which is our primary motivation). Thus, an interesting extension of our model would be to add a time component to *both the GP and WP kernels*.\n\nWe welcome further discussion from the reviewer, in case we have misinterpreted their suggestion.\n\n> ***it is unclear what new insights this method provides***\n\nWe agree that we can do more to spell out the scientific insights our model helps enable. We would like to point the reviewer to some existing results in the paper, as well as some new experiments that we’ve performed in the rebuttal phase.\n\n*Results in the paper:* (a) smooth covariance interpolation in Fig. 3D  (b) Generating covariances using a single trial per condition in Supplemental Fig. 6 (c) drawing a full posterior distribution over means and covariances attaching uncertainties to covariance estimates.\n\n*New results:* (a) assessing Fisher Information in Rebuttal Fig 1D (b) modeling covariances for Poisson distributed observations in Rebuttal Fig 1E (c)  linear vs. quadratic discriminant analysis in Rebuttal Fig 1C.\n\nFinally, the introduction section of our paper cites multiple neuroscience papers where covariance plays an important quantitative role. Our model is broadly applicable, and so opens up the possibility of many new scientific insights across these cited works.\n\nWe hope that the reviewer can understand it is challenging to validate a new form of analysis while simultaneously providing a scientific breakthrough. The primary focus of our paper is methodological, but we have done our best to outline scientific future directions (see our \"general rebuttal\").\n\n> ***Details of the variational approximation in the supplement are somewhat sparse.***\n\nWe regret that we were not clear and understand the importance of fixing these details. By “full GP inference” we interpret the reviewer to mean the closed form solution to GP regression with Gaussian observation noise. We do not use this anywhere since the WP posterior does not admit a closed form solution. Thus, to approximate the posterior we perform joint variational inference over the GP and WP parameters (they do not factorize as pointed out by Reviewer mr49).  Each parameter is modeled with a Gaussian with a learnable mean and variance, and no correlations are modeled across parameters. That is, we perform standard mean field variational inference (e.g. Blei et al. 2017, “Variational Inference: A Review for Statisticians”).\n\nWe are happy to answer more queries in the discussion period.\n\n\n\n"}, {"review_id": "P5xETQ4p0H", "reviewer": "Reviewer_ruJu", "summary": "The goal of this paper is to compute the covariance of recorded neurons in a given stimulus condition with a low number of samples. Although the conditions are different, some aspects are shared which justifies the fact that the covariance for a specific condition should depend on the covariance in other conditions. This idea is implemented using a Bayesian model which imposes smoothness of the means and models the covariance as $L (U U^\\top + \\Lambda) L^\\top$ where $L$ does not depend on the condition and $U$ and $\\Lambda$ evolve smoothly with the condition.", "questions": "- Can you comment on the absence of a method mixing population and condition-specific covariance ?\n- Is there a way to design a predictive task based on your data instead of relying on visual inspection?\n- Can you comment on the computation time of your method ?\n- How is the cross-validation made? Is there an internal cross-validation to select the hyper-parameters?", "limitations": "N.A", "rating": 7, "confidence": 4, "soundness": 2, "presentation": 4, "contribution": 3, "strengths": "Overall, the paper is very well-written and addresses a significant problem in neuroscience.\nThe use of smoothness in the means and covariance estimates across conditions is new as far as I know. The authors validate their approach by experimenting with synthetic and real data.\n", "weaknesses": "While the smoothness part is new, exploiting similar conditions to yield a better covariance estimate is not. The example I have in mind is [1] where condition-specific covariances are biased towards the population covariance but since this idea is very natural, I assume more examples exist. \n\n[1] Rahim, Mehdi, Bertrand Thirion, and Gaël Varoquaux. \"Population shrinkage of covariance (PoSCE) for better individual brain functional-connectivity estimation.\" Medical image analysis 54 (2019): 138-148.\n\nIn particular, I find the experiments a bit biased toward the success of the method presented by the authors. The absence of any competitor that uses both empirical and population covariance estimates makes the comparison a bit unfair (obviously models that do not do this will not perform well by design of the generative mechanism). \n\nThe other major drawback is the lack of clear metrics of performance on real data that demonstrate predictive power (e.g. a regression or classification task). The current metrics on the real data experiment rely mostly on visual inspection. As a side note, the visualization of principal components may hide some of the signal especially if the generated data is non-linear. Predictive tasks are better in that they give a clearer metric of performance. Is there a way to predict something related to the subjects (their age, their performance for instance) so that we have a proxy for the quality of the estimated covariance?\n\nA potential weakness (although I am not sure for this one), is the lack of internal cross-validation. Unless I am wrong, I believe the hyper-parameters are optimized on the same cross-validation loop as the one used for testing. If I am wrong, please clarify in the text by clearly stating that there are two nested cross-validation loops (an internal CV on the train set and an outer one to select the train and test splits).\n\nLastly, compared to Ledoit-Wolf or others, I would assume that the model presented by the author is much slower to fit. If this is the case, this should be clearly stated (a figure with the running time in function of N would be great).\n\nThere are small parts in the submission that I did not understand:\n- How do you optimize L?\n- Is $\\bar{\\mu}$ in equation (4) the average over samples in all conditions ? I believe this is never defined. \n- When you compute the log-likelihood (for instance to produce Figure 3 B) you need a model of the data. Is this model the one described by (3) ?\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> **While the smoothness part is new, exploiting similar conditions [...].**\n\nThank you for the pointer to the PoSCE paper. We had overlooked it since it is a neighboring field (none of us have experience with fMRI data).\n\nRoughly speaking, PoSCE can be adapted to our problem by estimating the covariance in condition $c$ as a weighted combination of the empirical/sample covariance in condition $c$ and the grand covariance across conditions. Our original work already compared to the two extremes – sample covariance and grand covariance were included as baselines – but it is a good idea to include a weighted combination of these two estimates as a stronger baseline.\n\nIn actuality, PoSCE posits a probabilistic model over the tangent space of the PSD manifold and uses a geometric mean to perform the averaging. However, when we applied their code to our data it spits out NaNs. Given time constraints, and since conceptually simpler baselines may be preferable to a reader anyways, we decided to take a weighted linear combination of the per-condition sample covariance and the grand covariance. We feel this is equally justified as the Riemannian distance used by PoSCE, since sample covariances are already the arithmetic mean of the second moments: $\\Sigma = (1/n)\\sum_i \\mathbf{x_i x_i^\\top}$.\n\nWe find that the Wishart process outperforms this modified PoSCE baseline on real and synthetic data (see “Weighted Avg” results in Rebuttal Fig F).\n\nImportantly, even if the PoSCE baseline was competitive in terms of performance, we still think our work would stand as an important contribution. In particular, PoSCE is not able to capture smoothness across conditions (as the reviewer mentions). It therefore cannot extrapolate or interpolate between conditions, which is one of our central motivations. See, e.g., the application to Fisher Information in the general rebuttal.\n\n> **The other major drawback is the lack of clear metrics [...].**\n\nThank you for this suggestion. Inspired by this comment, we’ve performed a new analysis to demonstrate how improved covariance estimation may lead to improved decoding of the experimental condition. See Rebuttal Fig C and our general response to all reviewers for more details about this experiment.\n\n> **A potential weakness [...] is the lack of internal cross-validation [...].**\n\nOur original submission used an initial cross-validation run to show how sensitive the log-likelihood was to hyperparameters such as smoothness and number of components in Fig. 4B. We then fixed these hyperparameters and ran them over 30 additional train-test folds. Since each of these 30 folds was a fully randomized split, we think it is unlikely that we would have overfit our small number of hyperparameters.\n\nHowever, we agree with the reviewer that the most rigorous way to do this is to have randomized cross-validation train-validation-test folds, and to select new hyperparameters for each fold on the basis of the validation set. We have now done this on the monkey data and show our results in Rebuttal Fig F. The results are effectively unchanged, but we will swap this new analysis into the final revision because we agree that it is more rigorous.\n\nWe will also clearly explain the cross-validation process in our methods – in particular, for each train-validation-test fold we do a randomized hyperparameter search over the kernel smoothness parameters and select the best model on the validation set for testing. We welcome any additional requests for details.\n\n> **Lastly, [...] the model presented by the author is much slower to fit [...].**\n\nPlease check our response to R1 for the wall clock time comparisons. Although the computational complexity of the model is tractable (as mentioned in supplementary B.1.1) our model takes longer than compared methods to run. Notice that this is due to a multitude of reasons. First, instead of a single estimate of the covariance, we infer a whole distribution over means and covariances per condition which can be used to compute uncertainties associated with inferred means and covariances and also to sample from them. In addition, since our model is GP, WP based, we uncover an underlying continuous space of means and covariances which can be used for interpolating test conditions and computing properties such as curvature and gradients wrt conditions. Developing full-blown Bayesian statistical models like these usually require more computational resources as opposed to ad-hoc point-estimate or frequentist counterparts. Although our implementation is relatively fast, we recognize that there’s still a lot of room for improving the run time of the algorithm using the latest developments in the GP literature. We leave this for future work.\n\n> **How do you optimize L?**\n\nThe matrix $L$ is initialized according to the Cholesky factorization of the grand empirical covariance and it’s optimized to maximize the evidence lower bound (ELBO). Note that $L$ is not a latent random variable: it is a parameter that we optimize as a point estimate. In contrast, we optimize a distribution over the latent variables of the model (e.g. covariance factors, denoted $U$ in the paper). These are jointly optimized with parameters like $L$. We regret that our description was too brief and caused confusion; we will improve the revision to make the paper more self-contained.\n\n> **Is $\\bar{\\mu}$ in equation (4) [...].**\n\nWe apologize for this oversight. $\\bar{\\mu}$ is the mean function for the prior GP, in all of our experiments we set it to the constant function zero, we will include both of these missing details in the revised manuscript.\n\n> **When you compute the log-likelihood [...].**\n\nYes. In all of our original analyses F(.) is a multivariate normal distribution as described on line 104. Since our original submission, we have also extended the model to be compatible with Poisson noise (see Rebuttal Fig E).\n\nAll the questions in the **Questions** section are addressed above."}]

Author rebuttal:
# General Response
We thank the reviewers for their thoughtful and productive critiques which did not identify any major technical errors. We have done our best to incorporate all reviewer feedback and requests for more details (see individual responses to each reviewer).

Below we summarize the two most important considerations that we hope the reviewers and AC take into account during their final deliberation.

**Technical contribution.** The reviewers largely agreed that our paper addresses an important open problem, is well-written, and is germane to the NeurIPS audience. One potential exception is Reviewer yodS, who writes “this is a very straightforward application of an existing method, which may be novel to neuroscientists, but does not represent much of new conceptual contribution.” We would like to highlight the following:

- The NeurIPS call for papers invites papers on “Applications.” It is especially common for statistical neuroscientists to publish methodological papers in the “Neuroscience and Cognitive Science” conference track. By adapting an existing method to neural data, we believe our paper fits squarely within this tradition.

- Even outside of neuroscience, Wishart processes are a relatively exotic model. Prior applications of the method appear mostly confined to quantitative finance.

- We were not able to identify an implementation of this model that could be used off-the-shelf on our data. Indeed, the Wilson and Ghahramani paper cited by the reviewer used a different inference method (MCMC). We found it necessary to use a more recent variational inference procedure that appeared at NeurIPS in 2019 by Heaukulani & van der Wilk. We found it necessary to implement their model from scratch. A deliverable outcome of our paper will be a self-contained code package that is customized to handle neural data (see supplement for a lightweight version of this package).

- Most importantly, we found that an extension of the model with low-rank-plus-diagonal structure outperformed the full-rank Wishart process. We have not seen such an extension in prior work, even though it has multiple benefits: the learnable diagonal term ensures that matrices are well-conditioned during inference, has fewer parameters (discourages overfitting), and is a reasonable prior for a “spiked covariance” model. If the reviewers agree this is novel, we will revise the text of our paper to emphasize this technical innovation.

- Since our original submission, we have also obtained good results using a Wishart process with Poisson-distributed observations (see Rebuttal Fig E). This is an important consideration for neural data, as evidenced by many prior submissions to the neuroscience track at NeurIPS (e.g. Zhou & Wei 2020 “pi-VAE”). Again, we are unaware of any paper that describes this variation of a Wishart process.

In summary, we respectfully disagree with the notion that our work is merely a “straightforward application.” It involved implementing the model from scratch, adapting the model (e.g. with low-rank-plus-diagonal structure), and benchmarking the model against multiple baselines and across two datasets spanning different species and experimental paradigms. Publishing this effort is necessary to spark an interest in this class of under-appreciated models. Indeed, every colleague in computational neuroscience that we’ve asked is unaware of Wishart process models.

**Scientific contribution.** We agree with Reviewer mr49’s suggestion that “it would help the significance somewhat if this estimation method could help discover (or refine) a scientific conclusion.” Along similar lines, Reviewer ruJu inquired whether it is possible to “demonstrate predictive power (e.g. a regression or classification task).”

We have thus far pursued two scientific applications with promising results.

First, by inferring condition-specific covariance matrices, our method enables decoding conditions by quadratic discriminant analysis (QDA). Linear discriminant analysis (LDA) is a more common model in neuroscience, but quadratic decoders are biologically plausible and of interest to the field (Pagan et al., 2018, Neural Computation). In Rebuttal Fig C, we show that the QDA models enabled by our method outperform linear decoders. This directly answers ruJu’s query.

Second, an even more exciting direction is to use Wishart processes to quantify information-limiting noise correlations (Moreno-Bote et al. 2014, Nature Neurosci). Our method enables continuous interpolation of the mean and covariance across conditions which is necessary to compute the Fisher information. Furthermore, when a Gaussian process is differentiable, its derivative is also a Gaussian process. Thus, we can compute a posterior to quantify uncertainty in our estimate of Fisher information. We showcase our method in Rebuttal Fig D.

Together, the two results above provide a more concrete glimpse into how our methods can be used for scientific discovery. We hope the reviewers are sympathetic to the challenges of providing a complete scientific story in the same paper as an important methodological advance.

**Gameplan.** We will update Figure 2 with the missing lines pointed out by reviewer sHno (see Rebuttal Fig A-B). We will add Rebuttal Figs C and D to Figure 2 of the paper to demonstrate scientific applications and significance. We will describe the Poisson model in a supplemental section (see Rebuttal Fig E). We will swap Rebuttal Fig F in for Figure 3B to address a comment by reviewer ruJu. We will add the graphical model schematic (Rebuttal Fig G) to Figure 1. We will edit the text of the paper as described in the individual reviewer responses. We hope that the reviewers and AC will agree that these edits will substantially improve the manuscript without fundamentally changing the main story and results.


Meta-review: The reviewers were unanimous in their appraisal of this paper's contributions as above the bar for acceptance to NeurIPS, and I congratulate the authors on their detailed and comprehensive rebuttals, which were critical in leading some reviewers to raise their scores.  I'm pleased to report that this paper has been accepted to NIPS.  Congratulations!  Please revise the manuscript to address all reviewer comments and questions.
Decision: Accept (poster)