Venue: NeurIPS 2021 Conference

Reviews:

final_score: rating: 6: Marginally above the acceptance threshold
review_content: 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.
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.

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

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

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

While 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.
limitations_and_societal_impact: Yes
review_title: null

final_score: rating: 5: Marginally below the acceptance threshold
review_content: summary: 
The 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. 

main_review: 
I 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.

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. 

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

3. 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.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: 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.
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.

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

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

Overall, 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.
limitations_and_societal_impact: Yes
review_title: null

final_score: rating: 7: Good paper, accept
review_content: 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.
main_review: ######

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

######

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

For 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^*) . 
 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). 

Moreover, 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. 

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

It would be nice to have some discussion as to how this method's inference compares to other approaches to GPFA with Poisson noise. 

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

Papers references: 

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

^Zhao, Yuan, et al. "Stimulus-choice (mis) alignment in primate area MT." PLoS computational biology 16.5 (2020): e1007614.

^*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.

limitations_and_societal_impact: None
review_title: null

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)