Scalable Bayesian GPFA with automatic relevance determination and discrete noise models
Authors: Kristopher T Jensen
Ta-Chu Kao
Jasmine Talia Stone
Guillaume Hennequin
Venue: NeurIPS 2021 Poster (2021)
Keywords: neuroscience
GPFA
latent variable models
motor control
neural dynamics
neural data analysis
variational inference

Abstract: Latent variable models are ubiquitous in the exploratory analysis of neural population recordings, where they allow researchers to summarize the activity of large populations of neurons in lower dimensional ‘latent’ spaces. Existing methods can generally be categorized into (i) Bayesian methods that facilitate flexible incorporation of prior knowledge and uncertainty estimation, but which typically do not scale to large datasets; and (ii) highly parameterized methods without explicit priors that scale better but often struggle in the low-data regime. Here, we bridge this gap by developing a fully Bayesian yet scalable version of Gaussian process factor analysis (bGPFA), which models neural data as arising from a set of inferred latent processes with a prior that encourages smoothness over time. Additionally, bGPFA uses automatic relevance determination to infer the dimensionality of neural activity directly from the training data during optimization. To enable the analysis of continuous recordings without trial structure, we introduce a novel variational inference strategy that scales near-linearly in time and also allows for non-Gaussian noise models appropriate for electrophysiological recordings. We apply bGPFA to continuous recordings spanning 30 minutes with over 14 million data points from primate motor and somatosensory cortices during a self-paced reaching task. We show that neural activity progresses from an initial state at target onset to a reach- specific preparatory state well before movement onset. The distance between these initial and preparatory latent states is predictive of reaction times across reaches, suggesting that such preparatory dynamics have behavioral relevance despite the lack of externally imposed delay periods. Additionally, bGPFA discovers latent processes that evolve over slow timescales on the order of several seconds and contain complementary information about reaction time. These timescales are longer than those revealed by methods which focus on individual movement epochs and may reflect fluctuations in e.g. task engagement.

Scalable Bayesian GPFA with automatic relevance determination and discrete noise models
Kristopher T. Jensen"
Ta-Chu Kao" Jasmine Stone Guillaume Hennequin
Department of Engineering
University of Cambridge (ktj21, tck29, jts58, gjeh2)
Abstract
Latent variable models are ubiquitous in the exploratory analysis of neural popula- tion recordings, where they allow researchers to summarize the activity of large populations of neurons in lower dimensional 'latent' spaces. Existing methods can generally be categorized into (i) Bayesian methods that facilitate flexible incorpora- tion of prior knowledge and uncertainty estimation, but which typically do not scale large datasets; and (ii) highly parameterized methods without explicit priors that scale better but often struggle in the low-data regime. Here, we bridge this gap by developing fully Bayesian yet scalable version of Gaussian process factor analysis (bGPFA), which models neural data as arising from a set of inferred latent processes with a prior that encourages smoothness over time. Additionally, bGPFA uses automatic relevance determination to infer the dimensionality of neural ac- tivity directly from the training data during optimization. To enable the analysis of continuous recordings without trial structure, we introduce novel variational inference strategy that scales near-linearly in time and also allows for non-Gaussian noise models appropriate for diectrophysiological recordings. We apply bGPFA to continuous recordings spanning 30 minutes with over million data points from primate motor and somatosensory cortices during a self self-paced reaching task. We show that neural activity progresses from an initial state at target onset to reach- specific preparatory state well before movement onset. The distance between these initial and preparatory latent states is predictive of reaction times across reaches, suggesting that such preparatory dynamics have behavioral relevance despite the lack of externally imposed delay periods. Additionally, bGPFA discovers latent processes that evolve over slow timescales on the order of several seconds and contain complementary information about reaction time. These timescales are longer than those revealed by methods which focus individual movement epochs and may reflect fluctuations in e.g. task engagement.
1 Introduction
The adult human brain contains upwards of 100 billion neurons [3]. Yet many of our day-to-day behaviors such as navigation, motor control, and decision making can be described in much lower dimensional spaces. Accordingly, recent studies across range of cognitive and motor tasks have shown that neural population activity can often be accurately summarised by the dynamics of "latent state" evolving in a low-dimensional space [7, 8, 13, 38, 43]. Inferring and investigating these latent processes can therefore help us understand the underlying representations and computations implemented by the brain [20]. To this end, numerous latent variable models have been developed and used to analyze the activity of populations simultancously recorded neurons. These models range from simple linear projections such as PCA to sophisticated non-linear and temporally correlated models [10, 16, 22, 43, 50].
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

prior
posterior
latent state
latent state
neuron
time
time
activity (neuron i)
activity (neuron I)
signal noise
time
latent state
latent state
Figure 1: Bayesian GPFA schematic. Bayesian GPFA places a Gaussian Process prior over the latent states in each dimension function of time (p( X [t); top left) well as a linear prior over neural activity as function of each latent dimension (p( F X); bottom left). Together with a stochastic noise process p( ()). which can be discrete for electrophysiological recordings, this forms a generative model that gives rise to observations Y (middle). From the data and priors, bGPFA infers posterior latent states for each latent dimension (p(X) Y); top right) as well as a posterior predictive observation model for each neuron (p(Yeestl Xtest; Y); bottom right). When combined with automatic relevance determination, the model learns automatically discard superfluous latent dimensions by maximizing the log marginal likelihood of the data (right, black vs. blue).
A popular latent variable model for neural data analysis is Gaussian process factor analysis (GPFA), which has yielded insights into neural computations ranging from time tracking to movement preparation and execution [1, 48, 49, 51]. However, fitting GPFA comes with a computational complexity of ((T)) and a memory footprint (TT)) for T time bins. This prohibits the application of GPFA to time series longer than a few hundred time bins without artificially chunking such data into "pseudo-triale trialss and treating these as independent samples. Additionally, canonical GPFA assumes a Gaussian noise model while recent work has suggested that non non-Gaussian models often perform better on neural data [12, 25, 60]. Here, we address these challenges by formulating a scalable and fully Bayesian version of GPFA (bGPFA; Figure 1) with a computational complexity of O log T') and memory cost of OTT). To do this, we introduce an efficiently parameterized variational inference strategy that ensures scalability to long recordings while also supporting non- Gaussin noise models. Additionally, the Bayesian formulation provides framework for principled model selection based on approximate marginal likelihoods [53]. This allows us perform automatic relevance determination and thus fit a single model without prior assumptions about the underlying dimensionality, which is instead inferred from the data itself [5, 41].
We validate our method on synthetic and biological data, where bGPFA exhibits superior performance to GPFA and Poisson GPFA with increased scalability and without requiring cross-valions voss-validation to select the latent dimensionality. We then apply bGPFA to longitudinal, multi-area recordings from primary motor (MI) and sensory (S1) areas during a monkey self-paced reaching task spanning 30 minutes. bGPFA readily scales to such datasets, and the inferred latent trajectories improve decoding of kinematic variables compared to the raw data. This decoding improves further when taking into account the temporal offset between motor planning encoded by MI and feedback encoded by S1. We also show that the latent trajectories for M1 converge consistent regions of state space for a given reach direction at the onset of each individual reach. Importantly, the distance in latent space to this preparatory state from the state target onset predictive of reaction times across reaches, similar to previous results in a task that includes an explicit 'motor preparation epoch' where the subject is not allowed to move [1]. This illustrates the functional relevance of such preparatory activity and suggests that motor preparation takes place even when the task lacks well-defined trial structure and externally imposed delay periods, consistent with findings by Lara al. [32] and Zimnik and Churchland [62]. Finally, we analyze the task relevance of slow latent processes identified by bGPFA which evolve on timescales of seconds; longer than the millisecond timescales that can be resolved by methods designed for trial-structured data. We find that some of these slow processes are also predictive of reaction time across reaches, and we hypothesize that they reflect task engagement which varies over the course of several reaches.
2

---

2 Method
In the following, we use the notation A refer to the matrix with elements dij. We use a to refer to the kth row or column of A with an index running from 1 to K, represented column vector.
2.1 Generative model
Latent variable models for neural recordings typically model the neural activity RNNT of N neurons at times t E R as arising from shared fluctuations in D latent variables X E DRT Specifically, the probability of a given recording can be written as
p (
(1)
where F E NRT are intermediate, neuron-specific variables that can often be thought of firing rates or similar notion of noise-free activity. For example, GPFA [59] specifies
F( = II (unti Sni, x
(2)
12,
X 6(F CX)
(3)
p(x = II (a O, K.))
with K₂ = ka(t, t))
(4)
That is, the prior over the thh latent function (t) is Gaussian process [45] with covariance function ka(:s) (usually a radial basis function), and the observation model p(Y (X) is given by a parametric linear transformation with independent Gaussian noise.
In this work, additionally introduce prior distribution the mixing matrix C E RNND with hyperparameters specific to each latent dimension. This allows us to learn an appropriate latent dimensionality for a given dataset using automatic relevance determination (ARD) similar previous work in Bayesian PCA (Appendix I; 5) rather than relying on cross-validation or ad-hoc thresholds of variance explained. Unlike in standard GPFA, the log marginal likelihood (Equation 1) becomes intractable with this prior. We therefore develop novel variational inference strategy [55] which also (i) provides a scalable implementation appropriate for long continuous neural recordings, and (ii) extends the model to general non- Gaussian likelihoods better suited for discrete spike counts.
In this new framework, which we call Bayesian GPFA (bGPFA), C' of the form Gnd ~ (0, 2)). where 8a a scale parameter associated with latent dimension d. Integrating C out in Equation then yields the following observation model:
Gaussian prior over
p(F O,
with S diag(s si.gg
(5)
Moreover, we use a general noise model p(Y = II P(umelfre) where p(ynelfrs) is any distribution for which we can evaluate its density.
2.2 Variational inference and learning
To train the model and infer both X and F from the data Y, we use nested variational approach. It intractable to compute log p( Y (t) (Equation 1) analytically for bGPFA, and we therefore introduce lower bound on log p( Y tt at the outer level and another These lower bounds are constructed from approximations to the posterior distributions over latents (X) and noise-free activity (F) respectively.
log p Y (Y|X) at the inner level.
Distribution over latents At the outer level, we introduce a variational distribution q( (X) over la- tents and construct an evidence lower bound (ELBO; 55) on the log marginal likelihood of Equation 1:
log (Ylt) > L == E [log [)] [gX) (xX())))))
Conveniently, maximizing this lower bound is equivalent to minimizing KL (x(X)))) and thus also yields an approximation to the posterior over latents the form of q( X). We estimate the first term of the ELBO using Monte Carlo samples from (XX) and compute the KL term analytically.
3

---

Here, we use so-called whitened parameterization of ((X) [19] that is both expressive and scalable large datasets;
((x) II Ma with Nd Mu= d=1
Kiva and Ea Bj=KAAATKK
(7)
where K2 / is any square root of the prior covariance matrix Ka. and va R is vector of variational parameters to be optimized. Ar RTxT is positive semi-definite variational matrix whose structure chosen carefully so that its squared Frobenius norm, log determinant, and matrix-vector products can all be computed efficiently which facilitates the evaluation of Equations 8 and 9. This whitened parameterization has several advantages. First, it does not place probability mass where the prior itself does not. In addition to stabilizing learning [39], this also guarantees that the posterior is temporally smooth for a smooth prior. Second, the KL term in Equation 6 simplifies to
1 Σ (1A 2 log Aal + || -T)-
XL) ( [0]
(8)
2
d
Third, ((X can be sampled efficiently via a differentiable transform (i.e. the reparameterization trick) provided that fast differentiable K3 v and ary products are available for any vector w:
dd (im) K3 (vz Aana) with ni I),
where and ~ q(xa)) This is important to form a Monte Carlo estimate of By(x) [log Y X)].
To avoid the challenging computation of KJ u for general K₂ [2], we directly parameterize K the positive definite square root of K, which implicitly defines the prior covariance function ka(;-)) In
this work we an RBF kernel for Ka and give the expression for K² in Appendix F. Additionally, T)
use Toeplitz acceleration methods to compute K j products in OT log time and with ((T) memory cost [48, 57].
We implement and compare different choices of Aa in Appendix F. For the experiments in this work, we use the parameterization A, = Li' aaCa. where a is diagonal with positive entries and Ca is circulant, symmetric, and positive definite. This parameterization enables cheap computation of KL divergences and matrix-vector products while maintaining sufficient expressiveness (Appendix F). All results are qualitatively similar when instead using simple diagonal parameterization Aa V.
Distribution over neural activity Evaluating log p( Y|x) En, p(y X) for each sample
drawn from ((X) is intractable for general noise models. Thus, we further lower-bound the ELBO of Equation by introducing an approximation qffn| X)) to the posterior p(fn) yn. X):
loy X) 2 (x((,) [logr(y> [x)] - KL [g(fn) X |p( f X)]
(10)
We repeat the whitened variational strategy described at the outer level by writing
x(fn|X)= = =N(fn: p., Σn,) with f.. k à and Σ, K L, (K.))"" where in E R is a neuron- meeron-specific vector of variational parameters to be optimized along with a power-triangular matrix L,, E Dx and K denotes the covariance matrix of (X), whose square root Ki = x S follows from Equation 5. The low-rank structure of K enables cheap matrix-vector products and KL divergences:
(11)
(p) [X] X)] = (I 2 [L,] 1221l=-D1
(12)
Note that the KL divergence does not depend on X in this whitened parameterization (Appendix H). Moreover, x in Equation 11 has the form of the exact posterior when the noise model is Gaussian (Appendix G), and it is equivalent to a stochastic variational inducing point approximation [18) for general noise models (Appendix H).
Finally, we need to compute the first term in Equation 10:
(xf [log p(y, ((y.lff)) = Σ ExGtally [log plynlllll
(13)
Fach term in this sum is simply 1-dimensional Gaussian expectation which can be computed analytically in the case of Gaussian or Poisson noise (with an exponential link function), and otherwise approximated efficiently using Gauss-Hermin quadrature (Appendix K; 18).

---

2.3 Summary of the algorithm
Putting Section 2.1 and Section 2.2 together, optimization proceeds at each iteration by drawing M
Monte Carlo samples {Xm H1 from ((X) and estimating the overall ELBO as:
1 M
Σ
Σ
Ll (ogg(ssssss) m,t
X Xmmaq(X
Σ KL
Σ KL KL (((#a))lo(aa))
(14)
where the expectation over qnf n X) is evaluated analytically or using Heruss- quadrature depending on the noise model (Appendix K). We maximize & using stochastic gradient ascent with Adam [29]. This has a total computational time complexity of (MMTT + M DT T) and memory complexity of (M where N is the number of neurons. T'' the number of time points, and D the latent dimensionality. For large datasets such as the monkey reaching data in Section 3.2, we compute gradients using mini-batches across time to mitigate the memory cost. That is, gradients for the sum over in Equation 14 are computed in multiple passes. The algorithm is described in pseudocode with further implementation and computational details Appendix L. The model learned by bGPFA can subsequently be used for predictions on held-out data by conditioning on partial observations as used for cross-validation in Section 3. and discussed in Appendix M. Latent dimensions that have been 'discarded' by automatic relevance determination will automatically have negligible contributions to the resulting posterior predictive distribution since the prior scale parameters sa are approximately zero for these dimensions (sce Appendix I for details).
3 Experiments and results
3.1 Synthetic data
We first generated an example dataset from the GPFA generative model (Equations 2-4) with true latent dimensionality 3. We proceeded to fit both factor analysis (FA), GPFA, and bGPFA with different latent dimensionalities D E [1, 10]. Here, we fitted bGPFA without automatic relevance determination such that sa = s Vd. As expected, the marginal likelihoods increased monotonically with D for both FA and GPFA (Figure 2a; Appendix 1). In contrast, the bGPFA ELBO reached its optimum value at the true latent dimensionality D = 3. This is a manifestation of "Occam's 's razor", whereby fully Bayesian approaches favor the simplest model that adequately explains the data Y [36]. When instead considering the cross-validated predictive performance of each method, performance deteriorated rapidly for D > 3 for FA and GPFA, while Bayesian GPFA was more robust to overfitting (Figure 2b). Notably, the introduction of ARD parameters {sa} in bGPFA allowed us to fit a single model with large D 10. This recovered the maximum ELBO of bGPFA without ARD and the minimum test eтror across GPFA and bGPFA without ARD (Figure 2a and b. black) without priori assumptions about the latent dimensionality or the need to perform extensive cross-validation. Consistent with the ground truth generative process, only 3 of the scale parameters sa remained well above zero after training (Figure 2b, inset). Similar this illustrative example with Gaussian data, bGPFA with ARD and Poisson noise also exceeded the optimal performance of Poisson GPFA when applied to both synthetic and experimental spike count data (Appendix E).
We then proceeded to apply bGPFA (D- = 10) to an example dataset drawn using Equations 4 and 5 with ground truth dimensionality D 2, and either Gaussian, Poisson, or negative binomial noise. For all three datasets, the learned parameters clustered into a group of two latent dimensions with high information content (Appendix J) and group of eight uninformative dimensions, consistent with the generative process (Figure 2c). In each case, we extracted the inferred latent trajectories corresponding to the informative dimensions and found that they recapitulated the ground truth up to linear transformation (Figure 2d). Fitting flexible noise models such as the negative binomial model important because neural firing patterns are known to be overdispersed in many contexts [4, 15, 54]. However, it is often unclear how much of that overdispersion should be attributed to common fluctuations in hidden latent variables (X our model) compared to private noise processes in single neurons [35]. In our synthetic data with negative binomial noise, we could accurately recover the single-neuron overdispersion parameters (Figure 2e; Appendix K), suggesting that such unsupervised models have the capacity to resolve overdispersion due to private and shared processes.
5

---

DGPFA (no ARD) DEPPA (ARD)
b
GPFA
-1.4
2.4
LUELBO
0.0
-6
lo
-1.6
1.6
-1.7
1.2
2
6
10
2
6
10
latent dimensionality
latent dimensionality
C
d
ground truth
Goussian
Poisson
latent dim 2
[lwllf
0.2
NegBinom
Kor N
5.0
-2.5
latent dim
log Sa
log
Figure 2: Bayesian GPFA applied to synthetic data. (a) Log likelihoods of factor analysis (yellow) & GPFA (green) and ELBO of Bayesian GPFA without ARD (blue) fitted to synthetic data with ground truth dimensionality of three for different model dimensionalities. bGPFA with ARD recovered a diree-dimensional latent space as well as the optimum ELBO of bGPFA without ARD (black dashed line). (b) cross-validated prediction errors for the models in (a) (Appendix M). bGPFA with ARD recovered the performance of the optimal GPFA and bGPFA models without requiring search over latent dimensionalities. Inspection of the learned prior scales {sa) and posterior mean parameters ||al|| (inset) indicates that ARD retained only D' right) and discarded the other 7 dimensions (bottom left). Shadings in (a) and (b) indicate +2 stdev. across model fits. (c) Learned parameters bGPFA with ARD and either Gaussian, Poisson or negative binomial noise models fitted to two-dimensional synthetic datasets with observations drawn from the corresponding noise models (Appendix K). The parameters clustered into two groups of informative (top right) and non-informative (bottom left) dimensions (Appendix J). (d) Latent trajectory in the space of the two most informative dimensions (c.f. (c)) for each model with the ground truth shown in black. (e) The overdispersion parameter Kn, for each neuron learned in the negative binomial model, plotted against the ground truth (Appendix K). Solid line indicates y = æ; corresponds to Poisson noise model.
informative dimensions (top
note that 6.
In summary, bGPFA provides flexible method for inferring both latent dimensionalities, latent trajectories, and heterogeneous single-neuron parameters in an unsupervised manner. In the next section, show that the scalability of the model and its interpretable parameters also facilitate the analysis of large neural population recordings.
3.2 Primate recordings
In this section, apply bGPFA to biological data recorded from a rhesus macaque during self- paced reaching task with continuous recordings spanning minutes (37, 42; Figure 3a). The continuous nature of these recordings as one long trial makes it challenging dataset for existing analysis methods that explicitly require availability many trials per experimental condition [43], and poses computational challenges to Gaussian process-based methods that cannot handle long time series [59]. While the ad-hoc division continuous recordings into surrogate trials can still enable the use of these methods [28]. here show that formulation of bGPFA readily applies to long continuous recordings. We fitted bGPFA with a negative binomial noise model to recordings from both primary motor cortex (M1) and primary somatosensory cortex (S1). For all analyses, a single recording session (indy_20 2ei_200104266 as in 28), excluded neurons with overall firing rates below 2 Hz, and binned data at 25 resolution. This resulted in data array Y E IR x 70482 (130 MI neurons and 70 S1 neurons).
used
We first fitted bGPFA independently to the M1 and S1 sub-populations with D = 25 latent dimensions. In this case, ARD retained 16 (MI) and (S1) dimensions (Figure 3b). We then proceeded to train a linear decoder to predict hand kinematics in the form of and y hand velocities from either

---

the inferred firing rates or the raw data convolved with 50 ms Gaussian kernel (28; Appendix M). We found that the model learned by bGPFA predicted kinematies better than the convolved spike trains, suggesting that (i) the latent space accurately captures kinematic representations, and (ii) the denoising and data-sharing across time in bGPFA aids decodability beyond simple smoothing of neural activity. Interestingly, by repeating this decoding analysis with an artificially imposed delay between neural activity and decoded behavior, we found that neurons in SI predominantly encoded current behavior while neurons in MI encoded a motor plan that predicted kinematics 100-150 ms into the future (Figure 3b). This is consistent with the motor neuroscience literature suggesting that MI functions as dynamical system driving behavior via downstream effectors [8].
We then fitted bGPFA to the entire dataset including both MI and S1 neurons. In this case, bGPFA retained 19 dimensions (Appendix D), and kinematic predictions improved over individual M1- and S1-based predictions (Figure 3b). In this analysis, the decoding performance as a function of delay between neural activity and behavior exhibited a broader peak than for the single-region decoding. We hypothesized that this broad peak reflects the fact that these neural populations encode both current behavior in SI as well as fitture behavior in M1 (Figure 3c). Indeed, when we took this offset into account by shifting all MI spike times by + 100 ms and retraining the model, decoding performance increased from 68.56% + 0.09 to 69. 89. + 0.06 (mean + sem variance explained across ten model fits; Appendix M). Additionally, the shifted data exhibited a narrower decoding peak attained for near-zero delay between kinematics and latent trajectories (Figure 3d). Consistent with the improved kinematic decoding, we also found that shifting the MI spikes by 100 ms increased the ELBO per neuron (-34,637.0 = 0.7 to 34,631.1 1.1 + 0.6) and decreased the dimensionality of the data (Appendix D; 46). These results suggest that M1 and S1 contain both overlapping but also non-redundant information, and that the most parsimonious description of the neural data is recovered by taking into account the different biological properties of M1 and $1.
We next wondered if bGPFA could be used to reveal putative motor preparation processes, which is non-trivial due to the lack of trial structure and well-defined preparatory epochs. We partitioned the data post-hoc into individual 'reaches', each consisting of a period of time where the target location remained constant. For these analyses, we only considered 'successful' reaches, where the monkey eventually moved to the target location, and we defined movement onset as the first time during reach where the cursor speed exceeded a low threshold (Appendix A). We began by visualizing the latent processes inferred by bGPFA as they unfolded prior to movement onset in each reach epoch. For visualization purposes, we ranked the latent dimensions based on their learned prior scales (a measure of variance explained; Appendix J) and selected the first two. Prior to movement onset, the latent trajectories tended to progress from their initial location at target onset towards reach-specific regions of state space (see example trials in Figure 3e for leftward and rightward reaches). To quantify this phenomenon, we computed pairwise similarities between latent states across all 762 reaches, during (i) stimulus onset and (ii) 75 ms before movement onset (chosen such that it is well before any detectable movement: Appendix A). We defined similarity as the negative Euclidean distance between latent states and restricted the analysis to 'fast' latent dimensions with timescales smaller than 200 ms to study this putatively fast process. When plotted as a function of reach direction, the latent similarities at target onset showed little discernable structure (Figure 3f, left). contrast, the pairwise similarities became strongly structured 75 ms before movement onset where neighboring reach directions were associated with similar preparatory latent states (Figure 3f, right). Similar albeit noisier results were found when using factor analysis or GPFA instead of bGPFA (Appendix A). These findings are consistent with previous reports of monkey MI partitioning preparatory and movement-related activity into distinct subspaces [14, 32), as well as with the analogous finding that a 'relative target' subspace is active before a 'movement subspace' in previous analyses of this particular dataset [28].
Previous work on delayed reaches has shown that monkeys start reaching earlier when the neural state attained at the time of the go cue which marks the end of delay period with a known reach direction - is close to an "opttmal subspace" [1, 23]. We wondered if a similar effect takes place during continuous, self-initiated reaching in the absence of explicit delay periods. Based on Figure 3e, we hypothesized that the monkey should start moving earlier if, at the time the next target is presented, its latent state is already close to the mean preparatory state for the required next movement direction. To test this, we extracted the mean preparatory state 75 ms prior to movement onset (as above) for each reach direction in the dataset. We found that the distance between the latent state at target onset and the corresponding mean preparatory state was strongly predictive of reaction time (Figure 3g,
7

---

a
task structure
b
51
varance captured
51
verlance captured
0.4
100ms shift
0.2
0.3
0.3
&
100ms shift
0.0
0.1
0.1
200
200
Sg
delay
delay
target anset
pre-movement
425
read timre (ma)
neach direction
labent dim
right left
125
direction
prep
Figure 3: Bayesian GPFA applied to primate data. (a) Schematic illustration of the self-paced reaching task. When target a 17x8 grid is reached (arrows; 8x8 shown for clarity), new target lights up on the screen (colours), selected at random from the remaining targets. In several analyses, we classify movements according to reach angle measured relative to horizontal (01. 62). (b) Learned mean and scale parameters for the bGPFA models. Small prior scales sy and posterior mean parameters (/]aall²) (2)) indicate uninformative dimensions (Appendix J). (c) We applied bGPFA to monkey M1 and S1 data during the task and trained a linear model to decode kinematics from firing rates predicted from the inferred latent trajectories with different delays between latent states and kinematics. Neural activity was most predictive future behavior in M1 (black) and current behavior in S1 (blue). Dashed lines indicate decoding from the data convolved with Gaussian filter. (d) Decoding from bGPFA applied to the combined MI and SI data (cyan). Performance improved further when decoding from latent trajectories inferred from data where MI activity was shifted by 100 ms relative to S1 activity (green). (e) Example trajectories in the two most informative latent dimensions for five rightward reaches (grey) and five leftward reaches (red). Trajectories are plotted from the appearance of the stimulus until movement onset (circles). During 'movement preparation', the latent trajectories move towards a consistent region of latent state space for each reach direction. (f) Similarity matrix of the latent state at stimulus onset showing no obvious structure (left) and 75 ms prior to movement onset showing modulation by reach direction (right). (g) Reaction time plotted against Euclidean distance between the latent state at target onset and the mean preparatory state for the corresponding reach direction (p- 0.45).
Pearson p = 0.45, = x 10 10-36). Such a correlation was also weakly present with factor analysis (p 0.21, 1.1 10-8) but not detectable in the raw data 0.002, p 0.95). We also verified that the strong correlation found with bGPFA was not an artifact of the temporal correlations introduced by the prior (Appendix C). Taken together, an important part of reaching movements even in an unconstrained self-paced task. Additionally, we showed that bGPFA captures such behaviorally relevant latent dynamics better than simpler alternatives, and our scalable implementation enables its use on the large continuous reaching dataset analysed here.
results suggest that motor preparation is
3.3 Long-timescalt latent processes
Some latent dimensions inferred by bGPFA also had long timescales on the order of 1.5 seconds, which is similar to the timescale of individual reaches (1-2 seconds; Appendix C). We hypothesized that these slow dynamics might reflect motivation or task engagement. Consistent with this hypothesis, we found that one of the slow latent processes (r 1.4 s) was strongly correlated with reaction 0.40, p 3.4 x 10-2s). Interestingly, the information
time during successful reaches (Pearson contained about reaction time this long timescale latent dimension was largely complementary to that encoded by the distance to preparatory states in the 'fast' dimensions (Appendix C), suggesting

---

0.8
time (minutes)
Speed (rui)
0.6
10-
latent state
15-
0.4
20-
0.2
0.0
20 25
10
1.5 (minutes)
20
25
time e (minutes)
time (minutes)
Figure 4: Analysis of period without task participation. (a) Cursor speed over the course of the recording session. Blue horizontal lines indicate the last successful trial before and first successful trial after period with no active task participation (blue shading). (b) Latent similarity matrix as a function of time during the task. The latent dynamics during task participation occur in a largely orthogonal subspace to the dynamics during the period with no active task participation. (c) Plot of latent state over time for a long- timescale latent dimension strongly correlated with reaction time (r 1.4 s).
that motor preparation and task engagement are orthogonal processes both contributing task performance.
The experimental recordings were also characterized by period of approximately five minutes towards the end of the recording session during which the monkey did not participate actively in the task and the cursor velocity was near-constant at zero (Figure 4a). When analysing neural activity across the periods with and without task participation, we found that neural dynamics moved to different subspace the monkey stopped engaging with the task (Figure 4b). Importantly, were able to simultaneously capture these context-depende and preparatory dynamics (Section 3.2) by fitting single model to the full 30 minute dataset. This suggests that bGPFA can capture behaviorally relevant dynamics within individual contexts even when trained richer datasets with changing contexts.
changes as well as movement-speciff
Finally, we wondered how the neural activity patterns during periods with and without task participa- tion were related to the long-timen timescaecale latent dimensions predictive of task engagement. Here found that the slow latent process considered above also exhibited a prominent change to a different state the monkey stopped participating in the task (Figure 4). This is consistent with our hypothesis that this latent process captures a feature related to task engagement which slowly deteriorated during the first 24 minutes of the task followed by a discrete switch to state with engagement in the task. During the period of active task participation, this latent dimension was also correlated with time within the session. Indeed, reach number and latent state were both predictive of reaction time, but with the latent trajectory exhibiting a slightly stronger correlation (Pearson p = 0.40 vs. = 0.37). It is surprising that task engagement decreases with time, and is in this case difficult to tease apart how motivation and time are differentially represented in such latent processes. However, based on the strong and abrupt modulation by task participation, this latent dimension appears represent an aspect of engagement with the task beyond the passing of time.
Taken together, thus find that bGPFA is capable of capturing not only single- dynamics and preparatory activity but also complementary processes evolving over longer timescales, which would be difficult to identify with methods designed for the analysis of many shorter trials.
4 Discussion
Related work The generative model of bGPFA can be considered an extension of the canonical GPFA model proposed by Yu al. [59] to include Gaussian prior over the loading matrix C (Section 2.1). In this view, bGPFA is to GPFA what Bayesian PCA is to PCA [5]; in particular, it facilitates automatic relevance determination to infer the dimensionality of the latent space from data [5, 41, Similar to previous work in the field, we also use variational inference to facilitate arbitrary observation noise models, including non-Gaussian models more appropriate for electrophy recordings [12, 25, 34, 50, 60, 61]. While variational inference has proven useful framework for such non-conjugate likelihood models, alternative approaches exist including the use of polynomial approximations to the non-linear terms in the likelihood [24]. Another major challenge in the

---

development of GP-based latent variable models such as bGPFA is to ensure scalability for longer time series. In this work, we utilize advances in variational inference [30, 47] to facilitate scalability to the large datasets recorded in modern neuroscience. In particular, we contribute a new circulant variational GP posterior expressed partly in the Fourier domain that both accurate and scalable. This is similar to Keeley et al. [25], who address the problem of scalability by assuming independence across Fourier features and formulating variational inference in the Fourier domain. However, we instead perform inference in the time domain and include additional factors in our variational posterior that ensure smoothness over time and allow for non-stationary posterior covariances. In contrast to these approaches, Zhao and Park [60] rely a low rank approximation to the prior covariance for inference and temporal subsamples for hyperparameter optimization to overcome the computational cost of model training. A conceptually similar approach employed by Duncker and Sahani [12] is the use of inducing points which has been studied extensively in the Gaussian process literature [17, 18, 52]. However, such low rank approximations can perform poorly on long time series where the number of inducing points needed is proportional to the recording duration [6].
bGPFA is also closely related to Gaussian process latent variable models (GPLVMs) [33, 53] which have recently found use in the neuroscience literature as a way of modelling flexible, nonlinear tuning curves [22, 34, 58]. This is because integrating out the loading matrix C' in p Y X) with Gaussian prior gives rise to a Gaussian process with a linear kernel. The low-rank structure of this linear kernel yields computationally cheap likelihoods, and our variational approach to estimating log p X is in fact equivalent to the sparse inducing point approximation used in the stochastic variational GP (SVGP) framework [17, 18]. In particular,
variational posterior is the same as that which would arise in SVGP with at least D inducing points irrespective of where those inducing points are placed (Appendix H). We also note that for Gaussian noise model, the resulting low-rank Gaussian posterior is the form of the exact posterior distribution (Appendix G). Additionally, since the bGPFA observation model and prior over latents are both GPs, bGPFA is an example of deep GP [1]] with two layers- the first with an RBF kernel and the second with a linear kernel. Finally, our parameterizations of the posteriors q (ara) and g(Jn) can be viewed as variants of the 'whitening' approach introduced by Hensman et al. [19] which both facilitates efficient computation of the KL terms in the ELBOs and also stabilizes training (Section 2.2).
Conclusion & impact In summary, bGPFA is an extension of the popular GPFA model in neu- roscience that allows for regularized, scalable inference and automatic determination of the latent dimensionality as well as the use of non non- Gaussian noise models appropriate for neural recordings. Importantly, the hyperparameters of bGPFA are efficiently optimized based on the ELBO on training data, which alleviates the need for cross- cross-valiattinn or complicated algorithms otherwise used for hyperparameter optimization in overparameterized models [16, 22, 27, 28, 58, 59]. Our approach can also be extended in several ways to make it more useful to the neuroscience community. For example, replacing the spike count-based noise models with a point process model would provide higher temporal resolution [12], and facilitate inference of optimal temporal delays across neural populations [31]. This will likely be useful as multi-region recordings become more prevalent in neuroscience [26]. Additionally, by substituting the linear kernel in Y (X) for an RBF kemnel in Euclidean space [58] or on a non-Euclicn Bucliucide manifold [22], we can recover scalable versions of recent GPLVM-based tools for neural data analyses with automatic relevance determination.
While such latent variable models have already played major role in systems neuroscience, progress in the field is still hampered by the relative inaccessability of these methods to the general neuroscience community. Here, we provide a ready-to-use Python package with GPU implementations of not only bGPFA with ARD, but also standard GPFA and Factor Analysis with Gaussian and non-Gaussian noise models. We hope that this implementation will help bridge the gap between methods development and the practical use of these methods by the community, similar to other recent efforts to democratize and compare latent variable models for neuroscience [44]. In future work, will also be interesting to extend bGPFA to non-linear observation models p(Y XX) and provide a general-purpont package for latent variable modelling in neuroscience. As progress in neuroscience continues to accelerate, these methods are likely to find extensive use in medical and industrial settings with applications e.g. to brain-computer interfaces in human subjects [21, 40, 56). With the rise of such technology transfer, the research community should strive maintain research transparency by providing open source code for methods development, and to remain aware of ethical issues associated with the increasingly narrow gap between research and applications in humans and other animals [9].
10

---

Acknowledgements
We are grateful to O'Doherty et al. [42] for making their data publicly available and Marine Schimel and David Liu for insightful discussions. We thank Marine Schimel, Yashar Ahmadian, Peter Stone, and Jonathan So for helpful comments on the manuscript. We thank David Liu for contributions to the codebase used for our analyses. K.T.J. was funded by a Gates Cambridge scholarship and J.T.S. by a Churchill scholarship.
References
[1] Afshar, A., Santhanam, G., Byron, M. Y., Ryu, S. L., Sahani, M., and Shenoy, K. V. (2011). Single-trial neural correlates of arm movement preparation. Neuron, 71(3):555-5644
[2] Allen, E., Baglama, J., and Boyd, S. (2000). Numerical approximation of the product of the square root of matrix with vector. Linear Algebra and its Applications, 310(1-3):167-118)
[3] Azevedo, F. A., Carvalho, L. R., Grinberg, L. T., Farfel, M., Ferretti, R. E., Leite, E., Filho, W. J., Lent, R., and Herculano-Houzell, S. (2009). Equal numbers of neuronal and nonneuronal cells make the human brain an isometrically scaled-up primate brain. Journal of Comparative Neurology, 513(5):532-544
[4] Azouz, R. and Gray, C. M. (1999). Cellular mechanisms contributing response variability of cortical neurons in vivo. Neurosci., 19(6):220--21) -223.
[5] Bishop, C. M. (1999). Bayesian PCA. Advances in neural information processing systems, pages 382-388.
[6] Chang, P. E., Wilkinson, W. J., Khan, M. E., and Solin, A. (2020). Fast variational learning in state-space Gaussian process models. In 2020 IEEE 30th International Workshop on Machine Learning for Signal Processing (MLSP), pages 1 1-6. IEEE.
[7] Chaudhuri, R., Gerçek, B., Pandey, B., Peyrache, A., and Fiete, (2019). The intrinsic attractor manifold and population dynamics of a canonical cognitive circuit across waking and sleep. Nature neuroscience, 2219):15221
[8] Churchland, M. M., Cunningham, J. P., Kaufman, M. T., Foster, J. D., Nuyujukian, P., Ryu, S. I., and Shenoy, K. V. (2012). Neural population dynamics during reaching. Nature, 487(7405):11-56
[9] Clausen, (2009). Man, machine and in between. Nature, 457(7233):0800-102
[10] Cunningham, J. P. and Byron, M. Y. (2014). Dimensionality reduction for large-scale neural recordings. Nature neuroscience, 17(1 10:1500-15
[11] Damianou, A. and Lawrence, N. D. (2013). Deep Gaussian processes. Artificial intelligence and statistics, pages 207-215. PMLR.
[12]
Duncker, and Sahani, M. (2018). Temporal alignment and latent Gaussian process factor inference in population spike trains. In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bae Assa-Bianchi, N., and Garnett, R., editors, Advances Neural Information Processing Systems, volume 31.
[13] Ecker, A., Berens, P., Cotton, R., Subramaniyan, M., Denfield, G., Cadwell, C., Smirnakis, S., Bethge, M., and Tolias, A. 2). State Dependence of Noise Correlations Macaque Primary Visual Cortex. Neuron, 12:235-248.
[14] Elsayed, G. F., Lara, A. H., Kaufman, M. T., Churchland, M. M., and Cunningham, J. P. (2016). Reorganization between preparatory and movement population responses in motor cortex. Nat. Commun., 7:1-15.
[15] Fenton, A. A. and Muller, R. U. (1998). Place cell discharge is extremely variable during individual passes of the rat through the firing field. Proceedings of the National Academy of Sciences, 318633188

---

[16] Gao, Y., Archer, E. W., Paninski, L., and Cunningham, J. P. (2016). Linear dynamical neural population models through nonlinear embeddings. In Lee, D., Sugiyama. M., Luxburg, U., Guyon, I,, and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 29.
[17] Hensman, J., Fusi, N., and Lawrence, N. D. (2013). Gaussian processes for big data. arXiv preprint 1rXiv:1309.68..
[18] Hensman, J., Matthews, A., and Ghahramani, Z. (2015a). Scalable variational Gaussian process classification. In Artificial Intelligence and Statistics, pages 351-360. PMLR.
[19] Hensman, J., Matthews, A. G., Filippone, M., and Ghahramani, Z. (2015b). MCMC for variationally sparse Gaussian processes. In Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 28.
[20] Humphries, M. D. (2020). Strong and weak principles of neural dimension reduction. arXiv preprint arXiv:2011.08088
[21] Jarosiewicz, B., Sarma, A. A., Bacher, D., Masse, N. Y., Simeral, J. D., Sorice, B., Oakley, M., Blabe, C., Pandarinath, C., Gilja, V., et al. (2015). Virtual typing by people with tetraplegia using a self-calibrating intracortical brain-computer interface. Science translational medicine, 7(313):31 33w1199-11110
[22] Jensen, K., Kao, T.-C., Tripodi, M., and Hennequin, G. (2020). Manifold GPLVMs for discovering non-euclidean latent structure in neural data. In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M. F,, and Lin, H., editors, Advances in Neural Information Processing Systems, volume 33, pages 225992
[23] Kao, T.-C., Sadabadi. M. S., and Hennequin, G. (2021). Optimal anticipatory control as a theory motor preparation: A thalamo-cortice circuit model. Neuron, 109:1567-1581.
[24] Keeley, S., Zoltowski, D., Yu, Y., Smith, S., and Pillow, J. (2020a). Efficient non-conjugate Gaus- sian process factor models for spike count data using polynomial approximations. In International Conference on Machine Learning, pages 5177-5186. PMLR.
[25] Keeley, S. L.,, Aoi, M. C., Yu, Y., Smith, S. L., and Pillow, W. (2020b). Identifying signal and noise structure in neural population activity with Gaussian process factor models. bioRxiv.
[26] Keeley, S. L., Zoltowski, D. M.. Aoi, M. C., and Pillow, J. W. (2020c). Modeling statistical dependencies multi-region spike train data. Current Opinion in Neurobiology.
[27] Keshtkaran, M. R. and Pandarinath, C. (2019). Enabling hyperparameter optimization in sequential autoencoders for spiking neural data. arXiv preprint arXiv: 1oXiv:1908.078%.
[28] Keshtkaran, M. R., Sedler, A. R., Chowdhury, R. H., Tandon, R., Basrai, D., Nguyen, S. L., Sohn, H., Jazayeri, M., Miller, L. E., and Pandarinath, C. (2021). A large-scale neural network training framework for generalized estimation of single-trial population dynamics. bioRxiv.
[29] Kingma, D. P. and Ba, J. (2015). Adam: A method for stochastic optimization. In Bengio, Y. and LeCun, Y., editors, 3rd International Conference on Learning Representations, ICLR 2015, San Diego CA, USA, May 7-9, 2015, Conference Track Proceedings.
[30] Kingma, D. P. and Welling, M. (2014). Auto-encoding variational bayes. In Bengio, Y. and LeCun, Y., editors, 2nd International Conference on Learning Representations, ICLR 2014, Banff AB, Canada, April 14-16, 2014, Conference Track Proceedings.
[31] Lakshmanan, K. C., Sadtler, P. T., Tyler-Kabara, E. C., Batista, A. P., and Yu, B. M. (2015). Extracting dineninnnll latent structure from time series in the presence of delays. Neural computation, 1825-11888886
[32] Lara, A. H., Elsayed, G. F., Zimnik, A. J., Cunningham, J. P., and Churchland, M. M. (2018). Conservation of preparatory neural events in monkey motor cortex regardless of how movement is initiated. eLife, 7:e31826.
[33] Lawrence, N. and Hyvärinen, A. (2005). Probabilistic non-linear principal component analysis with Gaussian process latent variable models. Journal of machine learning research, 6(11).
12

---

[34] Liu. D. and Lengyel, M. (2021). A universal probabilistic spike count model reveals ongoing modulation of neural variability. bioRxiv.
[35] Low, R. J., Lewallen, S., Aronov, D., Nevers, R., and Tank, D. W. (2018). Probing variability in a cognitive map using manifold inference from neural dynamics. BioRxiv, page 418939.
[36] MacKay, D. J. (2003). Information theory, inference and learning algorithms. Cambridge university press.
[37] Makin, J. G., O'Doherty, J. E., Cardoso, M. M., and Sabes, N. (2018). Superior amm- movement decoding from cortex with a new, unsupervised-leaminn algorithm. Journal of neural engineering, 15(2):026010.
[38] Minxha, J., Adolphs, R., Fusi, S., Mamelak, A. N., and Rutishauser, U. (2020). Flexible recruitment of memory-based choice representations by the human medial frontal cortex. Science, 368(6498).
[39] Murray, I. and Adams, (2010). Slice sampling covariance hyperparameters of latent Gaussian models. arXiv preprint arXiv: 0868.
[40] Musk, E. et al. (2019). An integrated brain-machine interface platform with thousands of channels. Journal of medical Internet research, 11610))11119.
[41] Neal, M. (2012). Bayesian learning for neural networks, volume 118. Springer Science Business Media.
[42] O' O'Doherty, J. E., Cardoso, M., Makin, J., and Sabes, P. (2017). Nonhuman primate reaching with multichannel sensorimotor cortex electrophysiology Zenodo http://doi. org/10.5281/zenodo, 583331.
[43] Pandarinath, C., O' Shea, J.J. J., Collins, J., Jozefowicz, R., Stavisky, S. D., Kao, J. C., Trautmann, M., Kaufman, M. T., Ryu, S. Hochberg, L. R., et al. (2018). Inferring single-trial neural population dynamics using sequential auto-encoders. Nature methods, 15(10):805-85 815.
[44] Pei, F., Ye, J., Zoltowski, D., Wu, A., Chowdhury, R. H., Sohn, H., O'Doherty, J. E., Shenoy, K. V., Kaufman, M. T., Churchland, M., Jazayeri, M., Miller, I. E., Pillow, J., Park, M., Dyer, E. L., and Pandarinath, C. (2021). Neural latents benchmark '21: Evaluating latent variable models neural population activity, arXiv preprint arXiv:2109.04463
[45] Rasmussen, C. E. and Williams, C. K. (1996). Gaussian processes for regression. MIT.
[46] Recanatesi, S., Ocker, G. K., Buice, M. A., and Shea-Brown, E. (2019). Dimensionality in recurrent spiking networks: global trends in activity and local origins in connectivity. PLoS computational biology, 1507)xe0006446.
[47] Rezende, J., Mohamed, S., and Wierstra, D. (2014). Stochastic backpropagation and approximate inference in deep generative models. In International conference on machine learning, pages 1278-1286. PMLR.
[48] Rutten, V., Bernacchia, A., Sahani, M., and Hennequin, G. (2020). Non-reversible Gaussian processes for identifying latent dynamical structure in neural data. Advances in Neural Information Processing Systems, 33.
[49] Sauerbrei, A., Guo, J.-Z., Cohen, J. D., Mischiati, M., Guo, W., Kabra, M., Verma, N., Mensh, B., Branson, K., and Hantman, A. W. (2020). Cortical pattern generation during dexterous movement input-driven. Nature, 5777799)
[50] Schimel, M., Kao, T.-C., Jensen, K. T., and Hennequin, G. (2021). iLQR- control-based learning of input-driven dynamics with applications to neural data. bioRxiv.
[51] Sohn, H., Narain, D., Meirhaeghe, N., and Jazayeri, M. (2019), Bayesian computation through cortical latent dynamics. Neuron, 103(5);3449947.
[52] Titsias, M. (2009). Variational learning of inducing variables in sparse Gaussian processes. In Artificial intelligence and statistics, pages 567-574. PMLR.
13

---

[53] Titsias, M. and Lawrence, N. D. (2010). Bayesian Gaussian process latent variable model. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statisties, pages 844-851. JMLR Workshop and Conference Proceedings.
[54] Tomko, G. J. and Crapper, D. R. (1974), Neuronal variability: non-stationary responses to identical visual stimuli. Brain research, 79(3):405 4188
[55] Wainwright, M. J. and Jordan, M. I. (2008). Graphical models, exponential families,, and variational inference. Now Publishers Inc.
[56] Willett, F. R., Avansino, D. T.,, Hochberg, L. R., Henderson, M., and Shenoy, K. V. (2021). pegh-ppeformance brain-to-text communication via handwriting, Nature, 593(7858):249-554
[57] Wilson, A. G., Dann, C., and Nickisch, H. (2015). Thoughts on massively scalable Gaussian processes. arXiv preprint arXiv: 151101187
[58] Wu, A., Roy, N. A., Keeley, S., and Pillow, W. (2017). Gaussian process based nonlinear latent structure discovery in multivariate spike train data. Advances in neural information processing systems, 30:3496.
[59] Yu, B. M., Cunningham, P., Santhanam, G., Ryu, S. I., Shenoy, K. V., and Sahani, M. (2009). Gaussian-process factor analysis for low-dimensional single-trial analysis of neural population activity. Journal of Neurophysiology, 102(1):614
[60] Zhao, Y. and Park, I. M. (2017). Variational latent Gaussian process for recovering single-trial dynamics from population spike trains. Neural computation, (9)5):293-11166
[61] Zhao, Y., Yates, J. I.., Levi, alignment in primate area mt. PLoS computational biology, 1605)11007
Huk, A. C.,, and Park, I. M. (2020). Stimulus-choice (mis)
[62] Zimnik, A. J. and Churchland. M. M. (2021). Independent generation of sequence elements by motor cortex. Nature Neuroscience, 24:412-424
14


Reviews:
[{"review_id": "lz6q7EFOgel", "reviewer": "Reviewer_Lhbz", "summary": "The authors propose a scalable inference algorithm for the GP factor analysis method and apply it for analysis of neural activity. The main advantages of the method are very fast execution, support for non-Gaussian likelihoods, and automatic selection of dimensionality for multivariate cases.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The GPFA model is relevant for the field, especially for the application of neural activity that fits well within the traditional scope of NeurIPS. The main contribution is also very clear, providing tangible improvement for the computational cost for its inference and hence supporting inference for longer temporal time scales. The paper also includes additional contributions related to the noise model and automatic dimensionality estimation, but these are fairly standard elements that are easy to incorporate as part of models like this. Overall, the method has potential for becoming a standard tool for the application.\n\nThe technical solution, however, is quite straightforward and limits the originality of the work; I do not see significant contributions for the more broad inference literature. The ARD prior used for estimating the dimensionality is the standard choice most researchers would start from, given its prevalent use in other FA methods. The variational approximation is based on recent literature and uses well-justified non-trivial elements like the whitened parameterization of [18] and suitable parameterization for fast computation as in Eq. (10), but it still remains more as a case study of how to best to VI for this particular model. Nevertheless, the algorithm itself seems very reasonable with carefully chosen details for efficient computation.\n\nThe empirical experiments are good, including both a clear study on synthetic data and a highly detailed application in modeling neural activity. While the application case would be already a bit too detailed for most machine learning venues, it does fit NeurIPS well. The graphical illustrations are well designed.\n\nOverall, the paper is well written, technically sound and non-trivial, and makes a clear contribution for a model of practical relevance. Even though the significance of the technical development outside the GPFA is a bit limited, I cannot think of any particular flaws in the presentation or execution.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "0.75"}, "extra_text": {}}, {"review_id": "lAJojReZVhj", "reviewer": "Reviewer_sBjz", "summary": "The authors introduce a Bayesian GPFA model, and extend this formulation to non-Gaussian noise models. A Gaussian prior over the loading matrix is proposed, and the paper employs a variational inference strategy for the inference of the model. The inference is shown to be efficient in the number of data points and helps in automatic relevance determination of the observations, as well as learning the appropriate number of latent variables to use. The method is used on a synthetic dataset to recover the dimensionality and the structure of the latent dimensions. Additionally, the method is applied to neural data from primates during a self-paced reaching task, with recordings from the primary somatosensory and motor cortices.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The authors have a clear presentation of their ideas and the methodology. The implemented variational inference method performs well on synthetic data and on spiking neural data from primates. It also leads to interesting take-aways for the neural data during the different phases of the task.\n\nPlease include an explanation on why GPFA has a lower MSE than bGPFA on synthetic data in Figure 2b with latent dimensionality = 2 (and for most other values for latent dimensionality). If the algorithm is indeed more efficient in practice, please show this.\n\nThe analysis of primate recordings using comparable methods includes the comparison between the similarity matrices of preparatory dynamics computed using raw data and FA – please comment on these same similarity matrices computed using GPFA, perhaps fit on a smaller number of trials due to computational considerations. Please also comment on the dimensionality found in the neural data using bGPFA vs. other comparable methods (for example, GPFA using the relevant noise models). Finally, please comment on the structure in the latent dimensions (as in Figure 3e) using other comparable methods.\n\nThe authors stress the possibility to perform inference on discrete (spike count) data using the proposed method. However, they do not use this capability for the electrophysiological neural data, choosing to bin the data and analyze the firing rates instead. Please explain.\n\nWhile the discoveries that are made on the neural data are interesting, further analysis may be needed to show the utility of the method against comparable methods.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "h2E2sGLsgNn", "reviewer": "Reviewer_6k1i", "summary": "\nThe authors propose a scalable inference method for Bayesian GPFA and show extensive results applying bGPFA to neural recordings from M1 and S1. They show interesting scientific findings of primate motor cortex with one reaching dataset. \n", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\nI think the last two points \"motor preparation is important in the unconstrained self-paces task\" and \"slow dynamics reflect motivation or task engagement\" are two novel findings in the neuroscience community (correct me if I'm wrong given my limited knowledge). But the motor preparation is only presented with a weak correlation, and the long timescale part has everything in the appendix. That kind of makes me doubt the true significance of these results. All other results are consistent with existing findings which are not strong scientific contributions. The modeling part is a bit weak given that neither bGPFA nor the scalable inference is novel, thus leading to an incremental contribution to the field of computational neuroscience.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "1. The data is from paper [25] which proposed a method called AutoLFADS. Their method is more advanced than GPFA and showed an extensive comparison with GPFA in [25]. I'm wondering whether bGPFA could compete against it to establish its modeling contribution. I think it might be hard but it's still worth showing. Cuz applying a well-known model to a well-known dataset and uncovering existing knowledge is not significant enough. \n\n2. The main technical novelty in this paper is scalable inference. But the inference tricks exploited in this paper are fairly common and standard in the field of GP, e.g., considering eq 5 as a GPLVM with a linear kernel. \n\n3. In the second experiment, the authors apply bGPFA to both M1 and S1 regions. The only conclusion is \"this broad peak reflects the fact that these neural populations encode both current behavior in S1 as well as future behavior in M1\". This seems to be a bit straightforward since the latent is found using both neural populations. What would be more interesting to show is what the latent would look like, how many latent dimensionalities are found using ARD (this is mentioned for individual fitting)? How the latent compares with the latent from M1 and S1 only? I think for latent dynamic/factor models for neural populations, it's more interesting to visualize more latents.\n\n\n\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2 hours"}, "extra_text": {}}, {"review_id": "-ZQVoy_mi91", "reviewer": "Reviewer_eNxg", "summary": "The authors introduce a new variant of GPFA for neural data that includes an ARD prior over the mixing matrix to automatically learn dimensionality. The authors also develop specific inference methods for this work as well as exploit structure in the covariance matrix for stability/speed. The authors test their model on synthetic data to demonstrate that the bGPFA model correctly identifies true latent dimensionality and performs well in situations with Poisson and Negative Binomial noise.  The authors also show model performance on multi-region reaching data in monkey cortex and explore the presence of time delays in the model.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "######\n\nThe authors have thoughtfully and thoroughly addressed all of my concerns. With the proposed edits, I think this paper is very appropriate for publication in this year's NeurIPS.\n\n######\n\nThis paper is clearly written and the contributions, in particular the the ARD prior as well as specifying an additional variational distribution over the noise observation , is an important advancement in this class of models used in neuroscience settings. The author's work exploring the monkey data is an interesting scientific contribution as well. However, the authors are overstating the advancements on a few points in the manuscript, and this should be corrected.\n\nFor one thing, in the abstract the authors state that their approach \"also allows for non-Gaussian noise models more appropriate for electrophysiological recordings\", and in the intro the authors state \"bGPFA also naturally extends to non-Gaussian data where it recovers ground truth parameters and latent trajectories.\" A reader here is left with the impression that a central contribution here is that GPFA has not been used in non-Gaussian settings, but the reality is that there are very many modern instances of GPFA with non-Gaussian noise. In fact, it is now more common to use Poisson GPFA than canonical GPFA for neuroscience. The authors cite two of these instances (Zhao and Park 2016 and Duncker and Sahani 2018), but there are recent inference advancements for GPFA for three discrete observation models -- Poisson, Negative-Binomial, and Binomial ( Keeley et al 2020 ICML* ), and GPFA with Poisson observations has been used in recent analysis for motor behavior (Zhao et al 2020^) and also recently used to dissect signal and noise representations in neural data (Keeley et al 2020 Neurips^*) . \n A more appropriate comparison for the evaluations in fig 2 would be one of these Poisson GPFA models, not just canonical GPFA, as we are dealing with neural data. (see Zhao and Park 2016 vLGP code here, which is quite user-friendly: https://github.com/catniplab/vlgp). \n\nMoreover, the Fourier trick used with the Toeplitz structured covariance matrix is quite similar to what is done in Keeley et al 2020 Neurips^* which represented of the variational distribution and prior covariance matrix in the Fourier domain for scalability for Poisson-observation GPFA. \n\nI find the rest of the manuscript, including the evaluation on real data  and the thorough appendix and available code, very practical and useful. The automatic selection of dimensionality is helpful, and the discussion comparing to GPLVMs quite interesting to the machine learning community. The new inference with whitened parametrization is a nice contribution, as well the proof of standard GPFA equivalence for their variational setting. \n\nIt would be nice to have some discussion as to how this method's inference compares to other approaches to GPFA with Poisson noise. \n\nIf this authors can better contextualize their achievements, and perhaps show some results with a comparison to a Poisson GPFA model (perhaps even their own, as they stated they have in their code?), I will gladly raise my score.\n\nPapers references: \n\n*Keeley, SL, et al. \"Efficient non-conjugate Gaussian process factor models for spike count data using polynomial approximations.\" International Conference on Machine Learning. PMLR, 2020.\n\n^Zhao, Yuan, et al. \"Stimulus-choice (mis) alignment in primate area MT.\" PLoS computational biology 16.5 (2020): e1007614.\n\n^*Keeley SL, Aoi MC, Yu Y, Smith SL, BR & Pillow JW (2020). Identifying signal and noise structure in neural population activity with Gaussian process factor models. Advances in Neural Information Processing Systems (NeurIPS) 33: 13795-13805.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "None", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: Two reviewers recommend acceptance and, after the discussion period, the other two reviewers don't oppose. I read the paper and I agree that although individual parts of the model and inference are known, it is their combination that brings significance to the work, especially by allowing its application to large scale data in neuroscience. I recommend Accept.
Decision: Accept (Poster)