Scalable Bayesian GPFA with automatic relevance determination and discrete noise models
Year: 2022

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 reachspecific 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.

heading: 1 Introduction
text: 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 a range of cognitive and motor tasks have shown that neural population activity can often be accurately summarised by the dynamics of a “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 of simultaneously 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 (NeurIPS 2021).
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 O(T 3) and a memory footprint O(T 2) 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-trials” and treating these as independent samples. Additionally, canonical GPFA assumes a Gaussian noise model while recent work has suggested that 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(T log T ) and a memory cost of O(T ). To do this, we introduce an efficiently parameterized variational inference strategy that ensures scalability to long recordings while also supporting non-Gaussian noise models. Additionally, the Bayesian formulation provides a framework for principled model selection based on approximate marginal likelihoods [53]. This allows us to 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-validation to select the latent dimensionality. We then apply bGPFA to longitudinal, multi-area recordings from primary motor (M1) 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 M1 and feedback encoded by S1. We also show that the latent trajectories for M1 converge to 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 at target onset is 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 et 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.

heading: 2 Method
text: In the following, we use the notationA to refer to the matrix with elements aij . We use ak to refer to the kth row or column ofA with an index running from 1 to K, represented as a column vector.

heading: 2.1 Generative model
text: Latent variable models for neural recordings typically model the neural activity Y ∈ RN×T of N neurons at times t ∈ RT as arising from shared fluctuations in D latent variables X ∈ RD×T . Specifically, the probability of a given recording can be written as
p(Y |t) = ∫ p(Y |F ) p(F |X) p(X|t) dF dX, (1)
where F ∈ RN×T are intermediate, neuron-specific variables that can often be thought of as firing rates or a similar notion of noise-free activity. For example, GPFA [59] specifies
p(Y |F ) = ∏ n,t N (ynt; fnt, σ2n) (2)
p(F |X) = δ(F −CX) (3) p(X|t) = ∏ d N (xd; 0,Kd) withKd = kd(t, t)) (4)
That is, the prior over the dth latent function xd(t) is a Gaussian process [45] with covariance function kd(·, ·) (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, we additionally introduce a prior distribution over the mixing matrix C ∈ RN×D 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 to 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 a 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), we use a Gaussian prior over C of the form cnd ∼ N (0, s2d), where sd is a scale parameter associated with latent dimension d. Integrating C out in Equation 3 then yields the following observation model:
p(F |X) = ∏ n N (fn; 0,XTS2X), with S = diag(s1, . . . , sD). (5)
Moreover, we use a general noise model p(Y |F ) = ∏ n,t p(ynt|fnt) where p(ynt|fnt) is any distribution for which we can evaluate its density.

heading: 2.2 Variational inference and learning
text: To train the model and infer bothX and F from the data Y , we use a nested variational approach. It is intractable to compute log p(Y |t) (Equation 1) analytically for bGPFA, and we therefore introduce a lower bound on log p(Y |t) at the outer level and another one on log p(Y |X) at the inner level. These lower bounds are constructed from approximations to the posterior distributions over latents (X) and noise-free activity (F ) respectively.
Distribution over latents At the outer level, we introduce a variational distribution q(X) over latents and construct an evidence lower bound (ELBO; 55) on the log marginal likelihood of Equation 1:
log p(Y |t) ≥ L := Eq(X) [log p(Y |X)]− KL [q(X)||p(X|t)] . (6) Conveniently, maximizing this lower bound is equivalent to minimizing KL [q(X)||p(X|Y )] and thus also yields an approximation to the posterior over latents in the form of q(X). We estimate the first term of the ELBO using Monte Carlo samples from q(X) and compute the KL term analytically.
Here, we use a so-called whitened parameterization of q(X) [19] that is both expressive and scalable to large datasets:
q(X) = D∏ d=1 N (xd;µd,Σd) with µd = K 1 2 d νd and Σd = K 1 2 d ΛdΛ T dK 1 2 d T , (7)
whereK 1 2
d is any square root of the prior covariance matrixKd, and νd ∈ RT is a vector of variational parameters to be optimized. Λd ∈ RT×T is a positive semi-definite variational matrix whose structure is 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
KL[q(X)||p(X|t)] = 1 2 ∑ d ( ‖Λd‖2F − 2 log |Λd|+ ||νd||2 − T ) . (8)
Third, q(X) can be sampled efficiently via a differentiable transform (i.e. the reparameterization trick) provided that fast differentiableK 1 2
d v and Λdv products are available for any vector v:
x (m) d = K
1 2
d (νd + Λdηd) with ηd ∼ N (0, I), (9) where x(m)d ∼ q(xd). This is important to form a Monte Carlo estimate of Eq(X) [log p(Y |X)].
To avoid the challenging computation ofK 1 2 d v for generalKd [2], we directly parameterizeK 1 2
d , the positive definite square root ofK, which implicitly defines the prior covariance function kd(·, ·). In this work we use an RBF kernel forKd and give the expression forK 1 2 d in Appendix F. Additionally, we use Toeplitz acceleration methods to computeK 1 2
d v products in O(T log T ) time and with O(T ) memory cost [48, 57].
We implement and compare different choices of Λd in Appendix F. For the experiments in this work, we use the parameterization Λd = ΨdCd, where Ψd is diagonal with positive entries and Cd 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 a simple diagonal parameterization Λd = Ψd.
Distribution over neural activity Evaluating log p(Y |X) = ∑ n log p(yn|X) for each sample drawn from q(X) is intractable for general noise models. Thus, we further lower-bound the ELBO of Equation 6 by introducing an approximation q(fn|X) to the posterior p(fn|yn,X): log p(yn|X) ≥ Eq(fn|X) [log p(yn|fn)]− KL [q(fn|X)||p(fn|X)] . (10) We repeat the whitened variational strategy described at the outer level by writing
q(fn|X) = N (fn; µ̂n, Σ̂n) with µ̂n = K̂ 1 2 ν̂n and Σ̂n = K̂ 1 2LnL T n (K̂ 1 2 )T , (11) where ν̂n ∈ RD is a neuron-specific vector of variational parameters to be optimized along with a lower-triangular matrix Ln ∈ RD×D; and K̂ denotes the covariance matrix of p(f |X), whose square root K̂ 1 2 = XTS follows from Equation 5. The low-rank structure of K̂ enables cheap matrix-vector products and KL divergences:
KL[q(fn|X)||p(fn|X)] = 1
2
( ‖Ln‖2F − 2 log |Ln|+ ||ν̂d||2 −D ) . (12)
Note that the KL divergence does not depend onX in this whitened parameterization (Appendix H). Moreover, q(fn|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: Eq(fn|X) [log p(yn|fn)] = ∑ t Eq(fnt|X) [log p(ynt|fnt)] . (13) Each term in this sum is simply a 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-Hermite quadrature (Appendix K; 18).

heading: 2.3 Summary of the algorithm
text: Putting Section 2.1 and Section 2.2 together, optimization proceeds at each iteration by drawing M Monte Carlo samples {Xm}M1 from q(X) and estimating the overall ELBO as:
L = 1 M ∑ Xm∼q(X) [∑ n,t Eq(fnt|Xm) [log p(ynt|fnt)] ]
− ∑ n KL [q(fn)||p(fn)]− ∑ d KL [q(xd)||p(xd)] , (14)
where the expectation over q(fnt|X) is evaluated analytically or using Gauss-Hermite quadrature depending on the noise model (Appendix K). We maximize L using stochastic gradient ascent with Adam [29]. This has a total computational time complexity of O(MNTD2 + MDT log T ) and memory complexity of O(MNTD2) 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 t in Equation 14 are computed in multiple passes. The algorithm is described in pseudocode with further implementation and computational details in 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.1 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 sd are approximately zero for these dimensions (see Appendix I for details).

heading: 3 Experiments and results

heading: 3.1 Synthetic data
text: We first generated an example dataset from the GPFA generative model (Equations 2-4) with a true latent dimensionality of 3. We proceeded to fit both factor analysis (FA), GPFA, and bGPFA with different latent dimensionalities D ∈ [1, 10]. Here, we fitted bGPFA without automatic relevance determination such that sd = s∀d. As expected, the marginal likelihoods increased monotonically with D for both FA and GPFA (Figure 2a; Appendix I). In contrast, the bGPFA ELBO reached its optimum value at the true latent dimensionality D? = 3. This is a manifestation of “Occam’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 {sd} 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 error across GPFA and bGPFA without ARD (Figure 2a and b, black) without a 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 sd remained well above zero after training (Figure 2b, inset). Similar to 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 a 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 a 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 a linear transformation (Figure 2d). Fitting flexible noise models such as the negative binomial model is 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 in 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.
In summary, bGPFA provides a flexible method for inferring both latent dimensionalities, latent trajectories, and heterogeneous single-neuron parameters in an unsupervised manner. In the next section, we show that the scalability of the model and its interpretable parameters also facilitate the analysis of large neural population recordings.

heading: 3.2 Primate recordings
text: In this section, we apply bGPFA to biological data recorded from a rhesus macaque during a selfpaced reaching task with continuous recordings spanning 30 minutes (37, 42; Figure 3a). The continuous nature of these recordings as one long trial makes it a challenging dataset for existing analysis methods that explicitly require the availability of 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 of continuous recordings into surrogate trials can still enable the use of these methods [28], here we show that our 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, we used a single recording session (indy_20160426, as in 28), excluded neurons with overall firing rates below 2 Hz, and binned data at 25 ms resolution. This resulted in a data array Y ∈ R200×70482 (130 M1 neurons and 70 S1 neurons).
We first fitted bGPFA independently to the M1 and S1 sub-populations withD = 25 latent dimensions. In this case, ARD retained 16 (M1) and 12 (S1) dimensions (Figure 3b). We then proceeded to train a linear decoder to predict hand kinematics in the form of x and y hand velocities from either
the inferred firing rates or the raw data convolved with a 50 ms Gaussian kernel (28; Appendix M). We found that the model learned by bGPFA predicted kinematics 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 S1 predominantly encoded current behavior while neurons in M1 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 M1 functions as a dynamical system driving behavior via downstream effectors [8].
We then fitted bGPFA to the entire dataset including both M1 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 S1 as well as future behavior in M1 (Figure 3c). Indeed, when we took this offset into account by shifting all M1 spike times by +100 ms and retraining the model, decoding performance increased from 68.56% ± 0.09 to 69.81% ± 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 M1 spikes by 100 ms increased the ELBO per neuron (−34,637.0± 0.7 to −34,631.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 S1.
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 a 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). In 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 M1 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 a delay period with a known reach direction – is close to an “optimal 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,
Pearson ρ = 0.45, p = 4× 10−36). Such a correlation was also weakly present with factor analysis (ρ = 0.21, p = 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, our results suggest that motor preparation is 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.

heading: 3.3 Long-timescale latent processes
text: 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 (τ = 1.4 s) was strongly correlated with reaction time during successful reaches (Pearson ρ = 0.40, p = 3.4× 10−28). Interestingly, the information contained about reaction time in this long timescale latent dimension was largely complementary to that encoded by the distance to preparatory states in the ‘fast’ dimensions (Appendix C), suggesting
that motor preparation and task engagement are orthogonal processes both contributing to task performance.
The experimental recordings were also characterized by a 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 a different subspace as the monkey stopped engaging with the task (Figure 4b). Importantly, we were able to simultaneously capture these context-dependent changes as well as movement-specific and preparatory dynamics (Section 3.2) by fitting a single model to the full 30 minute dataset. This suggests that bGPFA can capture behaviorally relevant dynamics within individual contexts even when trained on richer datasets with changing contexts.
Finally, we wondered how the neural activity patterns during periods with and without task participation were related to the long-timescale latent dimensions predictive of task engagement. Here we found that the slow latent process considered above also exhibited a prominent change to a different state as 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 a state with no 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 ρ = 0.40 vs. ρ = 0.37). It is not surprising that task engagement decreases with time, and it 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 to represent an aspect of engagement with the task beyond the passing of time.
Taken together, we thus find that bGPFA is capable of capturing not only single-reach 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.

heading: 4 Discussion
text: Related work The generative model of bGPFA can be considered an extension of the canonical GPFA model proposed by Yu et al. [59] to include a 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, 53]. 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 electrophysiological recordings [12, 25, 34, 50, 60, 61]. While variational inference has proven a 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 is 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 on 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 a 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(Y |X) is in fact equivalent to the sparse inducing point approximation used in the stochastic variational GP (SVGP) framework [17, 18]. In particular, our 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 a 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 a deep GP [11] with two layers – the first with an RBF kernel and the second with a linear kernel. Finally, our parameterizations of the posteriors q(xd) and q(fn) 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 neuroscience that allows for regularized, scalable inference and automatic determination of the latent dimensionality as well as the use of 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-validation 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 p(Y |X) for an RBF kernel in Euclidean space [58] or on a non-Euclidean 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 a 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, it will also be interesting to extend bGPFA to non-linear observation models p(Y |X) and provide a general-purpose 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 to 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].

heading: Acknowledgements
text: We are grateful to O’Doherty et al. [42] for making their data publicly available and to 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.

Review summaries:
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.

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.