Modeling Latent Neural Dynamics with Gaussian Process Switching Linear Dynamical Systems
Authors: Amber Hu
David M. Zoltowski
Aditya Nair
David Anderson
Lea Duncker
Scott Linderman
Venue: NeurIPS 2024 poster (2024)
Keywords: gaussian process
switching
slds
neural
neuroscience
dynamics
probabilistic
time series
Primary area: probabilistic_methods
TL;DR: Gaussian Process Switching Linear Dynamical Systems maintain the locally linear interpretability of rSLDS models while naturally capturing uncertainty in dynamics.

Abstract: Understanding how the collective activity of neural populations relates to computation and ultimately behavior is a key goal in neuroscience. To this end, statistical methods which describe high-dimensional neural time series in terms of low-dimensional latent dynamics have played a fundamental role in characterizing neural systems. Yet, what constitutes a successful method involves two opposing criteria: (1) methods should be expressive enough to capture complex nonlinear dynamics, and (2) they should maintain a notion of interpretability often only warranted by simpler linear models. In this paper, we develop an approach that balances these two objectives: the Gaussian Process Switching Linear Dynamical System (gpSLDS). Our method builds on previous work modeling the latent state evolution via a stochastic differential equation whose nonlinear dynamics are described by a Gaussian process (GP-SDEs). We propose a novel kernel function which enforces smoothly interpolated locally linear dynamics, and therefore expresses flexible -- yet interpretable -- dynamics akin to those of recurrent switching linear dynamical systems (rSLDS). Our approach resolves key limitations of the rSLDS such as artifactual oscillations in dynamics near discrete state boundaries, while also providing posterior uncertainty estimates of the dynamics. To fit our models, we leverage a modified learning objective which improves the estimation accuracy of kernel hyperparameters compared to previous GP-SDE fitting approaches. We apply our method to synthetic data and data recorded in two neuroscience experiments and demonstrate favorable performance in comparison to the rSLDS.

Modeling Latent Neural Dynamics with Gaussian Process Switching Linear Dynamical Systems
Amber Hu
David Zoltowski Stanford University dzoltow@stanford.c edu
Stanford University auberhu@stanford.c
Aditya Nair
David Anderson
Caltech & Howard Hughes Medical Institute adi nair@caltech. edu
Caltech & Howard Hughes Medical Institute wuwei@caltees edu
Lea Duncker* Columbia University Id3149@columbi edu
Scott Linderman* Stanford University swl 1el1f@ttnfoff ord. erd
Abstract
Understanding how the collective activity of neural populations relates to computa- tion and ultimately behavior is key goal in neuroscience. To this end, statistical methods which describe nigh-dimensione. neural time series in terms of low- dimensional latent dynamics have played a fundamental role in characterizing neural systems. Yet, what constitutes successful method involves two opposing criteria: (1) methods should be expressive enough capture complex nonlinear dynamics, and (2) they should maintain a notion of interpretability often only war- ranted by simpler linear models. In this paper, we develop an approach that balances these two objectives: the Gaussian Process Switching Linear Dynamical System (gpSLDS). Our method builds on previous work modeling the latent state evolution via stochastic differential equation whose nonlinear dynamics are described by a Gaussian process (GP-SDEs). We propose a novel kernel function which enforces smoothly interpolated locally linear dynamics. and therefore expresses flexible yet interpretable dynamics akin those of recurrent switching linear dynamical systems (rSLDS).. Our approach resolves key limitations of the rSL.DS such as arti- factual oscillations dynamics near discrete state boundaries, while also providing posterior uncertainty estimates of the dynamics. To fit our models, we leverage a modified learning objective which improves the estimation accuracy kernel hyperparameters compared to previous GP-SDE fitting approaches. We apply our method to synthetic data and data recorded in two neuroscience experiments and demonstrate favorable performance in comparison to the rSLDS.
Introduction
Computations in the brain thought be implemented through the dynamical evolution of neural activity. Such computations are typically studied in a controlled experimental setup, where animal engaged in behavioral task with relatively few relevant variables. Consistent with this, empirical neural activity has been reported to exhibit many fewer degrees of freedom than there are neurons in the measured sample during such simple tasks [1]. These observations have driven the use of latent
"Corresponding authors.
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

variable models to characterize low- -wwdeeenniooall structure in high- sigh-dimessional neural population activity [2. 3]. In this setting, neural activity is often modeled in terms of a www-dimensionll latent state that evolves with Markovian dynamies [4-14]. It is thought that the latent state evolution is related the computation of the system, and therefore, insights into how this evolution shaped through a dynamical system can help us understand the mechanisms underlying computation [15-21].
In practice, choosing an appropriate modeling approach for a given task requires balancing two key criteria. First, statistical models should be expressive enough to capture potentially complex and nonlinear dynamics required to carry out a particular computation. On the other hand, these models should also be interpretable and allow for straightforward post-hoc analyses of dynamics. One model class that strikes this balance is the recurrent switching linear dynamical system (rSLDS) [8]. The rSLDS approximates arbitrary nonlinear dynamics by switching between a fimite number of linear dynamical systems. This leads to powerful and expressive model which maintains the interpretability of linear systems. Because of their flexibility and interpretability, variants of the iLLS have been used many neuroscience applications [19, 22-28] and are examples of general set of models aiming to understand nonlinear dynamics using compact and interpretable components [29, 30].
However, rSLDS models suffer from several limitations. First, while the rSLDS is probabilistic model, typical use cases do not capture posterior uncertainty over inferred dynamics. This makes it difficult to judge the extent to which particular features of fitted model should be relied upon when making inferences about their role in neural computation. Second, the rSLDS often suffers from producing oscillatory dynamics in regions of high uncertainty in the latent space, such as boundaries between linear dynamical regimes. This artifactual behavior can significantly impact the interpretability and predictive performance of the rSLDS. Lastly, the rSLDS does not impose smoothness or continuity assumptions on the dynamics due to its discrete switching formulation. Such assumptions are often natural and useful in the context of modeling realistic neural systems.
In this paper, we improve upon the rSLDS by introducing the Gaussian Process Switching Linear Dynamical System (gpSLDS). Our method extends prior work on the Gaussian process stochastic differential equation (GP-SDE) model. a continuous-time method that places Gaussian process (GP) prior on latent dynamics. By developing a novel GP kernel function, we enforce locally linear, interpretable structure dynamics akin to that of the rSL.DS. Our framework addresses the afore- mentioned modeling limitations of the iSLDS and contributes new class of priors in the GP-SDE model class. Our paper organized as follows. Section provides background on GP-SDE and iSLDS models. Section 3 presents our new gpSLDS model and an inference and learning algorithm for fitting these models. In Section we apply the gpSLDS to a synthetic dataset and two datasets from real neuroscience experiments to demonstrate its practical use and competitive performance. We review related work in Section 5 and conclude our paper with a discussion in Section 6.
2 Background
2.1 Gaussian process stochastic differential equation models
Gaussian processes (GPs) define nonparametric distributions over functions. They are popular choice in machine learning due to their ability to capture nonlinearities and encode reasonable prior assumptions such as smoothness and continuity [31]. Here, we review the GP-SDE, Bayesian generative model that leverages the expressivity of GPs for inferring latent dynamics [10].
Generative model In GP-SDE, the evolution of the latent state R is modeled continuous-time SDE which underlies observed neural activity y(t.) E R at time-points t, E [0, T]. Mathematically, this is expressed as
dx = f(x) + E
(ly(())) xx] = g (Car(1:) + dd)
The drift function R + R describes the system dynamies, Σ a noise covariance matrix, and dw ~ (0, diT) is Wiener process increment. Parameters C' E DxK and d E RD define an affine mapping from latent to observed space, which then passed through a pre-specified inverse link function g()).

---

A GP prior is used to model each output dimension of the dynamics f()) independently. More formally, if f() ()]] then
(a(+) sid ~ n
for k: = 1, K,
(2)
where A. (..-) )) the kernel for the GP with hyperparameters O.
Interpretability GP-SDEs and their variants can infer complex nonlinear dynamics with posterior uncertainty estimates physical systems across a variety of applications [32-34]. However, one limitation of using this method with standard GP kernels, such as the radial basis function (RBF) kernel, is that its expressivity leads to dynamics that are often challenging to interpret. Duncker et al. [10], this was addressed by conditioning the GP prior of the dynamics f()) on fixed points f(x") O and their local Jacobians J(x") = 0 (<>) and subsequently learning the fixed-point locations x² and the locally-linearized dynamies Jxx) as model parameters. This approach allows for direct estimation of key features of f()). However, due to its flexibility, it is also prone to finding more fixed points than those included in the prior conditioning, which undermines its overall interpretability.
die
2.2 Recurrent switching linear dynamical systems
The rSLDS models nonlinear dynamics by switching between different sets of linear dynamics [8]. Accordingly, it retains the simplicity and interpretability of linear dynamical systems while providing much more expressive power. For these reasons, variations the rSLDS are commonly used to model neural dynamics [19, 22-29].
Generative model
The rSLDS a discrete-time generative model of the following form:
Cz (A, bar Qs,), Ely: æx] g(C,,, d)
where dynamics switch between J distinet linear systems with parameters (Aj, bj. Qj)y This is controlled by discrete state variable 8; modeled by a multiclass logistic regression,
J), which evolves via transition probabilities
p( 8t-1, 11- x exp(w=x,-11
x;-1 + T...))
The "recurrent" nature of this model comes from the dependence of eq. (4) on latent space locations 2y, As such, the rSLDS can be understood as learning a partition of the latent space into J linear dynamical regimes seprated by linear decision boundaries. This serves as important motivation for the parametrization of the gpSL.DS, as we describe later.
Interpretability While the rSLDS has been successfully used in many applications to model nonlinear dynamical systems, it suffers from few practical limitations. First, it often produces unnatural artifacts of modeling nonlinear dynamics with discrete switches between linear systems. For example, it may oscillate between discrete modes with different discontinuous dynamics when a trajectory is near a regime boundary. Next, common fitting techniques for rSLDS models with non-conjugate observations typically treat dynamics as learnable hyperparameters rather than as probabilistic quantities [26], which prevents the model from being able to capture posterior uncertainty over the learned dynamics. Inferring posterior distribution over dynamics is especially important in many neuroscience applications, where scientists often draw conclusions from discovering key features in latent dynamics, such as fixed points or line attractors.
3 Gaussian process switching linear dynamical systems
To address these limitations of the rSLDS, we propose new class of models called the Gaussian Process Switching Linear Dynamical System (gPSLDS). The gpSLDS combines the modeling advantages of the GP-SDE with the structured flexbility of the rSLDS. We achieve this balance by designing a novel GP kernel function that defines a smooth, locally linear prior on dynamics. While our main focus on providing an alternative to the ISLDS, the gpSLDS also contributes a new prior which allows for more interpretable learning of dynamics and fixed points than standard priors the GP-SDE framework ge. the RBF kernel). Our implementation of the gPSLDDS available at: https:/ /ttp///gttt// wtt.//oooo///o.ee http/dddddemmabb...... ab/ggsdds.
3

---

1D
B
2D
C
Generative Model
samples "in
sampled dynamics trajectories
f(x)
samples from (2)
f(x)
2.5
samples from (1)
0.5 f(x)
fi =5)
single-trial states
samples Agurri
f(4)
samples from
22
point process observations
((x)
Figure 1: SSL kernel and generative model. A. ID function samples, plotted in different colors, from GPs with five kernels: two linear kernels with different hyperparameters, partition kernels for each of the two regimes, and the SSL kernel. B. (top) An example т (x) in 2D and (bottom) a sample of dynamies from a SSL kernel in 2D with π (x: as hyperparameters. The æ11 and 22- directions the arrows are given by independent 1D samples of the kernl. C. Schematic of the generative model. Simulated trajectories follow the sampled dynamics. Each trajectory observed Poisson process or Gaussian observations.
3.1
The smoothly switching linear kernel
The core innovation of our method is novel GP kernel, which call the Smoothly Switching Linear (SSL. kemne.. The SSL kernel specifies GP prior over dynamics that maintains the switching linear structure of rSLDS models, while allowing dynamics smoothly interpolate between different linear regimes.
For every pair of locations x' E RRKK the SSL kernel with J linear regimes is defined
Essa (ar, x²)=
Σ ((# T cÿ) a TT (a))
(5)
Kigan (20))))
where e, E R M RxK is a diagonal positive semi-definite matrix, of R4, and n. (a: ) 2
with Σy=1 (x) = 1. To gain an intuitive understanding of the SSL kernel, we will separately analyze each of the two terms in the summands.
( (æx, x²), is standard linear kernel which defines GP distribution over linear The first term, Min
functions [31]. The superscript j denotes that it is the linear prior on the dynamics in regime j. M controls the variance of the function's 's slope in each input dimension, and ey is such that the variance of the function achieves its minimum value of of at a Cj. We expand on the relationship between the linear kernel and linear dynamical systems in Appendix A.
The second term is what we define as the partition kernel, (1) (x, x'), which gives the gpSLDS its (xx) (x)
(x)] as parametrizing a categorical dis- tribution over J linear regimes akin to the discrete switching variables in the rSLDS. We model π (xx as multiclass logistic regression with decision boundaries {w] ((x)= = 0 |j I,..., - 1}, where Φi xz) is any feature transformation of æ. This yields random functions which are locally constant and smoothly interpolate at the decision boundaries. More formally,
switching structure. We interpret π
exp(w) ((())/)()))))))
(x( = 1 + exp(ur](aa)rr) J f=
where w. = The hyperparameter R4 controls the smoothness of the decision boundaries.
As 0, 7T TT))) approaches a one-hot vector which produces piecewise constant functions with

---

sharp boundaries, and as 7 + oo the boundaries become more uniform. While we focus on the parametrization in eq. (6) for the experiments this paper, we note that in principle any classification method can be used, such as another GP or neural network.
The SSL kernel in eq. (5) naturally combines aspects of the linear and partition kernels via sums and products of kernels, which has an intuitive interpretation [35]. The product kernel (.())
Min (z, x "e)rra (z, x') enforces linearity in regions where тy (x(x) is close to 1. Summing over J regimes then enforces linearity in each the J regimes, leading to a prior on locally linear functions with knots determined by л (2:). We note that our kernel is reminiscent of the in Pfingsten et al. [36], which uses a GP classifier as a prior for π (a)) and applies their kernel to a GP regression setting. Here, work differs in that explicitly enforce linearity in each regime and draw a novel connection to switching models like the rSLDS.
Figure IA depicts ID samples from each kernel. Figure IB shows how SSL kernel with J 2 linear regimes separated by decision boundary zî + x2 = (top) produces a structured 2D flow field consisting of two linear systems, with 21- and zz- directions determined by 1D function samples (bottom).
3.2
The gpSLDS generative model
The full generative model for the gpSLDS incorporates the SSL kernel in eq. (5) into GP-SDE modeling framework. Instead of placing a GP prior with a standard kernel on the system dynamics as in eq. (2), simply plug in our SSL kernel that
fa(+) GP(O, (inl (). for k = 1.... K,
where the kernel hyperparameters are θ (M, o. {c,1 + r]. We then sample latent states and observations according the GP-SDE via eq. (I). schematic of the full generative model is depicted in fig. IC.
Incorporating inputs In many modern neuroscience applications, external inputs to the system, such as experimental stimuli, influence latent states. To this end, we also consider extension of the model in eq. (1) which incorporates additive inputs of the form,
often interested in how
dx = (f(x) + Bu(t))dt Eldw.
(8)
where vl() E Rd a time-varying, known input signal and e RR maps inputs linearly to the latent space. The latent path inference and learning approaches presented in the following section can naturally be extended to this setting, with updates for B available in closed form. Further details are provided in Appendix B.3-B.6.
3.3 Latent path inference and parameter learning
For inference and learning the gpSLDS, we apply and extend a variational expectation- maximization (VEM) algorithm for GP-SDEs from Duncker et al. [10]. In particular, we propose modification of this algorithm that dramatically improves the learning accuracy of kernel hyperpa- rameters, which are crucial to the interpretability of the gpSLDS. We outline the main ideas of the algorithm here, though full details can be found in Appendix B.
As in Duncker et al. [10], we consider a factorized variational approximation to the posterior,
q( f. q(x)
II p(fx uk. ()q(u)) k=1
(9)
where we have augmented the model with sparse inducing points to make inference of tractable [37]. { C R and take values Standard vEM maximizes the evidence lower bound (ELBO) to the marginal log likelihood log p(y )) by alternating between updating the variational posterior q and updating model hyperpa- rameters + [38]. Using the factorization in eq. (9), will denote the ELBO L(g(z),q(1)) θ).
The inducing points are located at {zm /
f = ug-
For inference of q(xr), we follow the approach first proposed by Archambeau et al. [39] and extended by Duncker et al. [10]. Computing the ELBO using this approach requires computing variational
5

---

expectations of the SSL kernel, which we approximate using are not available in closed form. Full derivations of this step are provided in Appendix B.3. For inference of q(u), we follow Duncker et al. [10] and choose Gaussian variational posterior for q(ux)) = Wuus m Shi)), which admits closed-form updates for the mean m and covariance S given ((x) and θ. Duncker et al. [10] perform these updates before updating O via gradient ascent on the ELBO in each vEM iteration.
quadrature as they
k+
In our setting, this did not work well practice. The gpSLDS often exhibits strong dependencies between u(2) and θ which makes standard coordinate- steps in vEM prone to getting stuck in local maxima. These dependencies arise due to the highly structured nature of our GP prior; small changes in the decision boundaries { can lead to large (adverse) changes in the prior, which prevents vEM from escaping suboptimal regions of parameter space, To overcome these difficulties. we propose a different approach for learning θ: instead of fixing q(u, (u)) and performing gradient ascent on the ELBO, we optimize θ by maximizing a partially optimized ELBO,
max L(a(z)) () 81 7(1)
o- arg max
(10)
Due to the conjugacy of the model, the inner maximization can be performed analytically. This approach circumvents local optima that plague coordinate ascent on the standard ELBO. While other similar approaches exploit model conjugacy for faster vEM convergence in sparse variational GPs [37, 40] and GP-SDEs [41], our approach the first to our knowledge that leverages this structure specifically for learning the latent dynamics of a GP-SDE model. We empirically demonstrate the superior performance of our learning algorithm in Appendix C.
3.4 Recovering predicted dynamics
It is straightforward to obtain the approximate posterior distribution over f := f(x") evaluated at any new location x". Under the assumption that f only depends on the inducing points, we can use the approximation (f = I P(f Ug O)q ux) du, which can be computed in closed-form using properties of conditional Gaussian distributions. For a batch of points {>r] Jimi, this can be computed in O(N time. The full derivation can be found in Appendix B.4.
This property highlights an appealing feature of the gpSLDS over the rSLDS. The gpSLDS infers a posterior distribution over dynamics at every point in latent space, even in regions of high uncertainty. Meanwhile, as we shall see later, the rSLDS expresses uncertainty by randomly oscillating between different sets of most-likely linear dynamics, which is much harder to interpret.
4 Results
4.1
Synthetic data
We begin by applying the gPSLDS to a synthetic dataset consisting of two linear rotational systems, one clockwise and one-coumterclockwise, which combine smoothly at x1 = 0 (fig. 2A). We simulate 30 trials of latent states from SDE as in eq. (1) and then generate Poisson process observations given these latent states for D = 50 output dimensions (i.e. neurons) over T = 2.5 seconds (Fig. 2B). To initialize C' and d, we fit Poisson LDS [4] data binned at 20ms with identity dynamics. For the rSLDS, we also bin the data at 20ms. We then fit the gpSLDS and rSLDS models with J = 2 linear states using 5 different random initializations for 100 VEM iterations, and choose the fits with the highest ELBOs in each model class.
We find that the gpSLDS accurately recovers the true latent trajectories (fig. 2C) well the true rotation dynamics and the decision boundary between them (fig. 2D). We determine this decision boundary by thresholding the learned π(x) at 0.5. In addition,
obtain estimates of fixed point locations by computing the posterior probability I P Po(r) (lffe())) < c) for a small > 0: the locations æ with high probability are shaded in purple. This reveals that the gpSLDS finds high-probaber fixed points which overlap significantly with the true fixed points, denoted by stars. In comparison, both the rSLDS and the GP-SDE with RBF kernel do not learn the correct decision boundary nor the fixed points as accurately (fig. 2F-F). Of particular note, the RBF kernel incorrectly extrapolates and finds a superfluous fixed point outside the region traversed by the true latent states.

---

A
True dynamics
B
Observed data
Inferred (gDSLDS)
G
Dynamics along sampled tr trajectory ypSLoS
20
BUDS
inforeed
10
0.0
2.5
5
-25
Inferred latents & dynamics (gDSLDS)
E
Inforred latents and dynamics (GP SDE RBF)
F
Inferred latents dynamics ((SLDS)
time
H
Quant. comparisons
Dysanee JANO VOC
20
20
50
number of tnals
Figure 2: Synthetic data results. A. True dynamics and latent states used to generate the dataset. Dynamies are clockwise and counterclockwise linear systems separated by ø1 0. Two latent trajectories are shown on top of a kernel density estimate of the latent states visited by all 30 trials. B. Poisson process observations from an example trial. C. True vs. inferred latent states for the gPSLDS and rSLDS, with 95% posterior credible intervals. D. Inferred dynamics (pink/green) and two inferred latent trajectories (gray) corresponding to those in Panel A from a gPSLDS fit with linear regimes. The model finds high-probability fixed points (purple) overlapping with true fixed points (stars). E. Analogous plot to for the GP-SDE model with RBF kernel. Note that this model does not provide partition of the dynamics. F. rLLDS inferred latents, dynamics, and fixed points (pink/green dots). G. (top) Sampled latents and corresponding dynamics from the gpSLDS, with 95% posterior credible intervals. (bottom) Same, but for the rSLDS. The pink/green trace represents the most likely dynamics at the sampled latents, colored by discrete switching variable. H. MSE between true and inferred latents and dynamics for gpSLDS, GP-SDE with RBF kernel, and rSLDS while varying the number of trials in the dataset. Error bars are +2EE over 5 random initializations.
Figure 2G illustrates the differences in how the gpSLDS and the rSLDS express uncertainty over dynamics. To the left of the dashed line, we sample latent states starting from To (7, 0) and plot the corresponding true dynamics. To the right, we simulate latent states Tsamp from the fitted model and plot the true dynamics (in gray) and the learned most likely dynamics (in color) at samp- For a well-fitting model, we would expect the true and learned dynamics to overlap. We see that the gpSLDS produces smooth simulated dynamics that match the true dynamics at sa (top). By contrast, the rSLDS expresses uncertainty by oscillating between the two linear dynamical systems, hence producing uninterpretable dynamics 7aaapp (bottom). This region of uncertainty overlaps with the x1 = 0 boundary, suggesting that the rSLDS fails to capture the smoothly interpolating dynamics present in the true system.
Next, we perform quantitative comparisons between the three competing methods (fig. 2H). We find that both continuous-time methods consistently outperform the rSLDS on both metrics, suggesting that these methods are likely more suitable for modeling Poisson process data. Moreover, the gDSLDS better recovers dynamics compared to the RBF kernel, illustrating that the correct inductive bias can yield performance gains over a more flexible prior, especially in a data-limited setting.
Lastly, we note that the gpSLDS can achieve more expressive power than the rSLDS by learning nonlinear decision boundaries between linear regimes, for instance by incorporating nonlinear features into ((x) in eq. (6). We demonstrate this feature for 2D limit cycle in Appendix D.
4.2 Application to hypothalamic neural population recordings during aggression
In this section, we revisit the analyses of Nair et al. [27], which applied dynamical systems models to neural recordings during aggressive behavior in mice. To do this, we reanalyze a dataset which consists of calcium imaging of ventromedial hypothalamus neurons from a mouse interacting with two consecutive intruders. The recording was collected from 104 neurons at 15 Hz over ~343
7

---

A
ISLDS latents & dynamics
B
goSLDS lotents dynamics
C
goSLDS attractor dynamics
D
Forward simulation
20
20
10
40
1602
0.35
/copp dinvlited forward
Figure 3: Results on hypothalamic data from Nair et al. [27]. In each of the panels A-C, flow field arrow widths are scaled by the magnitude of dynamies for clarity of visualization. A. rSLDS inferred latents and most likely dynamics. The presumed location of the line attractor from [27] is marked with red box. B. gpSLDS inferred latents and most likely dynamies in latent space. Background is colored by posterior standard deviation of dynamics averaged across latent dimensions, which adjusts relative to the presence of data in the latent space. C. Posterior probability of slow points in gpSLDS, which validates line-attractor like dynamics, marked by a red box. D. Comparison of in-sample forward simulation R² between gpSLDS, rSLDS, and GP-SDE with RBF kernel. To compute this, we choose initial conditions uniformly spaced 100 time bins apart in both trials, simulate latent states steps forward according to learned dynamics (with l ranging from 100-1500), and evaluate the R? between predicted and true observations in Nassar et al. [23]. Error bars are +SE over 5 different initializations.
seconds (i.e. 5140 time bins). Nair et al. [27] found that an rSLDS fit to this data learns dynamies that form an approximate line attractor corresponding to an aggressive internal state (fig. 3A). Here, we supplement this analysis by using the gPSLDS to directly assess model confidence about this finding.
For our experiments, we z-score and then split the data into two trials, one for each distinct intruder interacting with the mouse. Following Nair et al. [27], we choose J 4 linear regimes to compare the gpSLDS and rSLDS. We choose K = 2 latent dimensions to aid the visualization of the resulting model fits; we find that even in such low dimensions, both models still recover line attractor- like dynamics. We model the calcium imaging traces as Gaussian emissions on an evenly spaced time grid and initialize C' and d using principal component analysis. We fit models with 5 different initializations for 50 vEM iterations and display the runs with highest forward simulation accuracy (as described in the caption of fig. 3).
In fig. 3A-B, we find that both methods infer similar latent trajectories and find plausible flow fields that are parsed in terms of simpler linear components. We further demonstrate the ability of the gpSLDS to more precisely identify the line attractor from Nair et al. [27]. To do this, we use the learned (f)) to compute the posterior probability of slow dynamics on a dense (80 x 80) grid of points in the latent space using the procedure in Section 3.4. The gPSLDS finds a high-probability region of slow points corresponding to the approximate line attractor found in Nair et al. [27] (fig. 3C). This demonstrates key advantage of the gpSLDS over the rSLDS: by modeling dynamics probabilistically using structured GP prior, we can validate the finding of a line attractor with further statistical rigor. Finally, we compare the gPSLDDS, rSLDS, and GP-SDE with RBF kernel using an in-sample forward simulation metric (fig. 3D). All three methods retain similar predictive power 500 steps into the future. After that, the gPSLDS performs slightly worse than the RBF kernel; however, it gains interpretability by imposing piecewise linear structure while still outperforming the rSLDS.
4.3 Application to lateral intraparietal neural recodings during decision making
In neuroscience, there is considerable interest in understanding how neural dynamics during decision making tasks support the process of making a choice [26, 34, 42-47J. In this section, we use the gpSLDS to infer latent dynamics from spiking activity in the lateral intraparictal (LIP) area of monkeys reporting decisions about the direction of motion of a random moving dots stimulus with varying degrees of motion coherence [42, 43). The animal indicated its choice of net motion direction (either left or right) via saccade. On some trials, a 100ms pulse of weak motion, randomly oriented to the left or right, was also presented to the animal. Here, we model the dynamics of 58 neurons recorded across 50 trials consisting of net motion coherence strengths in {-.512, - -128, .128, .512). where the sign corresponds to the net movement direction of the stimulus. We only consider data 200ms after motion onset, corresponding to the start of decision-related activity.
8

---

B
0332
projection onto londed dr.
-0.38
0.028
time motion
D
-0.228
-0.32
time until soccade
Figure 4: Results on LIP spiking data from a decision decision-making making task in Stine et al. [42]. A. gPSLDS inferred latents colored by coherence, inferred dynamics with background colored by most likely linear regime, and the learned input-driven direction depicted by an orange arrow. B. Projection of latents onto the ID input-driven axis from Panel A, colored by coherence (top) and choice (bottom). C. Inferred latents with 95% credible intervals and corresponding 100ms pulse input for an example trial. D. Posterior variance of dynamics produced by the gpSL.DS.
To capture potential input-driven effects, we fit version of the gpSLDS described in eq. (8) with K latent dimensions and J 2 linear regimes over 50 vEM iterations. We encoded the input signal as +1 with sign corresponding to pulse direction. In Figure 4A, we find that not only does the gpSLDS capture distinct visual separation between trials by motion coherence, but it also learns a precise separating decision boundary between the two linear regimes. Our finding is consistent with previous work on a related task, which found that average LIP responses can represented by a 2D curved manifold [48], though here we take a dynamical systems perspective. Additionally, our model learns an input-driven effect which appears to define separating axis. To verify this, we project the inferred latent states onto the ID subspace spanned by the input effect vector. Figure 4B shows that the latent states separate by coherence (top) and by choice (bottom), further suggesting that the pulse input relates to meaningful variation in evidence accumulation for this task. Fig. 4C shows an example latent trajectory aligned with pulse input; during the pulse there is noticeable change the latent trajectory. Lastly, in Fig. 4D we find that the gpSLDS expresses high confidence in learned dynamics where latent trajectories are present and low confidence in areas further from this region.
5
Related work
There are several related approaches to learning nonlinear latent dynamics in discrete or continuous time. Gaussian process state-space models (GP-SSMs) [49-53] can be considered a discrete-ting time analogue to GP-SDEs. In a GP-SSM, observations are assumed to be regularly sampled and latent states evolve according to a discretized dynamical system with a GP prior. Wang et al. [49] and Turner et al. [50] learned the dynamics of a GP-SSM using maximum posteriori estimation. Frigola et al. [51] and Eleftheriadis et al. [52] employed a variational approximation with sparse inducing points to infer the latent states and dynamics in a fully Bayesian fashion. In our work, we use continuous-tine framework that more naturally handles irregularly sampled data, such as point-process observations commonly encountered in neural spiking data. Neural ODEs and SDEs [54, 55) use deep neural networks to parametrize the dynamics of a continuous time system, and have emerged as prominent tools for analyzing large datasets, including those in neuroscience [33, 56-58]. While these methods can represent flexible function classes, they are likely to overfit to low-data regimes and may be difficult to interpret. In addition, unlike the gpSLDS, neural ODEs and SDEs do not typically quantily uncertainty of the learned dynamics.
In the context of dynamical mixture models, Köhs et al. [59, 60] proposed continuous tittiuuuu---imeeeee switching model a GP-SDE framework. This model assumes a latent Markov jump process over time which controis the system dynamics by switching between different SDEs. The switching process models dependence on time, but not location latent space. contrast, the gPSLDD does not explicitly represent a latent switching process and rather models switching probabilities as part of the kernel

---

function. The dependence of the kernel on the location in latent space allows for the gpSLDS to partition the latent space into different linear regimes.
While our work has followed the inference approach of Archambeau et al. [61] and Duncker et al. [10], other methods for latent path inference in nonlinear SDEs have been proposed [32, 41, 62, 63]. Verma et al. [41] parameterized the posterior SDE path using an exponential family-based description. The resulting inference algorithm showed improved convergence of the E-step compared to Archambeau et al. [39]. Course and Nair [32, 62] proposed an amortization strategy that allows the variational update of the latent state to be parallelized over sequence length. In principle, any of these approaches could be applied to inference in the gpSLDS and would be an interesting direction for future work.
6 Discussion
In this paper, we introduced the gpSLDS to infer low-dimensional latent dynamical systems from high-dimensional, noisy observations. By developing a novel kernel for GP-SDE models that defines distributions over smooth locally-linear functions, we were able to relate GP-SDEs to rSLDS models and address key limitations of the rSLDS. Using both simulated and real neural datasets, demonstrated that the gpSLDS can accurately infer true generative parameters and performs favorably in comparison to rSLDS models and GP-SDEs with other kernel choices. Moreover, our real data examples illustrate the variety of potential practical uses of this method. On calcium imaging traces recorded during aggressive behavior, the gPSLDS reveals dynamics consistent with the hypothesis of a line attractor put forward based on previous rSLDS analyses [27]. On a decision making dataset, the gPSLDS finds latent trajectories and dynamics that clearly separate by motion coherence and choice, providing dynamical systems view consistent with prior studies [42, 48].
In our experiments, we demonstrated the ability of the gpSLDS to recover ground-truth dynamical systems and key dynamical features using fixed settings of hyperparameters: the latent dimensionality K and the number of regimes J. For simulated data, we set hyperparameters real data, we chose hyperparameters based on prior studies and did not further optimize these values, However, for most real neural datasets, we do not know the true underlying dimensionality optimal number of regimes. To tune these hyperparameters, we can resort to standard techniques for model comparison in the neural latent variable modeling literature. Two common evaluation metries forward simulation accuracy [23, 27] and co-smoothing performance [6, 64, 65].
their true values; for
While these results are promising, we acknowledge few limitations of the gpSLDS. First, the memory cost scales exponentially with the size of the latent dimension due to using quadrature methods to approximate expectations of the SSL kernel, which are not available in closed form. This computational limitation renders difficult to fit the gpSLDS with many (i.e. greater than 3) latent dimensions. One potential direction for future work would be to instead use Monte Carlo methods to approximate kernel expectations for models with larger latent dimensionality. In addition, while both the gPSLDS and rSLDS require choosing discretization timestep for solving dynamical systems, in practice we find that the gpSLDS requires smaller steps for stable model inference. This allows the gpSLDS to more accurately approximate dynamics with continuous-time likelihoods, the cost of allocating more time bins during inference. Finally, we acknowledge that traditional variational inference approaches such as those employed by the gPSLDS- tend to underestimate posterior variance due to the Kl-divergence-baat based objective [38]. Carefully assessing biases introduced by our variational approximation to the posterior would be an important topic for future work.
Overall, the gpSLDS provides a general modeling approach for discovering latent dynamics of noisy measurements in intepretable and fully probabilistic manner. We expect that our model will be useful addition to the rSLDS and related methods on future analyses of neural data.
Acknowledgments and Disclosure of Funding
This work was supported by grants from the NIH BRAIN Initiative (U19NS11320) ROINS131987, & RFIMM133778) and the NSF/NIH CRCNS Program (ROINSI (R01NS130789), as well as fellowships from the Simons Collaboration on the Global Brain, the Wu Tsai Neurosciences Institute, the Alfred P. Sloan Foundation, and the McKnight Foundation. We thank Barbara Engelhardt, Julia Palacios, and the members of the Linderman Lab for helpful feedback throughout this project. The authors have no competing interests to declare.
10

---

References
[1] Peiran Gao and Surya Ganguli. On simplicity and complexity in the brave new world of large-scale neuroscience. Current Opinion in Newrobiology. 32: 52:148-155 2015.
[2] John P Cunningham and Byron M Yu. Dimensionality reduction for large-scale neural record- ings. Nature Neuroscience, 17(11):100-1509, 2014.
[3] Shreya Saxena and John P Cunningham. Towards the neural population doctrine. Current Opinion in Neurobiology, 55:103-111, 2019.
[4] Anne C Smith and Emery N Brown. Estimating a space model from point process observations. Neural Computation, 15(5):965-99)) 2003.
[5] Liam Paninski, Yashar Ahmadian, Daniel Gil Ferreira, Shinsuke Koyama, Kamiar Rahnama Rad, Michael Vidne, Joshua Vogelstein, and Wei Wu. A new look state-space models for neural data. Journal of Computational Neuroscience, 29:107-126, 2010.
[6] Jakob H Macke, Lars Bucsing, John P Cunningham, Byron M Yu, Krishna V Shenoy, and Maneesh Sahani. Empirical models of spiking in neural populations. Advances in Neural Information Processing Systems, 24, 2011.
[7] Evan Archer, Memming Park, Lars Buesing, John Cunningham, and Liam Paninski. Black box variational inference for state space models. arXiv preprint arXiv:1511.0736 2015.
[8] Scott Linderman, Matthew Johnson, Andrew Miller, Ryan Adams, David Blei, and Liam Paninski. Bayesian learning and inference in recurrent switching linear dynamical systems. In Artificial Intelligence and Statistics, pages 914-922. PMLR, 2017.
[9] Chethan Pandarinath, Daniel J O' Shea, Jasmine Collins, Rafal Jozefowicz, Sergey Stavisky, Jonathan C Kao, Eric M Trautmann, Matthew T Kaufman, Stephen Ryu, Leigh R Hochberg, al. Inferring single-trial trial neural population dynamics using sequential auto-encoders. Nature Methods, 15(10)8005-81) 2018.
[10] Lea Duncker, Gergo Bohner, Julien Boussard, and Maneesh Sahani. Learning interpretable continuous-timm models of latent stochastic dynamical systems. In International Conference on Machine Learning, pages 1726-1734 PMLR, 2019.
[11] Marine Schimel, Ta-Chu Kao, Kristopher T Jensen, and Guillaume Hennequin. ilqr-vae: control- based learning of input-driven dynamics with applications neural data. In International Conference on Learning Representations, 2021.
[12] Mikhail Genkin, Owen Hughes, and Tatiana A Engel. Learning non-stationary langevin dynamics from stochastic observations of latent trajectories. Nature Communications, 12(1): 5986, 2021.
[13] Christopher Langdon and Tatiana A Engel. Latent eircuit inference from heterogeneous neural responses during cognitive tasks. bioRxiv, pages 2022-01, 2022.
[14] Matthew Dowling, Yuan Zhao, and II Memming Park. Large-scale variational gaussian state- space models. arXiv preprint ar arXiv:2403.01371, 2024.
(15) Saurabh Vyas, Matthew D Golub, David Sussillo, and Krishna V Shenoy. Computation through neural population dynamics. Annual Review of Neuroscience, 43:249-275, 2020.
[16] Lea Duncker and Maneesh Sahani. Dynamics on the manifold: Identifying computational dynamical activity from neural population recordings. Current Opinion in Neurobiology, 70: 163-170, 2021.
[17] David Sussillo and Omri Barak. Opening the black box: low-dimensional dynamics in high- dimensional recurrent neural networks. Neural Computation, (5(3))626-644 2013.
[18] Daniel O'Sa,a, Lea Duncker, Werapong Goo, Xulu Sun, Saurabh Vyas, Eric Trautmann, Ilka Diester, Charu Ramakrishnan, Karl Deisseroth, Maneesh Sahani, et al. Direct neural perturbations reveal dynamical mechanism for robust computation. bioRxiv, pages 2022-12, 2022.

---

[19] Joana Goldado-Magrane Valerio Mante, and Maneesh Sahani. Inferring context context-dependent computaaiions through linear approximations of prefrontal cortex dynamics. Science Advances, 10($1))ead14743, 2024.
[20] Amit Vinograd, Aditya Nair, Joseph Kim, Scott W Linderman, and David J Anderson. Causal evidence of a line attractor encoding an affective state. Nature, pages 1-3, 2024.
[21] George Mountoufaris, Aditya Nair, Bin Yang, Dong- Kim, Amit Vinograd, Samuel Kim, Scott W Linderman, and David Anderson. A line attractor encoding persistent internal state requires neuropeptide signaling. Cell, 187(21):5998-6015, 2024.
[22] Biljana Petreska, Byron M Yu, John P Cunningham, Gopal Santhanam, Stephen Ryu, Krishna V Shenoy, and Maneesh Sahani. Dynamical segmentation of single trials from population neural data. Advances in Neural Information Processing Systems, 24, 2011.
[23] Josue Nassar, Scott Linderman, Monica Bugallo, and II Memming Park. Free-structured recur- rent switching linear dynamical systems for multi-scale modeling. In International Conference on Learning Representations, 2018.
[24] Jalil Taghia, Weidong Cai, Srikanth Ryali, John Kochalka, Jonathan Nicholas, Tianwen Chen, and Vinod Menon. Uncovering hidden brain state dynamics that regulate performance and decision-make making during cognition. Nature Communications, 9(1):2505, 2018.
[25] Antonio C Costa, Tosif Ahamed, and Greg Stephens. Adaptive, locally linear models of complex dynamics. Proceedings of the National Academy of Sciences, 11615)1100--1510, 2019.
[26] David Zoltowski, Jonathan Pillow, and Scott Linderman. A general recurrent state space framework for modeling neural dynamics during decision decision-making. In International Conference on Machine Learning, pages 1180-111911 PMLR, 2020.
[27] Aditya Nair, Tomomi Karigo, Bin Yang, Surya Ganguli, Mark J Schnitzer, Scott W Linderman, David J Anderson, and Ann Kennedy, An approximate line attractor the hypothalamus encodes an aggressive state. Cell, 186(1):178-19)))) 2023.
[28] Mengyu Liu, Aditya Nair, Nestor Coria, Scott W Linderman, and David Anderson. Encoding of female mating dynamics by a hypothalamic line attractor. Nature, pages 1-3, 2024.
[29] Jimmy Smith. Scott Linderman, and David Sussillo. Reverse engineering recurrent neural networks with Jacobian switching linear dynamical systems. Advances in Neural Information Processing Systems, 34: 1700666671111 2021.
[30] Noga Mudrik, Yenho Chen, Eva Yezerets. Christopher J Rozell, and Adam S Charles. Decom- posed linear dynamical systems (dlds) for learning the latent components of neural dynamics. Journal of Machine Learning Research, 25(59): 2024.
[31] Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
[32] Kevin Course and Prasanth Nair. State estimation of a physical system with unknown governing equations. Nature, 622(7982):211-267, 2023.
[33] Timothy Doyeon Kim, Thomas Zhihao Luo, Tankut Can, Kamesh Krishnamurthy, Jonathan W Pillow, and Carlos D Brody. Flow-field inference from neural data using deep recurrent networks. bioRxiv, 2023.
[34] Thomas Zhihao Luo, Timothy Doyeon Kim, Diksha Gupta, Adrian G Bondy, Charles D Kopec, Verity A Elliot, Brian DePasquale, and Carlos D Brody. Transitions in dynamical regime and neural mode underlie perceptual decision decision-making making. bioRxiv, pages 2023- 10, 2023.
[35] David Duvenaud. Automatic model construction with Gaussian processes. PhD thesis, 2014.
[36] Tobias Pfingsten, Malte Kuss, and Carl Edward Rasmussen. Nonstationary gaussian process re- gression using a latent extension of the input space. Fighth World Meeting of the International Society for Bayesian Analysis (ISBA 2006), 2006.
12

---

[37] Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial Intelligence and Statisties, pages 567-574. PMLR, 2009.
[38] David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: review for statisticians. Journal of the American Statistical Association, 1 11((518):859 877, 2017.
[39] Cedric Archambeau, Dan Corford, Manfred Opper, and John Shawe-Taylor Gaussian process approximations of stochastic differential equations. In Gaussian Processes Practice, pages 1-16. PMLR, 2007.
[40] Vincent Adam, Paul Chang, Mohammad Emtiyaz E Khan, and Amn Solin. Dual parameteri- zation of sparse variational gaussian processes. Advances Neural Information Processing Systems, 34:1 1174-111484 2021.
[41] Prakhar Verma, Vincent Adam, and Arno Solin. Variational gaussian process diffusion processes. In International Conference on Artificial Intelligence and Statisticx, pages 1909-1917. PMLR, 2024.
[42] Gabriel M Stine, Eric M Trautmann, Danique Jeurissen, and Michael N Shadlen. A neural mechanism for terminating decisions. Neuron, 111(166)201-221,,, 2023.
[43] Jamie D Roitman and Michael Shadlen. Response of neurons in the lateral intraparietal area during combined visual discrimination reaction time task. Journal of Neuroscience, 22(21): 9475 1115-9949,220 2002.
[44] Jochen Ditterich. Stochastic models of decisions about motion direction: behavior and physiol- ogy. Neural Networks, 19(8):9111012, 2006.
[45] Joshua I Gold and Michael N Shadlen. The neural basis of decision making. Annu. Rev. Neurosci., 30:535-574, 2007.
[46] Kenneth W Latimer, Jacob Yates, Miriam LR Meister, Alexander C Huk, and Jonathan W Pillow. Single-trial spike trains in parietal cortex reveal discrete steps during decision-making. Science, 34966244)1184-887, 2015.
[47] Mikhail Genkin, Krishna Shenoy, Chandramouli Chandrasekaran, and Tatiana A Engel. The dynamics and geometry of choice in premotor cortex. bioRxiv, 2023.
[48] Gouki Okazawa, Christina E Hatch, Allan Mancoo, Christian K Machens, and Roozbeh Kiani. Representational geometry perceptual decisions the monkey parietal cortex. Cell, 184(14); 3748 3761, 2021.
[49] Jack Wang, Aaron Hertzmann, and David J Fleet. Gaussian process dynamical models. Advances Neural Information Processing Systems, 18, 2005.
[50] Ryan Turner, Mare Deisenroth, and Carl Rasmussen. State-space inference and learning with gaussian processes. In Artificial Intelligence and Statistics, pages 868-875 JMLR Workshop and Conference Proceedings, 2010.
[51] Roger Frigola, Yutian Chen, and Carl Edward Rasmussen. Variational gaussian process state- space models. Advances in Neural Information Processing Systems, 27, 2014.
[52] Stefanos Eleftheriadis, Tom Nicholson, Mare Deisenroth, and James Hensman. Identification of gaussian process state space models. Advances in Neural Information Processing Systems, 30, 2017.
[53] Gergo Bohner and Maneesh Sahani. Empirical fixed point bifurcation analysis. arXiv preprint arXiv:1807.01486. 2018
[54] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in Neural Information Processing Systems, 31, 2018.
[55] Yulia Rubanova, Ricky TQ Chen, and David K Duvenaud, Latent ordinary differential equations for irreuularl-samppled time series. Advances in Neural Information Processing Systems, 32, 2019.
13

---

[56] Timothy D Kim, Thomas Z Luo. Jonathan W Pillow, and Carlos D Brody. Inferring latent dynamics underlying neural population activity via neural differential equations. In International Conference Machine Learning, pages 5551-5561. PMLR, 2021.
[57] Timothy Doyeon Kim, Tankut Can, and Kamesh Krishnamurthy. Trainability, expressivity and interpretability in gated neural odes. In International Conference on Machine Learning, pages 16393-16423, 2023.
[58] Christopher Versteeg, Andrew R Sedler, Jonathan D McCart, and Chethan Pandarinath. Ex- pressive dynamics models with nonlinear injective readouts enable reliable recovery of latent features from neural activity. arXiv preprini urXiv:2309.06402 2023.
[59] Lukas Köhs, Bastian Alt, and Heinz Koeppl. Variational inference for continuous-time switching dynamical systems. Advances in Neural Information Processing Systems, 0-220545-20555 2021.
[60] Lukas Köhs, Bastian Alt, and Heinz Koeppl. Markov chain monte carlo for continuous time switching dynamical systems. In International Conference on Machine Learning, pages 1430- 1 11454, PMLR, 2022.
[61] Cédric Archambeau, Manfred Opper, Yuan Shen, Dan Cornford, and John Shawe-Taylow Variational inference for diffusion processes. Advances in Neural Information Processing Systems, 20, 2007.
[62] Kevin Course and Prasanth Nair. Amortized reparametrization: Efficient and scalable variational inference for latent sdes. Advances in Neural Information Processing Systems, 36, 2024.
[63] Botond Cseke, David Schnoerr, Manfred Opper, and Guido Sanguinetti. Expectation propagation for continuous time stochastic processes. Journal of Physics A: Mathematical and Theoretical, 49(49):494002 2016
[64] Anqi Wu, Stan Pashkovski, Sandeep R Datta, and Jonathan W Pillow. Learning a latent manifold of odor representations from neural responses in piriform cortex. Advances in Neural Information Processing Systems, 31, 2018.
[65] Stephen Keeley, Mikio Aoi, Yiyi Yu, Spencer Smith, and Jonathan W Pillow. Identifying signal and noise structure in neural population activity with gaussian process factor models. Advances in Neural Information Processing Systems, 33:13795-13805 2020.
14

---

Supplementary Material
Table of Contents
Appendix A: Relationship between linear kernel and linear dynamical systems
Appendix B: Inference and learning
-Appendix C: Empirical results for new learning objective
Apppendix D: Additional synthetic data results
15

---

A Relationship between linear kernel and linear dynamical systems
Here, we draw mathematical connection between functioos sampled from a GP with linear kernel, Kiin (I. x') (x - c) MM c) o
f(-) GP(O, (in(-:)))
(11)
and functions of the form
f(z) - GT Bg-
(12)
This connection is useful for understanding the relationship between our model formulation and the typical linear dynamical systems formulation as in the rSLDS. In particular, we will show that the linear kernel in eq. (11) equivalently places prior on / and fo eq. (12), in similar manner to Bayesian linear regression.
By definition of GP, for any finite of input locations {0;
have
(f()) ((x2)
[(xN)] ~0, 4)
(13)
where Pij = (a;, m.y). Equivalently, for every pair i. I....,,
Cov(J(as); ((e,)) = (z. c) M (a; - c) c) +o0 =x; Mæj c. M(; + æz)) + c c.c+0o +0
(14)
(15)
Under eq. (12), treating ,3 and Bo as random, this would become
Cov(BT x; + fa, (TT,, + 8o) x. C0v(3,33) + Cov(On, (x;; + mz)) + Var(3o)
(16)
Equation (15) and eq. (16) are equivalent and eq. (13) satisfied if and only
3 3n
~ N (0.
M
MM c/c og
-c
M
(17)
This shows that a GP with linear kernel is equivalent to Bayesian linear regression model with prior on coefficients the form eq. (17).
B Inference and learning
We present the full inference and learning details for the gpSLDS with additive inputs,
dx (f(x) + dv( Σ dw,,
E(y(4) æ] = g (Cx(t,) d).
(18)
Our approach primarily follows that of Duncker et al. [10]. We the following notation from Duncker al. [10] throughout this section:
(h(-))pp) denotes the expectation of h(-)) with respect to the distribution p(-).
Knz E R denotes the covariance matrix defined by kernel (=at(;;)) for two batches of points {x, ) and { Nz Specifically, [Kzl (as zy). If one of these batches has only one point, c.g. N1 1, we denote the corresponding covariance vector as k., E R. IxNN Note, these quantities depend on O through the kerne;; omit writing this dependence for brevity.
B.1 Augmenting the generative model
Following Duncker et al. [10], we first augment our generative model with sparse inducing points
{u C A at input locations (=m LM C DRD Inducing points can be seen as pseudo-
observations of f()) at input locations {zm J They allow for tractable inference of f at any new batch of N input locations by reducing computational complexity from O(N3) to O M²), where typically M is chosen to be much smaller than [37]. After this augmentation, the joint likelihood of our model becomes
K
ply. x. f.u () = p(y x)p(x f) II p(fa ux. O)p(us ().
(19)
16

---

Treating uy as pseudo-observations we assume the following augmented prior:
p )) = (u 0, 5.)).
(20)
Then view [x(-) uty.
GP conditioned uks
fi(()) uk GP( (.), Mislue (::)),
(21)
where
Pfauue (a)) = kig K2 21 fuluus (z,, x²) = Kust Kasa xx') kuz K2g Rzn.
B.2
Variational lower bound
As in Duncker et al. [10], we consider variational approximation the posterior of the form
K
(ax,, f, = q(a)
II
p(/x uAx. ()q(u))
(22)
k=1
Using this factorization, we derive the ELBO to the marginal log-likelihood of our model. By Jensen's inequality,
p ()) = log / p(y x)p(< f)p(f u, p)p(u O)dxdfdu
/ g(ar, f. (u log p(y x))(<r ff {f u. O)p(u +) fuu
((a, f, u)
/ g(xx, f. u) log p(y x)p(< II p(ux ())
fxd
((x) II ((ux)
K
(logp(y) ()))(())) (KL|g(=)]lb(() | C(g(z), {()), θ),
Σ k=1
KL (x))l((uu [))]
where
(f) II / p(fa un O)q(un)dux=
(23)
B.3 Inference latent paths
To perform inference over the posterior of latent paths g(x:), we follow a method first proposed in Archambeau al. [39] and extended by Duncker al. [10].
As in Archambeau et al. [39], we choose a posterior distribution g(x) characterized by Markov Gaussian process,
a(x)) {dæ (-A(1)a(()) + b(()) di + Eidw. ~ ~~(m
(24)
=f,(w)
This distribution satisfies posterior marginals g(xx:) = -Nax; m, St) which satisfy the differential equations
dm A(1)m(() b(1)
(25)
dl
dS((t) -A())((()) S(()A(() + Σ
(26)
dl
Archambeau al. [39] maximize the ELBO with respect ((x) subject the constraints in eq. (25) and eq. (26) using the method of Lagrange multipliers. They show that after applying integration by parts, the Lagrangian becomes
Ż= = L(g(()) θ) Ci C2
(27)
17

---

where
C c= ( b(t) dA(t) m(t) dt A(T) m(T) X(0) m(0) a(0)mm(0)
(28)
T
(1)(t)
C2 = (()(((())())) S()) Σ)
dt (t)
(29)
(T)) (r)
(30)
As in Archambeau al. [39]. we assume that X(T) 0 and (T)
To find the stationary points of the Lagrangian, we first take derivatives of E with respect to m0),, S()), m()), ((), ((t, b(t) and set them to 0. The derivatives with respect to m()) and (
lead to the updates
(0) (0) ((O)((0),
(()) (2 (0) (())-
(31)
where we assume prior on initial conditions p(0o) N(xx) (()), (
The derivatives with respect to m(t) and (t) lead to the stationary equations
dA(t) A(t) X(t) A(t)
ac
(32)
dl
dm((1
d (t) = A(t) T 1 (t))((t)
Oc 0s(t)
O
(33)
dt
with P.y = for j and 1 otherwise. The inclusion of P was proposed by Duncker et al. [10] to adjust for taking derivatives with respect to symmetric matrix.
To take derivatives with respect A(t) and b(t), first extend result from Appendix A of Archambeau et al. [39] to our affine inputs model. This allows us to rewrite the KL-divergence term between the posterior and prior latent paths in the ELBO
(KL[g(x)]lp(A ff)) = 2 ('' (If + ((f+Bu(t) fq) + Bu(t) fa)) g(a)) df
L. ((Bu(t) Af) +f)) (()))1(1) dt
T
J.
((AF) (())))()) de
(34)
+ ''"' u(t) BT (Af) (Af) dt + (' u(t) B Bu(t)dt
where Af: f - -fq- The integrand of the first term in eq. (34) can be expanded as
0s
((A)) ((0))011)))))) ( f) (()))))))))) + 2Tr
A(I)
S(L)
(35)
dx
((z).o(f)
+ Tr ((t)((())) m())m(()) 2t) A()) (f) (f))())))))))))
+ b Tb(t) (f)g(a))g(s) - 26 A())m(t)
where we have used the identity (ff())) oir) (x -m))) g(w) derived from Stein's 's lemma. Note that computing eq. (35) relies three quantities,
(ef
S. which can be
q()))a(f)))
0f
(g ( (())))
Ox
v(m).v(f)
which can be written terms which depend on expectations of the kernel with respect to 9(z). We derive these as follows:
( (f)g()).015) = (kiz K--lu u)) q(())
(kia) g(m)) K K_ m
18

---

K
(f
(fa(z))) al(c))
g(m))
Var [fff + ( (((x))
x(x))
K
+ (/ar( plfa mu))) g(m))
((a)).((a)
Term
Term 2
(z))) w)) olac)
Term
(kax kx Ked kee) g(w)) + (Kee K S K koo) g(ar)
k=1
Kn s2 "m (m TT K kim g(20)
Torm
0f
=(m, K
Ukez die
dx
o(xc)
a(a))((())
m! K
0k.n ((x))
The above three function expectations can thus be expressed in terms the following four kernel expectations with respect to g(x)) N2 m, S):
@ri(z, x))
(m(#, ()) (<(2, 9)))
(x(Z),))<(z, x)),,,,)))))),1)))))))))))) 2())
dx
((x)
For our SSL kernel these kernel expectations not available in closed form, so in practice approximate them using Mauss-Hermite quadrature.
Next, differentiating Ž with respect to A() and b(l) yields the updates
0f
A(t)
2
(36)
Ox
a(t).a(f)
b(t) (f(())))(())) + A(t)m(t) + Bu(t) - -x(()
(37)
Note that these updates have one key difference from the previously derived updates in Archambeau Bu(() in eq. (37). Intuitively, this is because the posterior bias term b(t) is fully time-varying, so it captures changes in the latent states due to input-driven effects in the posterior.
al. [39] and Duncker et al. [10]: the input-dependent
In summary, the inference algorithm for updating ((x) is as follows. In each iteration of vEM, we repeat the following forward-backward style algorithm:
Solve for mt(t), S(t) forward in time starting from m()). S()) via eq. (25) and eq. (26). 2. Solve for X(t), T backward in time starting from A(T), (T) via eq. (32) and eq. (33).
3. Update A(I) and b(l) via eq. (36) and eq. (37).
After solving these stationary equations, we update m(0) and S(0) via eq. (31).
Computational details Solving for ml)), S(),, A()), and (i) requires integrating continuous- time ODEs. In practice we use Fuler integration with a small discretization step A/ relative to the sampling rate of the data, though principle any ODE solver can used. We found that the ELBO iterations.
usually converges within 20
19

---

The ODEs for solving A(t) and L (t) in eq. (32) and eq. (33) depend on evaluating gradients of the ELBO with respect (t) and Stt). We use modern autodifferentiation capabilities in JAX to compute these gradients.
B.4 Updating dynamics and hyperparameters with a modified learning objective
As we describe in Section 3.3, our gpSLDS inference algorithm uses modified objective for hyperparameter learning. In this section, we discuss this objective in detail and present closed-form updates for the inducing points given the hyperparameters.. Then, using the inducing points, we will derive the posterior distribution over f-)) at any location in the latent space.
After updating the latent paths g(x) described in Appendix B.3, update hyperparameters θ using a partially optimized ELBO. This update can be written as
o- arg max max L(g(2)) 01) 3 q(u)
(38)
Following Duncker et al. [10], we choose the variational posterior
((u)) = Nuss m. SW)
(39)
Given ((x) and θ this leads to the closed-form updates,
5h K: K + [' (keakkaz) K.:
(40)
(" ((bez))g(ce) (-A(t)m(()) + b(t) - Bu(())
m = Shia K-1 22
=S
(41)
( /Ōkgz
((I)A(1) de
ax
o((x)
In the above equation, m E M contains m E R M in each column. The inside maximization of eq. (38) can be computed analytically using these closed-form updates. Note that m and Sk* depend on θ through the prior kernel covariances K., and kzz- Therefore, eq. (38) can be understood as performing joint optimization of the ELBO with respect θ through m and Sh* as well as through the rest of the ELBO. In practice, this allows vEM to circumvent dependencies between (uu) and Ô, leading to more accurate estimation of both quantities. For our experiments, we use the Adam optimizer to solve eq. (38).
After obtaining @* in each vEM iteration, explicitly update q(u)) using eq. (41) and eq. (40) for the next iteration.
Recovering predicted dynamics Given (updated) variational parameters m and St. is straight- forward to compute the posterior distribution of f* :- I (x*) at any location x" in the latent space. Recall the variational approximation from eq. (23). If we apply this approximation to f", we have
K
K
(f
II
II k=1
plfk uy, ()q(ux)duu
(42)
To evaluate this analytically, we use properties of conditional Gaussian distributions. First note that by our augmented prior,
p(ft uk θ) = Nif Mu.:Kak K= uk, Uisalx" xx" ki K ker
(43)
Then, by conjugacy of Gaussian distributions,
a(f)) I N K Uks Kisal (x0 x²") k Bor keo N(ux m Shid
= =ft KarKK m (2s xx)))))))"))")" k K Keo + & K s K kew)). (44)
20

---

B.5
Learning observation model parameters
For the experiments in this paper, we considered two observation models: Gaussian observations and Poisson process observations.
Gaussian observations
We consider the observation model
p x) II (V((()) Cz((t) + d, R).
(45)
where R E RD is a diagonal covariance matrix. The expected log-likelihood is available in closed form and is given by
((gg)((()() C(() d, (T(()(
(logply m)) g(a)
(46)
Closed-form updates for C, d and R are also available:
C Σ (v(()) - d)m(())
[((()))m())m)
(47)
((v(t)) - C"m((t))
(48)
71.
R
(ua((t))2 2ya(ti)cmmlu (c + (T,))))))))
(49)
74.
where nt is the number of observed time points, R2 is the d-th entry R, and cd is the d-th row of C.
Poisson process observations observations of the form
The second observation model we consider is Poisson process
p(4) x) (() = PP()(CP((( dd),
(50)
where either g(a) exp(a) (exponential inverse link) g(a) log(1 + exp(a)) (softplus inverse link). For the exponential inverse link, the expected log-likelihood available in closed form and given by
(logp(((t) [x))) x))(())) = - f/ exp (mm(t) + ddigg(CS())CT) + ((Cm(()) + d)
4
(51)
For the softplus inverse link, the expected log-likelihood not available in fooed-formm but can be approximated by Gauss-Hermite quadrature second-order Taylor expansion around m(t).
For both Poisson process models, we update C' and d using gradient ascent the expected log- likelihood with the Adam optimizer.
B.6 Learning the input effect matrix
Here we derive closed- form update for B, which linearly maps external inputs to the latent space. The only term in the ELBO which depends B is (KL[g(z) [p(< fil)yor> We differentiate this term written in eq. (34) and arrive the update
B' (f- (A(aaaatsaa + A(t)m(tt ((t))o(t)Tal
(1"" e(t))(t)TTa
(52)
Note that the term (Io (())(())Tal can be pre-computed since v(l) is known.
21

---

Estimation error for decision boundary
A
C
Standard VEM
Modified VEM
1.25
leameed truc
leamedd
truc
1.00
e Worue
0.75
A
0.50
5
5
0.25
-5
X1 5
-5 X1
5
standard modified
Figure 5: Comparison between the standard vEM approach in Duncker et al. [10] and our modified vEM approach. A. Estimation error between the true and learned decision boundaries, computed as deseribed in Appendix C. For each vEM approach, we fit 5 gpSLDS models with different random initializations. The estimation errors are denoted in light blue/purple dots. The runs that display in the next two panels are denoted by a solid blue/purple dot. B. The standard vEM approach fails to learn the true decision boundary of 21 O. C. The modified vEM approach precisely learns this decision boundary.
C Empirical results for new learning objective
In this section, we empirically compare the standard vEM approach from Duncker et al. [10] to our modified vEM approach in which we learn kernel hyperparameters on a partially optimized ELBO, For this experiment, we use the same synthetic dataset from our main result in Section 4.1, We fit the gpSLDS with 5 different random initializations using both standard vEM and modified vEM. For these fits, we fix the values of and d throughout learning to ensure that the resulting models are anchored to the same latent subspace (in general, they are not guaranteed to end up in the same subspace due to rotational unidentifiabilityy. Each run was fit with 50 total vEM iterations; each iteration consisted of 15 hawaad-hackwaadd solves to update (xx) and 300 Adam gradient steps with a learning rate of 0.01 to update kernel hyperparameters.
To compare the quality of the learned hyperparameters, we quantitatively assess the error between the learned and true decision boundaries. In this simple example with 2, the decision boundary can be parametrized as wp + W + W2#2 0 for some w (wo, W), w2) The true decision boundary is characterized by Wiue = (0, 1) We denote the learned decision boundary as vb. Next, metric between the learned and true decision boundaries follows. We first normalize the learned decision boundary and define uporm
compute
sb
We do not need to do this for
Winee since it is already normalized. Then, we use the error metric
err( wù Wirne) = min @0 - Winel 1/,, 1
+ Wirue ll)
(53)
Including both terms in the minimum necessary due to unidentifiabillty of the signs of to.
Figure 5A compares this error metric across the model fits for each vEM method. is clear that the modified vEM approach consistently outperforms the standard vEM approach in terms of more accurately estimating the decision boundary. In addition, the error for standard vEM has much higher variance, since the algorithm is prone to getting stuck in local maxima of the ELBO. Figures 5B-C display the learned versus true decision boundaries in the latent space for two selected runs. We select the run from each vEM approach which achieved the lowest decision boundary error metric, as denoted by solid dots in fig. 5A. We find that the model fit with standard vEM learns a decision boundary which noticeably deviates from the true boundary. On the other hand, the model with with modified vEM recovers the true boundary almost perfectly. This illustrates that our modified approach dramatically improves kernel hyperparameter estimation in practice and enables the gPSLDD to be much more interpretable in the latent space.
22

---

true latents and dynamics
observed data
C
inferred latents
inferned
0
0.0
3.5
1.0
1.5
time in seconds
time seconds
learned partition
pasteriar var. of dynamics
10-1
R
0.5
10a
10a5
Figure 6: Additional synthetic data results on 2D limit cycle from a gpSLDS fit with quadratic decision boundaries. A. True dynamics and true latent trajectories on example trials used to generate the dataset. Dynamics are unstable rotation and stable rotation with fixed points at (0, 0), separated by of + z 4. B. Poisson process observations for an example trial. C. gpSLDS inferred latent trajectory with 95% posterior credible intervals for example trial. D. The learned T (x) accurately recovers the true circular boundary between the two sets of linear dynamics. E. The gpSLDS learned posterior variance on dynamics. The posterior variance is low in regions heavily traversed by the true latent paths, and is high in regions with little to no data.
D Additional synthetic data results
To further demonstrate the expressivity of the gpSLDS over the rSLDS, we apply the gpSLDS to synthetic dataset where the true decision boundary between linear regimes is nonlinear. The rSLDS can only model linear decision boundaries in order for its inference algorithm remain tractable.
For this example, we generate a synthetic dataset consisting of an unstable linear system and stable linear system separated by the decision boundary z + z 4. Both of the linear systems have fixed points at (0, 0). The smooth combination of these linear systems results in a 2D limit cycle (fig. 6A). We simulate 30 trials of Poisson process observations from D 50 neurons over T - 2 seconds (fig. 6B). We initialize the observation model parameters C and d using a Poisson LDS with data binned at 20ms. Then, we fit gpSLDS with J = regimes, and with π T(*) modeled using the feature transformation
Φ(x) =
x7
x3
(54)
The results of this experiment are shown in fig. 6C-E. In fig. 6C we find that the gpSLDS successfully recovers the the true latent trajectory with accurate posterior credible intervals for example trial. Furthermore, fig. 6D demonstrates that by using the quadratic feature transformation in eq. (54), the gpSLDS accurately learns the true flow field and true decision boundary =f + xx3 4. In addition, the values of π (xx) smoothly transition between 0 and near this boundary, highlighting the ability of our method to learn smooth dynamics if present. Lastly, in fig. 6E plot the inferred posterior variance of our method. We find that the gpSLDS is more confident in regions of the latent space with more data (e.g. at the decision boundary), and less confident in regions of latent space with little to no data.
E Computing resources
We fit all of our models on a NVIDIA A100 GPU on an internal computing cluster. A breakdown of approximate compute times for the main experiments in this paper includes:
23

---

Synthetic data results in Section 4.1: 1.5 hours per model fit, 440 hours for the entire experiment.
Real data results in Section 4.2: 1.5 hours per model fit, ~8 hours for the entire experiment.
Real data results in Section 4.3: hour per model fit, ~5 hours for the entire experiment,
We note that these estimates do not include the full set of experiments we performed while carrying out this project (such preliminary or failed experiments).
24

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: In our abstract and introduction, propose a new method for interpretable modeling of neural dynamics. In our paper, we apply our method and demonstrate its interpretability and competitive performance on a synthetic dataset and two real datasets.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
We discuss limitations of our work in the discussion section of the paper.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiatiinn asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
25

---

Answer: [NA]
Justification:
Our paper does not contain theoretical results.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: For experiments, provide detailed descriptions and references to datasets, as well as descriptions of the models we fit.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is new model (e-g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path reproducing verifying the results. 5. Open access to data and code
26

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification: We have made our codebase on GitHub publicly available. We provide both data and code for our main synthetic data result, and we also provide references to the dataset and instructions for fitting models to reproduce the second real data result. We do not provide data or code for the first real data result since that data has not been released to the public.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nnps.cc
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: / //iiss /tippscc/pppliu httpccpwwccpcpcueeeeeeeoeppppeee httgg/uddde//euuubnnnnnnne www.ibbiissoooooioic.. for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: We provide all details for the experimental setting that are needed to understand the results for all three sets of experiments.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: Yes, we report thorough quantitative comparisons with error bars for the experiments in Section 4.1 and 4.2. Throughout the paper, we fit multiple (i.c. at least 5) models with different initializations for experiments.
Guidelines:
The answer NA means that the paper does not include experiments.
27

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We provide information on computing resources in Appendix E.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https: / http//nnurirp.ppe httpps/ccpppplicecicccceiiiieeee
Answer: [Yes]
Justification:
We have reviewed and followed the Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: Our work develops a new statistical modeling tool for the analysis of neu- roscience data in constrained experiments. We expect that this work will not have broad socictal impacts, positive negative.
28

---

Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
This paper does not involve data or models which have high risk of misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: We cite the original papers which analyzed produced both of the real datasets that we consider in Sections 4.2 and 4.3.
Guidelines:
The answer NA means that the paper does not existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
29

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, paperswithcode com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
We provide our codebase on GitHub which comes with clear documentation.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moodl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
This paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
This paper does involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
30

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
31


Reviews:
[{"review_id": "8f9lbTSjZ3", "reviewer": "Reviewer_rveQ", "summary": "This paper introduces the Gaussian Process Switching Linear Dynamical System (gpSLDS), a novel approach for modeling latent neural dynamics. This model extends the Gaussian process stochastic differential equations framework by incorporating a new kernel function that supports smoothly interpolated locally linear dynamics. This innovation allows the gpSLDS to maintain the expressiveness needed for complex systems while enhancing interpretability, addressing the limitations of the rSLDS. The paper's contributions include the development of the gpSLDS model, introduction of a novel kernel function to balance expressiveness and interpretability, and a new learning algorithm that improves the accuracy of kernel hyperparameter estimation. The model's effectiveness is shown through applications to both synthetic and real neuroscience data, showing superior performance compared to existing methods like the rSLDS. This advancement provides a more robust framework for understanding neural computation and dynamics.", "questions": "N/A", "limitations": "The authors should explore and discuss the limitations of the gpSLDS, including specific scenarios or conditions under which the model may underperform or fail.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The gpSLDS introduces a novel kernel function within the Gaussian process framework that uniquely addresses the trade-off between model expressiveness and interpretability in the analysis of neural dynamics. This model contribution is original in the using of switching dynamics within Gaussian process, allowing for smoothly interpolated transitions between local linear regimes.\n\n2.  The methodology presented in the paper is well-developed. The paper proposes a GP-SDE model with a well-conceived kernel that facilitates a locally linear interpretation of complex dynamic behaviors.\n\n3. The paper is well-written with a clear structure that guides the reader through the problem motivation, model formulation, and experiments.", "weaknesses": "1. The proposed model effectively transforms a linear dynamical system into a Gaussian Process with a well-designed kernel. While this adaptation supports smoothly interpolated locally linear dynamics, it does compromise computational efficiency, shifting from linear time to cubic time cost.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for taking the time to read our submission and for noting the strengths of our work! We are especially glad to see the reviewer’s comments on the originality and soundness of our method, as well as the clarity of our submission. \n\nA large part of the review centers around questions of computational efficiency and limitations of the gpSLDS. We refer the reviewer to the General Response, in which we address both of those themes. In addition, we include more specific responses to your individual review here. \n\n## Weaknesses\n**Re. Modeling benefits vs. computational tradeoffs**\n> While this adaptation supports smoothly interpolated locally linear dynamics, it does compromise computational efficiency, shifting from linear time to cubic time cost.\n\nWe provide a general discussion of the computational scaling of our algorithm in General Response A, showing that our algorithm still scales linearly in sequence length. However, it is true that it incurs larger computational costs in the scaling with respect to the latent state dimensionality. \n\nWhile the gpSLDS does introduce more computational complexity, it brings several key modeling advantages: \n\n1. The gpSLDS can smoothly interpolate between locally linear dynamics. This maintains interpretability by allowing each component to be analyzed downstream using principles of linear systems, while achieving greater expressivity by being able to learn nonlinearities between these linear regimes. Crucially, this resolves problems commonly experienced in the rSLDS, such as artifactual oscillations of dynamics at regime boundaries (Fig. 2G). \n\n2. Our GP-based approach allows us to infer approximate posterior distributions over dynamics at any point in the latent space, whereas the rSLDS often infers uninterpretable dynamics at regime boundaries and does not explicitly treat dynamics parameters as probabilistic quantities. \n\nWe believe that these modeling advantages of the gpSLDS effectively address limitations in the rSLDS that hinder its interpretability in practice. In addition, our implementation of the gpSLDS and other GP-SDE based models leverages fast parallelization and automatic differentiation in JAX with GPU compatibility. This allows the gpSLDS to achieve more efficient runtimes on a per-timestep basis than the most widely-used rSLDS implementation (see General Response A & attached PDF), despite the additional complexity in GP inference. We will release a codebase with the paper if accepted, which we believe will provide a valuable and practical tool for practitioners.\n\n## Limitations\n> The authors should explore and discuss the limitations of the gpSLDS, including specific scenarios or conditions under which the model may underperform or fail.\n\nWe refer the reviewer to General Response B for a discussion on the limitations of the gpSLDS, as well as possible model extensions.\n\n**Thank you again for your positive comments and thoughtful response!** If you have any further questions, we would be happy to answer them.\n\n– The gpSLDS authors", "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": "", "extra_scores": {}}, {"review_id": "WJS6jPVGAb", "reviewer": "Reviewer_9ds8", "summary": "The paper explores latent state inference and parameter learning within a switching stochastic dynamical system. In this context, the dynamics are represented by a stochastic differential equation, with the drift function modeled as a Gaussian process. Notably, the paper introduces a novel kernel for this Gaussian process—a mixture of linear functions that captures the system’s switching behavior. The inference process employs a variational inference framework. To validate its effectiveness, the paper conducts evaluations using both synthetic and real neuroscience data.", "questions": "- Could the authors provide further details on how the presented approach could incorporate a multi-dimensional correlated Gaussian process prior?", "limitations": "-", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "The proposed method offers a fresh perspective on inference for switching dynamical systems. It introduces a novel Gaussian process (GP) kernel that captures the switching behavior between linear functions. The experimental results are convincing and fair, and the paper is well-written and easy to follow.\n\nAdditionally, the authors suggest a minor modification to the variational inference (VI) algorithm to achieve faster updates.", "weaknesses": "My primary concern with this work lies in how the latent dynamics, denoted as $f$, are modeled as independent draws from the same Gaussian process (GP) prior. In my view, this implies that the posterior over latent components should be independent in each component. Consequently, learning becomes challenging when only single components of the system are observed. This limitation could pose a substantial problem, especially in classical tracking scenarios where only specific state components are observable.\n\nAdditionally, I perceive the contribution as somewhat marginal. Essentially, the paper introduces a new GP kernel, but beyond that, the impact seems limited.\n\nChoice of Linear Kernel:\n\nI wonder about the rationale behind using a linear kernel. Why not explore switching between nonlinear kernels? \n\nInference Algorithm Analysis: \n\nLastly, a detailed examination of the variational inference algorithm would provide valuable insights. Additionally, a synthetic example for comparison could help demonstrate the algorithm’s performance.\n\nFor instance, the authors could select a linear stochastic differential equation (SDE) and perform closed-form parameter estimation using methods like the Expectation Maximization (EM) algorithm. In this scenario, latent state inference could be achieved through RTS smoothing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for taking the time to read our submission and for providing insightful feedback. We especially appreciated that the reviewer found our submission to be clearly written and a \"fresh perspective\" on SLDS models!\n\n## Weaknesses\n**Re. Prior independence assumption**\\\nWe thank the reviewer for this comment and agree it is important to consider how certain modeling choices might introduce estimation error or limit expressivity. Here, we provide a discussion on this point which we'll also include in the paper. We note that in our setting, we assume high-dimensional observations are driven by mixtures of latent components. Observing a lower-dimensional projection of the SDE where only single components are observed is not a setting we consider here, but would be an interesting extension.\n\nWe model $f$ using independent GPs, following a body of previous work [Eleftheriadis et al. 2017; Duncker et al. 2019; Fan et al. 2023]. We note that although the prior assumes independence, **this does not imply that the true posterior over $f$ is independent across dimensions**, as the likelihood depends on the latent state $x$ which combines $f$ across dimensions. We approximate the true posterior using a variational approximation that factorizes over latent dimensions. This enables a tractable inference algorithm which propagates posterior uncertainty to the inference and prediction of latent dynamics. While we could add covariance terms to the variational approximation, this would introduce more parameters which may complicate inference and learning. Nonetheless, we agree it is important to carefully assess biases that the variational approximation may induce in the recovered estimates (e.g. Turner and Sahani 2011) and this will be a topic of future work.\n\n**Re. Impact of contribution**\n> the paper introduces a new GP kernel, but beyond that, the impact seems limited\n\nWhile we do introduce a new GP kernel, we believe its broader implications are significant for both the ML and neuroscience communities. We identified key limitations of the rSLDS and drew a nontrivial connection to GP-SDEs via the design of a novel kernel, retaining many of the advantages of the rSLDS while addressing its limitations. The gpSLDS extends a rich line of work on rSLDS and related models which have made a significant impact on interpretable data analysis in neuroscience [Taghia et al. 2018; Costa et al. 2019; Nair et al. 2023; Liu et al. 2023, Vinograd et al. 2024]. We have also developed a fast and efficient JAX codebase implementing the gpSLDS and other GP-SDE models. We will release our codebase to the public upon acceptance, providing practitioners with a valuable tool for modeling neural dynamics.\n\n**Re. Why linear kernel?**\n> I wonder about the rationale behind using a linear kernel. Why not explore switching between nonlinear kernels?\n\nWe designed the gpSLDS to switch between linear kernels to impose interpretable structure on complex nonlinear dynamics so that they can be easily analyzed downstream. Typical analyses of nonlinear dynamics focus on linearized dynamics around fixed points [Duncker et al. 2019; Sussillo & Barak 2013], and often require second-stage analyses like fixed-point finding [Golub & Sussillo 2018; Smith et al. 2021]. In contrast, interpretable features are readily available in the gpSLDS due to its piecewise-linear structure. In neuroscience, dynamical motifs of linear systems are hypothesized to underlie various kinds of neural computations (e.g. line attractors for evidence integration, rotational dynamics for motor control) [Vyas et al. 2020]. Our goal is to extract these features from neural data in an interpretable way. \n\nWhile it is straightforward to extend our GP kernel to switch between nonlinear kernels instead, this added model flexibility would likely make it difficult to correctly learn regime boundaries. In principle, a sufficiently expressive nonlinear kernel may not need to switch at all. In our case, by allowing linear functions to switch as a learnable function of the latent state, we can capture complex nonlinearities in the dynamics while also adding structure to make parameter learning more feasible.\n\n**Re. Validating the variational inference alg.**\n> a synthetic example for comparison could help demonstrate the algorithm’s performance\n\nWe note that our experiment in Appendix C shows the improved performance of our modified vEM algorithm over standard vEM on synthetic data. For this experiment, we used the same synthetic dataset as in Figure 2 of the main text (two linear systems separated by a vertical boundary). \n> The authors could select a linear SDE and perform closed-form parameter estimation...\n\nThanks for this suggestion – while this would indeed allow us to compare our vEM algorithm to closed-form EM, our focus is on performing approximate posterior inference for nonlinear SDEs, for which closed-form updates are not available. Therefore, we performed a slightly different experiment in Appendix C: comparing our modified vEM algorithm to standard vEM on a dataset with simple, but nonlinear, dynamics. We chose this setting because the nonlinearity introduces complex dependencies between the dynamics and kernel hyperparameters that would not be present in the linear case.\n\n## Questions\n> …incorporate a multi-dimensional correlated Gaussian process prior?\n\nAs we discuss above, the prior independence assumption does not imply that the true posterior is independent across dimensions. In principle it is possible to incorporate correlation structure into the prior, but this would introduce $O(K^2)$ more model parameters that may be harder to learn. Since the true posterior is still correlated across dimensions even with an independent GP prior, we decided to stick with the simpler approach. \n\n**Thank you again for your positive comments and insightful response!** If you have any further questions, we would be happy to answer them.\n\n– The gpSLDS authors", "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": "", "extra_scores": {}}, {"review_id": "4w7gNj6Oey", "reviewer": "Reviewer_E7YR", "summary": "This paper proposes a new model called Gaussian Process Switching Linear Dynamical System (gpSLDS). The model is more interpretable and infers more stable latent compared with the alternative rSLDS. Particularly, the weird oscillations of the latent can be avoided due to the newly proposed Smoothly Switching Linear (SSL) kernel. Extensive experimental results on one synthetic and two real-world datasets shows the benefit of gpSLDS compared with the one with the traditional RBF kernel and the rSLDS.", "questions": "* Typo: Line 129, function's slope of\n* Figure 1A, if different colors are different samples (should be clearly state if so), then it is meaningless to show these samples unless kernel parameters are told.\n* Line 177, $\\boldsymbol z_m \\in \\mathbb{R}^D$? And it seems like this $\\boldsymbol z$ is different from the $z_j$.\n* What does the conjugacy of the model mean in line 193?\n* How about the effect of choosing $K$ and $J$. At least some explanations or guidance and the corresponding effects are necessary for them.\n* How about the algorithm complexity. Computing GP is time consuming, which might hinder the computing efficiency for long sequences.", "limitations": "/", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "* This proposed method is detailed in math and intuitive with some instructive explanations.\n* One synthetic and two real-world experiments are done, which is good.", "weaknesses": "* Some of the presentations can be improved.\n* Lacks some complexity analysis or results. This might be an important factor to be considered for long sequences.\n\nSee questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for taking the time to read our submission and for providing helpful feedback! We are especially pleased that the reviewer found our methodology to be intuitive and supported by solid experimental results.\n\n## Weaknesses\n**Re. Complexity analysis**\n> Lacks some complexity analysis or results. This might be an important factor to be considered for long sequences.\n\nWe refer the reviewer to General Response A for an analysis and discussion of the computational complexity of the gpSLDS. We highlight that while typical GP inference methods scale as $O(T^3)$ where $T$ is the number of time steps, we overcome this by employing inducing points as in Titsias (2009) and Duncker et al. (2019) to reduce this complexity to $O(TM^2)$ where $M$ is the number of inducing points.  It is true, however, that the gpSLDS incurs larger computational costs in scaling with respect to latent dimension due to using quadrature to approximate kernel expectations. For many real world datasets, such as the ones we explore in our paper, the key dynamical features of interest can be captured in low latent dimensions. That said, there will be some applications which require large latent dimensionalities; for these cases, one possible workaround would be to use Monte Carlo methods instead of quadrature. We also highlight that in practice, our implementation of the gpSLDS leverages fast parallelization and auto-differentiation in JAX, which yields faster runtimes per time bin than the rSLDS (see General Response A & attached PDF).\n\n## Questions\n**Re. Line 129**\n\nThank you for catching this typo! We will be sure to fix this in our paper.\n\n**Re. Figure 1A**\n\nIn Figure 1A, different colors do indeed represent different samples from the GP. We will make sure to clearly state this in the figure caption. We decided to show multiple samples for each GP kernel to provide intuition for what each distribution over functions looks like and how they build up to the SSL kernel. \n\n**Re. Line 177 notation**\n\nThanks for catching this typo as well! Line 177 should say $\\mathbf{z}_m \\in \\mathbb{R}^K$. Thank you also for pointing out the overloaded notation for $z$. We will keep the notation for the inducing inputs and will change the notation for the rSLDS discrete states to $s_j$. \n\n**Re. Model conjugacy**\n\nHere, model conjugacy refers to the fact that we can compute the distributions $q(\\mathbf{u}_k)$ which maximizes the ELBO in closed-form. This is due to the conjugacy between the Gaussian prior on inducing points $p(\\mathbf{u}_k \\mid \\Theta) = \\mathcal{N}(\\mathbf{u}\\_k \\mid 0, \\mathbf{K}\\_{zz})$ and its corresponding Gaussian variational posterior $q(\\mathbf{u}\\_k) = \\mathcal{N}(\\mathbf{u}\\_k \\mid \\mathbf{m}\\_u^{k}, \\mathbf{S}_u^{k})$. We provide more detail about this update step in Appendix B.4, and present the closed-form updates for $\\mathbf{m}\\_u^{k}$ and $\\mathbf{S}_u^{k}$ in Equations (40)-(41) of that section.\n\n**Re. Choosing $K$ and $J$**\n\nThanks for this question. In cases where we do not know the true latent dimensionality or number of linear regimes, we can use standard model comparison metrics in the neural latent variable modeling literature to choose these hyperparameters. For example, we could compare models with different $K$ and $J$ based on forward simulation accuracy, which measures how well we can predict neural activity if we sample future latent states from the fitted model [Nassar et al. 2019; Nair et al. 2023]. Another option would be to use a validation technique called co-smoothing, which involves re-running the inference step of a fitted model on held-out trials after withholding some neurons, and then evaluating the expected log-likelihood on those withheld neurons [Macke et al. 2011; Wu et al. 2018; Keeley et al. 2020]. We will add a discussion on choosing $K$ and $J$ to our paper.\n\n**Thank you again for your positive comments and insightful response!** If you have any further questions, we would be happy to answer them.\n\n– The gpSLDS authors", "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": "", "extra_scores": {}}, {"review_id": "PO1wOaSkQv", "reviewer": "Reviewer_3Mtv", "summary": "This paper introduces a new approach to model the low-dimensional latent dynamics of a collection of neurons over time. Their approach balances the two desiderate of capturing complex nonlinear dynamics while remaining interpretable. Specifically, they introduce the Gaussian Process Switching Linear Dynamical System, which models the latent state as a Gaussian process with a specific novel kernel function, which interpolates between different linear dynamics akin to a switching linear dynamical system. This mitigates some of the pathologies of rSLDS and also provides uncertainty estimates for the dynamics. The authors validate their method both on synthetic and experimental data from neural recordings.", "questions": "- Figure 3: GP with RBF kernel seems to perform better with respect to forward simulation, why do you think that is?\n- Could you comment on how expensive the different approaches are relative to each other in terms of training and prediction?", "limitations": "While the paper discusses limitations of their method throughout the paper, this is a bit sparse. I think the paper could be improved by a dedicated limitations paragraph in the discussion in Section 6.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "This paper tackles an important topic in neuroscience and provides interesting methodological advances continuing a line of research on inferring the latent dynamics of neuronal populations. \nThe paper is overall well-written and provides a concise and clear exposition.\nFor example, Figure 2G very clearly illustrates the problems that recurrent switching linear dynamical systems have.\nThe methodology seems sound to me and is motivated clearly in terms of the issues with rSLDS.\nThe experiments are well-executed and illustrative and importantly also use real data from two animal experiments tackling questions about behavior and decision-making. In addition, the authors provide code for their method and for two of the experiments.", "weaknesses": "The biggest weakness the paper currently has in my opinion is that it is hard to judge what the computational trade-offs are between the different methods. While the authors provide estimates for the overall compute they use in all experiments (see Appendix E), there is no comparison between different methods. I would have liked to see an analysis that compares GP-RBF, rSLDS and gpSLDS not only as a function of the number of trials or number of steps simulated forward, but also in terms of their predictive performance as a function of training time / budget. Otherwise, it is hard to judge whether these methods were compared fairly.\n\nOne aspect that would make the paper stronger, but is not strictly necessary in my opinion, is to add another baseline, which is different from the two approaches that were combined in gpSLDS, e.g. something like LFADS (Sussillo et al, 2016).\n\n- Sussillo, D., Jozefowicz, R., Abbott, L., & Pandarinath, C. (2016). LFADS - Latent Factor Analysis via Dynamical Systems. In Advances in Neural Information Processing Systems (NeurIPS)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for taking the time to read our submission and for noting the strengths of our work! We are especially pleased that the reviewer praised our work's clear motivation, solid experimental results, and relevance to the neuroscience community.\n\n## Weaknesses\n**Re. Runtime comparison**\n\nThank you for this suggestion – we refer the reviewer to General Response A for new experimental results and a discussion about comparing runtimes and accuracy between the three methods. To further address your specific concerns:\n\n> predictive performance as a function of training time / budget\n\nInstead of comparing the predictive performance as a function of training time, we decided to compare the accuracy of the recovered latent states and dynamics and corresponding runtimes for a fixed number of vEM iterations sufficient for convergence. We did this for the following reasons:\n\n1. Measuring how well the model recovers ground truth latent states and dynamics is the most direct way of evaluating the model’s performance on synthetic data, as we are primarily interested in how well the model captures the underlying latent variables.\n\n2. For each model class, our aim is to fit the best possible model (in terms of ground truth recovery) using our available computing resources. Therefore, we first fit all models to convergence and then compared their runtimes and performance. We believe that this approach best reflects how these models are applied in practice. While investigating the learning dynamics as a function of fitting time is certainly an interesting question, it is somewhat orthogonal to the primary scientific goal of understanding latent neural dynamics. \n\n> hard to judge whether these methods were compared fairly\n\nSince we fit methods using different discretization steps to ensure the numerical stability in each model class, we report total runtime as well as runtime normalized by number of time bins. This allows us to fairly compare runtimes, which we expect to scale linearly with time bins in all 3 models. \n\nFinally, we note that training time depends on many factors specific to each individual implementation, and this should be taken into account when comparing runtimes. In General Response A & the attached PDF, we show that our implementation of the gpSLDS and other GP-SDE models – which leverages fast and efficient parallelization and auto differentiation in JAX – is more efficient per time bin than the most widely-used rSLDS implementation. We believe that this is a valuable contribution to the machine learning and neuroscience communities in itself, and we plan to release a codebase to the public along with the paper.\n\n**Re. LFADS**\n\nThank you for the suggestion to compare our method with LFADS. We agree it would be interesting to explore this comparison in future work. On a high level, we expect LFADS to reconstruct neural data well, especially given sufficient data. However, the RNN dynamics of LFADS can be difficult to interpret after fitting and do not come with uncertainty estimates. Moreover, because LFADS is a deep learning model with many more parameters than the gpSLDS, we anticipate that it would struggle more to learn accurate dynamics in data-scarce settings, such as the hypothalamic dataset we analyze in our paper. In contrast, models with structured probabilistic priors, such as the gpSLDS, are better suited to correctly infer key dynamical motifs in these settings.\n\n## Questions\n**Re. Comparison to GP-SDE w/ RBF kernel**\n> Figure 3: GP with RBF kernel seems to perform better with respect to forward simulation, why do you think that is?\n\nThis is likely because the GP-SDE with RBF kernel is more expressive than the gpSLDS kernel, but it comes at the cost of being less interpretable. The RBF kernel is universal in the sense that it can approximate any smooth function with arbitrarily small error [Micchelli and Xu 2006]. However, because of its flexibility, key dynamical features such as fixed points are not readily available for downstream analysis. On the other hand, the gpSLDS finds the best (smoothly-interpolating) piecewise-linear approximation to a nonlinear system. This aids in interpretability by showing how nonlinear dynamics can be partitioned into linear regimes, each of which can be individually analyzed. Therefore, we expect that the gpSLDS finds dynamics which are less accurate but more representative of dynamical motifs of interest.\n\n**Re. Computational complexity and tradeoffs**\n> Could you comment on how expensive the different approaches are relative to each other in terms of training and prediction?\n\nWe provide a discussion on the computational complexity of the gpSLDS in General Response A. All 3 methods use inference algorithms that scale linearly in sequence length. With respect to latent dimension, the gpSLDS scales exponentially (due to using quadrature) while the rSLDS scales cubically. However, the exponential scaling in the gpSLDS could be overcome by using Monte Carlo approximations instead of quadrature for large latent dimensionalities. \n\nIn addition, all methods can predict dynamics at new locations efficiently once the models are fit. For the GP-SDE based models, predicting dynamics at a batch of $B$ new locations costs $O(KBM^2)$, where $M$ is the number of sparse inducing points. For the rSLDS, this operation costs $O(KB)$ by reading off the dynamics of the most likely discrete state. In practice, the main cost comes from fitting the model rather than from predicting dynamics. We will include a discussion about relative computational tradeoffs to our paper.\n\n**Re. Addressing limitations**\n\nThanks for this suggestion – we have written a discussion on limitations and possible model extensions in General Response B, which we plan to add as its own paragraph in the paper.\n\n**Thank you again for your positive comments and insightful response!** If you have any further questions, we would be happy to answer them.\n\n– The gpSLDS authors", "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": "", "extra_scores": {}}]

Author rebuttal:
We thank the reviewers for taking the time to read our submission and for providing thoughtful and insightful feedback! We were pleased that the reviewers unanimously supported our submission as a valuable contribution to the NeurIPS community, citing that it was **1) easy to follow, 2) clearly motivated, and 3) backed by convincing experimental results.**

Here, we address two main themes that were brought up in several of the reviews. We respond to individual reviewer concerns separately.

# (A) Computational complexity of gpSLDS
Reviewers rveQ, E7YR, and 3Mtv brought up questions about the computational complexity and efficiency of the gpSLDS. 

Performing inference and hyperparameter learning in the gpSLDS relies on computing expectations of $f(\cdot)$ with respect to the variational marginals $q(x(t)) \sim \mathcal{N}(m(t), S(t))$ and the approximate posterior GP $q(f)$. GP-based methods are computationally expensive for two reasons: 1) they typically scale cubically in the number of input points to the GP, and 2) evaluating posterior expectations with respect to distributions over GP inputs involve computing expectations of nonlinear kernel functions, which are typically not available in closed form. 

To overcome 1), we follow Titsias (2009) and Duncker et al. (2019) in using inducing points to perform inference of $q(f)$. This reduces the computational complexity of evaluating $f(\cdot)$ on a sequence of $T$ latent states from $O(KT^3)$ to $O(KTM^2)$ for $M$ inducing points (in our synthetic experiments we choose $M=16$ ($4 \times 4$ grid) for a latent dimensionality of $K=2$). For 2) we perform quadrature with $N$ nodes per latent dimension, so the total number of nodes needed to accurately approximate kernel expectations scales as $O(N^K)$ (for all experiments we use $N=6$). This represents the main computational bottleneck of the algorithm and can pose challenges for fitting the gpSLDS in settings with large latent dimensionality. In the real datasets we study in our paper, a small number of latent dimensions could sufficiently capture key dynamical features of interest, and we expect this to be the case in other applications as well. However, in cases that require larger latent states, it is possible to use Monte Carlo methods instead of quadrature to approximate kernel expectations. We will include a thorough discussion of computational complexity and possible model extensions to our paper.

Reviewer 3Mtv suggested doing a runtime comparison between gpSLDS, GP-SDE with RBF kernel, and rSLDS. **In the attached PDF**, we include additional experimental results comparing the runtime and accuracy for the three methods on the synthetic dataset from Fig. 2. For the gpSLDS and GP-SDE with RBF kernel, we choose a discretization step of 1ms, yielding $T = 2500$ time bins. For the rSLDS, we choose a discretization step of 20ms, yielding $T = 125$ time bins. We observed that using too large of a discretization step for the GP-SDE based methods and using too small of a discretization step for the rSLDS can in some cases lead to numerical instabilities, due to their respective continuous-time and discrete-time model formulations. Therefore, we chose these settings in order to fairly compare the three models. We ran each method with $K = 2$ for $100$ variational EM iterations, which allowed all three methods’ ELBOs to converge. To deal with the differing sequence lengths, we report both total runtime and runtime normalized by the sequence length. We report standard errors across 5 runs with different random initializations per model.

These results demonstrate a computational tradeoff between the GP-SDE based models and the rSLDS: while the GP-SDE based models require more time bins (via smaller discretization steps) to accurately approximate continuous-time dynamics, they also recover the true latent variables with a much higher degree of accuracy on this continuous-time point process dataset than the rSLDS. In addition, these results show that our own implementations of the gpSLDS and GP-SDE with RBF kernel are more efficient than the most widely-used rSLDS implementation on a per time-bin basis. That being said, we note that differences in runtime are impacted by differences in implementation, such as the discretization step size, GPU compatibility, and choice of optimizer, which vary across models. **We stress that the focus of our contribution is not to optimize for runtimes, but rather to find the best-fitting model possible in order to draw accurate and reliable scientific conclusions.**

# (B) Clearly addressing limitations
Reviewers rveQ and 3Mtv suggested more clearly presenting the limitations of the gpSLDS. We agree that this would be a valuable addition to the paper, and we will include a paragraph discussing the following points:

1. Quadrature and scalability: As mentioned above, because we use quadrature methods to compute kernel expectations that are not available in closed-form, it could be challenging to scale the gpSLDS to settings with large latent dimensionalities. One possible extension could be to explore using Monte Carlo methods to approximate these expectations in higher latent dimensions.

2. Inference algorithm alternatives: One area of potential improvement for the gpSLDS would be to incorporate more recent methods for latent state inference in GP-SDEs [Verma et al. (2024); Course and Nair (2024)]. In particular, Verma et al. (2024) proposed an algorithm for inferring $q(x)$ akin to gradient descent, which they show achieves faster and more stable convergence compared to the fixed-point iteration method we use in our paper [Archambeau et al. (2006)]. While Verma et al. did not originally consider inference over $q(f)$, their method could be directly plugged into the inference algorithm for the gpSLDS.

**Again, we thank all of the reviewers for their positive comments and constructive feedback!** 

– The gpSLDS authors

Meta-review: By taking a product of a linear kernel and a "partition" kernel, gpSLDS elegantly improves Duncker & Sahani's work on estimating locally linear latent dynamical systems. The paper is a key advance in the topic, and it is clearly written with detailed experiments.
Decision: Accept (poster)