Probabilistic Decomposed Linear Dynamical Systems for Robust Discovery of Latent Neural Dynamics
Authors: Yenho Chen
Noga Mudrik
Kyle A. Johnsen
Sankaraleengam Alagapan
Adam Shabti Charles
Christopher John Rozell
Venue: NeurIPS 2024 poster (2024)
Keywords: Computational Neuroscience
Probabilistic Modeling
State Space Models
Dynamical Systems
Primary area: neuroscience_and_cognitive_science

Abstract: Time-varying linear state-space models are powerful tools for obtaining mathematically interpretable representations of neural signals. For example, switching and decomposed models describe complex systems using latent variables that evolve according to simple locally linear dynamics. However, existing methods for latent variable estimation are not robust to dynamical noise and system nonlinearity due to noise-sensitive inference procedures and limited model formulations. This can lead to inconsistent results on signals with similar dynamics, limiting the model's ability to provide scientific insight. In this work, we address these limitations and propose a probabilistic approach to latent variable estimation in decomposed models that improves robustness against dynamical noise. Additionally, we introduce an extended latent dynamics model to improve robustness against system nonlinearities. We evaluate our approach on several synthetic dynamical systems, including an empirically-derived brain-computer interface experiment, and demonstrate more accurate latent variable inference in nonlinear systems with diverse noise conditions. Furthermore, we apply our method to a real-world clinical neurophysiology dataset, illustrating the ability to identify interpretable and coherent structure where previous models cannot.

Probabilistic Decomposed Linear Dynamical Systems for Robust Discovery of Latent Neural Dynamics
Yenho Chen Noga Mudrik', Kyle A. Johnsen' Sankaraleengam Alagapan?, Adam S. Charles', and Christopher J. Rozell
Machine Learning Center, Georgia Institute of Technology
2Schoo of Electrical and Computer Engineering, Georgia Institute of Technology
SCoulter Dept. of Biomedical Engincering, Emory University and Georgia Institute of Technology
'Department of Biomedical Engineering, Mathematical Institute for Data Science, Center for Imaging Science, Kavli Neuroscience Discovery Institute, Johns Hopkins University yenhoggatech, edu, nmudrik16jhu, edu, kjohnsen@gatech.c edu sankar. allggaaanggatech. edu, jhu edu, crozell@gatech.co edu
Abstract
Time-varying linear state-space models are powerful tools for obtaining mathemat- ically interpretable representations of neural signals. For example, switching and decomposed models describe complex systems using latent variables that evolve according to simple locally linear dynamics. However, existing methods for latent variable estimation are not robust dynamical noise and system nonlinearity due to noise-sensitive inference procedures and limited model formulations. This lead to inconsistent results on signals with similar dynamics, limiting the model's ' ability provide scientific insight. In this work, we address these limitations and propose a probabilistic approach to latent variable estimation decomposed models that improves robustness against dynamical noise. Additionally, introduce tended latent dynamics model to improve robustness against system nonlinearities. We evaluate our approach several synthetic dynamical systems, including empirically-deriin brain-computer interface experiment, and demonstrate more accurate latent variable inference in nonlinear systems with diverse noise condi- tions. Furthermore, we apply our method to real-world clinical neurophysiology dataset, illustrating the ability identify interpretable and coherent structure where previous models cannot.
Introduction
central goal computational neuroscience is to develop models capable of discovering latent structure within noisy, high-dimensional neural signals. By identifying hidden relationships within neural recordings, we can begin to understand, predict, and control the behaviors of the underlying systems. Modeling neural time-series is challenging due to the range of temporal dynamics present. For example, there may be gradual short-term fluctuations, abrupt shifts in response to external stimuli,
and long-term global drifts resulting from changes in baseline activity levels (36 43 39
Although black-box approaches based on deep learning are available 33 20 their complexity
often obscures the relationships learned from the data, making it difficult to extract scientific insights from these models. As result, practitioners may favor time-varying linear state-space models which offer mathematicalyy interpretable representations by approximating complex dynamics with simple
"Code is available at:
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

locally linear regimes [6]. However, obtaining latent variable estimates that are robust to dynamical noise and system nonlinearity in these state-space models is challenging. When applied to neural time-series, latent variable estimation may become unstable due to inflexible model formulations or sensi-sen inference procedures. This can incorrectly produce disparate results for signals generated from the same underlying system.
For example, switching linear dynamical systems (SLDS) and related models 42 44 segment time-series into discrete linear dynamical states, providing a piecewise linear approximation of the underlying system while highlighting coherent groups of activity. However, the assumption of discrete components can be a poor modeling choice for neural signals that contain continuous-vall fluctuations, such as gradual or random changes to the system speed as seen during neural ramping activity [I or Lévy walk dynamics in the cerebral cortex [27]. As demonstrated in [30] and our
experiments in Section 4.3] when applied to real- datasets, inference of the switching variables can result in rapid, random oscillations between the discrete modes, indicating that the model is unable to identify meaningful structure in the data.
To address the limitations of discrete states, the decomposed linear dynamical systems (dLDS) model 30) learns dictionary set of línear dynamical regimes, referred to as dynamic operators (DOs), that can be modified and combined through linear combination of sparse coefficients. By allowing coefficients be time-varying and continuuus-valued dLDS naturally captures both gradual changes by adjusting the coefficient magnitudes and abrupt shifts by changing the set of active DOs over time. Inference accomplished by optimizing over cost function that encourages data reconstruction while also constraining the structure of the dynamic coefficients to be sparse and temporally smooth.
Unfortunately, there are two critical shortcomings that prevent the robust inference of the latent variables in dLDS. First, the cost-based inference procedure sensitive to noise, because of the regularization term encouraging temporal smoothness. This term sequentially propagates errors from noisy coefficient estimates over the length of the time series. Consequently, the model may produce inconsistent coefficient estimates on similar signals and have poor multi-step inference performance, indicating that the learned dynamics are unable to generalize well beyond single time step. Second, the original latent dynamics model lacks method for accurately representing systems with multiple fixed points, causing DO coefficients to oscillate or switch arbitrarily in a way that may not align with the underlying process. To learn an effective decomposed model in practice, we require strategy that provides robust estimates of latent variables despite the presence of noise and system nonlinearity.
In this work, we address these limitations by introducing the probabilistic decomposed linear dy- namical systems (p-dLDS) model. Our approach improves robustness of latent variable estimation while maintaining the richness of a decomposed dynamical systems model. First, we propose a probabilistic inference procedure that reduces the model's sensitivity to temporal noise by accounting for uncertainty in the latent variable estimates over time. Namely, we introduce time-informed hierar- chical variables that encourage both sparse and smooth model coefficients. We devise a variational expectation maximization (vEM) procedure to perform inference and learning over this probabilistic structure. Second, we incorporate time varying offset term to model systems that orbit multiple fixed points. While we analytically identify model degeneracies with this formulation, we propose an additive decomposition strategy that prevents convergence to trivial solutions.
Through several synthetic examples, we demonstrate how these contributions lead to improved accu- racy and robustness of latent variable estimation despite difficult noise conditions. We extend these results to an empirically- -eriive brain-compute interface experiment [8] showcasing robustness to highly nonlinear observation functions and the ability to extract meaningful insights from the learned latent variables. Finally. we illustrate how our method effectively identifies interpretable and coherent structure in clinical neurophysiology dataset where previous models are unsuccessful.
Background and Related Work
State Space Models. Our goal is to accurately describe the evolution of high-dimensional time-series data ye E M with the following state-space equations,
y: = Da. + d + Eye cy ~ No Ey), (Observations)
(1)
x, = f( 1) + Cza, Cz. ~ No0, Σ,), (Dynamics)

---

where 2y E R is the latent state, fa(-) is the dynamics function, and D E R and d E RM describe linear observation function. Our work focuses on the case when N < M, which compresses high-dimensional signals into a low-dimensional latent space. By choosing fr to be time-varying linear operator, we can approximate complex nonlinear dynamics with simple locally- linear components, balancing expressivity with mathematical interpretabilityy. However, learning time-varying linear operator from data can be challenging, and typically requires additional constraints on the underlying generative model to identify meaningful representations.
Switching Linear Dynamical System (SLDS), SLDS approximates nonlinear systems by introducing a discrete switching variable 2 = (1..... } into the time-varying linear dynamics equation,
xt++ 2; + F2, ba, + Ex.-
At each time step t. the latent state D2 evolves according to the 2,-th linear regime defined by Fz, NNN and b., RNN while the switching variables evolve according to Markov matrix. Inference is performed through a vEM algorithm, where the approximate posterior of the latent variables is estimated through coordinate ascent updates over tractable subgraphs. There many extensions of SLDS, such as rSLDS [26] which modifies its generative behavior by informing the transitions of 2, with -t-1- 1. However, switching models are inherently limited when describing complex signals due to their discrete formulation. For instance, a switched representation is unable to learn that dynamic regime may exhibit a range of variations. neural systems, these variations may arise from random spiking processes [37 [4]] or systems with randomly distributed speeds In SLDS, each variation is learned as a separate discrete state, thus obscuring that the learned states are related. Furthermore, the switching formulation cannot adapt the learned system to unseen variations (i.e. different levels of random speeds). This can produce unstable inference behavior, where the switching state oscillates unpredictably or collapses to single uninformative state.
Decomposed Linear Dynamical Systems (dLDS). dLDS 30 relaxes the discrete formulation by approximating nonlinear and nonstationary signals with vime-varying mixture of linear dynamical systems (LDS) defined by the following equations and constraints,
L(+1 = It + Fiz; + Cz.:
F = E
k==1 fRCIK
s.t. Ct is sparse.
Every transition F; is decomposed as a linear combination of sparse coefficients c₁ E RK and a dictionary of K DOs fs E R R Figure shows the corresponding graphical model. Inference of the latent variables is accomplished by solving the Basis Pursuit Denoising with Dynamic Filtering (BPDN-DF) (6) objective sequentially for all and An, \1, X2 O,
Le. C. = arg min| ye Dæ, 112 + *ollx, - :.-1 - F 112 Aa/eeell + >alle, 2-1112
This produces a point estimate of x; and e, that matches the likelihood function resulting from Equations (1) and (2). In this objective, the dynamic coefficients are encouraged to be sparse through the L penalty and temporally smooth through the & penalty centered around the previous coefficient estimate. However, this approach is sensitive to noise because inference relies on propagating noisy point estimates of 2 over time. As a result, BPDN-DF may accumulate errors that can lead to significantly different coefficient estimates Furthermore, the lack of robustness to noise can degrade multi-step inference performance, causing the inferred system to quickly diverge from the true system. This suggests that the inferred latent variables only capture the local activity narrowly and are unable to accurately represent the dynamics beyond a single time-step. Another drawback of dLDS arises from the dynamics model in equation ( which implicitly assumes that the observed dynamics contain a single fixed point that revolves around the origin. This limits dLDS's ability to model systems that cannot be easily mean-centered such as those with multiple fixed points or nonstationary drifts.
signals sampled from the same generative process.
Sparse Bayesian Learning with Dynamic Filtering. Sparsity is achieved in probabilistic models through hierarchical scale-mixture priors [00] 2 4). To integrate dynamical information into prob- abilistic sparse signal inference, previous work [3]] proposes the Sparse Bayesian Learning with Dynamic Filtering (SBL-DF) framework where the following hierarchical model is defined,
p(x,, Cr, 7)) p(x)(c:) II
b..k). b.
(3)
3

---

A)
B)
Coefficient
Conficcintss
Dynamic Coefficients
Latent State
Fe-3
decomposition
Observed
Lanent
WY
8k-3
Data
Figure 1: (A) Graphical model of dLDS. (B) p-dLDS includes hierarchical variables for probabilistic sparse inference and reparameterizes the latent space include time-varying offset term.
Here, v(x,(c,) = ( >₂I) specifies the likelihood, where + is measurement matrix. Sparsity encouraged through the zero-mean Gaussian priors (t. (0, (t.)) independently placed on each element of the sparse vector. The variance parameters "/, are defined by an inverse gamma hyperprior v(Tt,x/ot,as (e,x)) Ig (or,k: (t,x)) with shape parameters 0r,k and 8r, When marginalizing over "r,,''s ss see that p(ct, dett,e 3r,x) becomes t-distribution known for its high kurtosis, which is essential for producing sparse solutions. To propogate dynamics information, the past estimate čt-1 informs the hyperprior parameters of the next estimate such that bt,k/at, ot 1
Probabilistic Decomposed Linear Dynamical Systems
We build upon dLDS and propose a probabilistic decomposed linear dynamical systems (p-dLDS) model. Rather than propagating noisy point estimates of the latent variables during inference, we improve robustness by marginalizing over uncertainty with respect to time. Additionally, we propose tractable method for extending the dynamics model to systems with multiple fixed points.
3.1 Time-narying offset term
Note that for any parameter setting, the dLDS local dynamics (eq. (2) reduces to linear dynamical system (LDS) that characterized by single fixed point centered around the origin. Yet, dynamical systems often consist of much more complicated behaviors. Nonlinearities can cause signal navigate through multiple fixed points throughout its trajectory. while nonstationarities may change the behavior of the system entirely with new fixed points emerging or disappearing. Simple preprocessing measures, such as mean-centering the data, are inadequate account for these behaviors. To enable robustness against these behaviors, we introduce time-varying offset term b. R into Equation (2) as flexible way account for dynamics not readily captured by the original dLDS latent dynamics model.
Let the transition between any two state vectors Dt, L141 1 E NR- be defined by the linear dynamics matrix F,- E -.NNN and the dynamics offset b R... For any X 0, the objective,
Lemma 1.
arg min | D. F,N, b + F| F.b.
is minimized when F. O and b, = 01+1 Tt:-
This result, proven in Appendix B.I reveals that introducing time-varying offset term makes inference of the dynamics degenerate problem. While the solution in Lemma [ minimizes the objective, fails to capture any meaningful structure in F, as result of b, being unconstrained. To prevent the convergence to these trivial solutions, we decompose the latent state space
4 b.
where L, captures fast dynamics and b, captures slow-varyine varying trend behavior. These latent variables follow the dynamic equations 4+1 = L, F,I, + 61. and b = b, 6b,11 where E,, Eg R represent noise sampled from EG N/ Es) and A E,) respectively.

---

3.2 Probabilistic Time-Informed Sparsity
In decomposed models, we aim to achieve two goals simultaneously: sparsity and smoothness of coefficients over time. Motivated by this, we incorporate dynamics-informed probabilistic structure. First, we assume that each coefficient evolves independently of the others. Second, we introduce hierarchical variance parameter ' that controls the sparsity for each Ct Moreover, we introduce dynamics information during sparse inference by encouraging similar active support set in consec- utive time slices through the variance hyperpriors. Put together, the resulting coefficient transition density in p-dLDS becomes,
II - k=1
(5)
We define the first term the right-hand side with the following functional form,
c.
p(ot, 1e) exp
N(
(6)
2yr,k
207-1,
This density captures the constraints of sparsity and smoothness for the inferred coefficients Cr,k- When the variance around zero "t, is small, this structure promotes sparsity by shrinking coefficient values towards zero. Conversely, when the variance around the previous time step - is small, encourages smooothness by shrinking coefficients towards the previous value. While the idea of combining two shrinkage effects in single density has been explored in previous works [15 5 24
those approaches generally require manual balancing of the two penalties. In contrast, we devise procedure in the following section that estimates these variance parameters automatically during inference and learning.
The second density the right-hand side from equation (5) is defined similarly to the hyperprior in IG(C, 1,)- where & weighs the influence the dynamics when estimating kt,)). The resulting graphical model is shown in Figure [ We note that since the value of the previous coefficient is squared, the overall prior placed on the inverse gamma density follows x2 distribution.
SBL-DF. (i.e., p(",,x/et-1..
3.3 Inference and Learning
The joint distribution of p-dLDS is given by,
T
p(a, cl,,),c,)00) p(x1)
II p(y mx))
TT-
K II
II
p( xc, ca)
p(c Ca, (441))
where we denote = #1.7 for brevity. Exact posterior inference is intractable due to the noncon- jugacy introduced by incorporating time-informed sparsity-inducing structure into the graphical model. As result, we devise a variational expectation maximization (yEM) procedure where the approximate posterior factorized
p(æ, c, 7,,, g())g(c,7).
Here, the parameters are given by 0 E [fi.K E. Ex. Ee). Our approach contrasts with BPDN-DF, which estimates latent variables through separate & problems at each point in time. Instead, we preserve the dime-dependence structure within each class of latent variables and leverage efficient inference algorithms that marginalize over uncertainty with respect to time. In general, we seek to maximize the variational lower bound.
(log:) pwy x, c, 7(0 - log (()))((c))))) 7)], (0) = (y(a)g(e77)
with coordinate ascent updates on the latent state posterior, the dynamics coefficients posterior, and the model parameters.
5

---

Updating Latent State Posterior. The optimal coordinate ascent variational update is given by, ((x) x exp (By(er) [log p(æ, y. c, y00)))
(8)
Our assumed decomposition in equation 4 allows us to define the latent state transition density as p( p Tr. Ct) p( 4+1 + bi+1) W4 + Fil, + b,). Substituting this into
equations and (8) we get that the optimal coordinate ascent approximate posterior becomes,
T II Ny: Da, + d, Mn)
T II
((x:) = (a01 E1)
(21:444 F1-14-1 + b,, Ex.
where, 14l and E₁ the mean and covariance of the initial state.
Lemma 2. Let L, b RN be independent random variables such that ~ p(l) and b ~ p(b). Their sum = + b is distributed according to N(p + ps. Eg Ex) when 1) p(b) = Nus, Ex) and p(t) = (M, Σi) and when 2) p(b) 5(b (16) and p(t) (N,, E, -b))
We leverage Lemma 2 reparameterize the state space trajectories into deterministic and stochastic component. The deterministic component captures the slow-moving offset density ((6)
8(b, b,), where b, is estimated using moving average of window size S which can be efficiently parallelized. The remaining dynamics are captured in the stochastic component. We define the family of variational distributions to be the class of linear Gaussian state space models, such that (l) W4,. Σn T W44-1 + 11-14-1, Ez). Conditioned on estimates b and samples of the optimal coordinate ascent variational update for g(l) is efficiently computed using the Kalman We provide a full derivation of the update rule and Lemma 2 in Appendix [3.2]
Smoother (35)
Updating Dynamics Coefficient Posterior. Sparse probabilistic representations introduce non- Gaussian factors which prevent closed-form message passing inference. Specifically, nonconjugacy arises from the inverse gamma term p( (//t+,,,/c. since it is parameterized by ct ~ x². Moreover, the posterior distribution over the coefficients is highly multi-moda modal as a result of the implicit t- distribution in hierarchical model. To update the coefficient posteriors, procedure, where we factorize q(c, 7) = g(c)g(7). First, we obtain an initial estimate of the variational distributions using SBL-DF. Second, we update g(c,,s) = (CT, 7e,k) using stochastic gradient descent (SGD) over
propose three-step
T
c² arg max
log logp(141 ct) + log p(ce41/c, 91+1) + log p(fe+1le))
(10)
where we have estimated the expectation in the optimal coordinate ascent update rule using samples from & ~ 4(7) and ~ q()) To retain coefficient sparsity, we only update coefficients within the active support set. In work, this is defined as coefficients that have an initial estimate of [ct,u] where 1) 10-4, Finally, we update g(>)) based closed form conjugacy rules.
Update Parameters. Given our updated posteriors of the latent variables, we proceed to update the model parameters based on the ELBO,
0" arg max))yta)((ehet [p(y, c, æ, y/0) log ((=))(())(()) = arg magmaxp(y. C, 710), (11)
where we estimate the expectation with samples from our variational distributions and drop terms not dependent on 6. We use SGD to update all model parameters, which possible when we assume that the covariance matrices have diagonal structure.
4 Results
We demonstrate P-dLDS in a variety of synthetic examples, highlighting improved robustness to noise and system nonlinearity. Additionally, we apply our model to clinical neurophysiology dataset, revealing interpretable patterns where previous methods fail. We compare our method against SLDS, iSLDS, and dLDS as described in Section 2 All datasets are split 50:50 for training and testing. Due to space constraints, we provide full deseriptions of the simulation setup in Appendix Cland metric definitions in Appendix

---

A)
NASCAR
B
latent Stare Space
D
State
V1o
Y10
Ubeeree y
Observed
ULDS
p-ILDS
2
Speeed
200
9m
500
Figure 2: Probabilistic model and offset term reduce estimation errors. (A) Example trial from the NASCAR experiment colored by the true switching labels (not provided during training). Each track segment has a random speed т. (B) Inferred state space, colored by discrete state or dominant coefficients. p-dI.DS identifies correct track segments. (C) Example trial from the Lorenz experiment. The speed ramps according to the time intervals År in an ODE solver. (D) Inferred state space, colored by the dominant coefficients. The time-varying offset term allows p-dLDD coefficients to switch according to the true speed and accurately model the two fixed points in the opposing lobes.
Table 1: Metrics for synthetic dynamical systems. Bold means hest performance. (1) indicates higher score is better while (J) indicates that lower is better. indicates that value diverged towards -o0. Switch events for decomposed models are defined times where the active set of DOs change from the previous time step.
NASCAR
Lorenz
Dynamics MSE (1)
Switch MSE (4)
100-step R² (t)
Dynamics MSE (t)
Switch MSE (1)
100-step
Model
(x10-*)
(x10-3)
R² (t)
SL.DS
0.0995
12.89
0.184
0.431
0.0204
-3.47
rSLDS
0.1065
13.17
0.238
0.304
0.0208
-11.54
dLDS
123.19
13.28
X
1.123
0.1529
X
p-dLDS (ours)
0.033
7.34
0.450
0.141
0.0137
0.418
4.1 Synthetic Dynamical Systems
NASCAR with Random Speeds. We evaluate our inference procedure on the NASCAR dataset [26 Since this system easily mean-centered, we can isolate the effect of our proposed inference procedure as offset terms are not necessary. To make this dataset more realistic, we introduce speed variability into the ground truth dynamics opposed having a perfect constant speed all time. Specifically, whenever a trajectory enters new segment of the track, the system experiences random change in speed. We trained all models on 30 trials, each consisting of 1000 time steps (Fig. A with randomly sampled initial points, and randomly constructed 10-dimensional linear observation matrix (see Appendix C. for full details). Performance evaluated on held-out trials where the ground truth switching states are defined by the different segments of the track. In all models, we set the M = 10, N = and K = 4 for DOs or switching states. our experiments, define the "discrete states" for decomposed models as the DO state with the largest coefficient magnitude.
32 23).
Figure 2 shows that changes in the system speed mask the true transition behavior between segments of the track iSLDS. Moreover, dLDS identifies coherent segments, but inappropriately learns different switching pattern for outer and inner edges of the track. In contrast, p-dLDS identifies switching pattern most consistent with the true track segments despite the presence noise and randomness in the system's 's speed. We note that while there are four true segments, decomposed models form more parsimonious representation by identifying similar behaviors in different track segments such as in both edges and curves.
7

---

Table summarizzs our quantitative evaluations on three metrics: 1) the mean squared enror (MSE) between the learned and ground truth latent dynamics, 2) the MSE between the inferred and true switch rate to determine agreement of the discrete switching behavior, and 3) the 100-step inference R2 demonstrate that the learned system generalizes beyond a single step on held-out data, (See Appendix Dorr mathematical definitions). We see that p-dL.DS broadly outperforms existing methods in all metrics and significantly improves inference for decomposed models.
Lorenz System with Random Ramping. Next, we consider the Lorenz system, a chaotic nonlinear system with multiple fixed points, exploring the effect of the offset term. The system described by the differential equation, å = [z æ1), æ1 (p- xz) - 22, 21 - [za] where the parameters 28, 3 8/3, and 10 define a chaotic attractor with two opposing lobes (Fig. 2D). We introduce continuous fluctuations in the underlying dynamics by randomly ramping the system's 's speed throughout each trajectory. This is accomplished by adjusting the evaluation time intervals given to an ODE solver. Similar to before, we randomly construct a linear observation function with M 10 and train models on 30 randomly constructed trials with 1000 time points each (see Appendix C.2). Furthermore, we define the ground truth switch events as the time points when the signal transitions between the two lobes in addition to the moments when a ramping period concludes. All models are trained with a latent space of N = 3 and K = 4 states or DOs.
In figure 2p we see that rSLDS does distinguish between the different speeds along the outer and inner sections of the attractor. Instead, the discrete states obscure the continuum of speeds by incorrectly grouping all activity in each lobe into single regime. Furthermore, we observe that dLDS is limited without an offset term, unable to accurately represent multiple fixed points. Instead of aligning with the two attractor lobes, transitions in the dominant coefficients occur radially relative to the origin and fail to reconstruct the two orbiting fixed points, Conversely, p-dLDS's offset term enables learning a system where coefficients better match the true geometry. This representation correctly recovers differences between the outer and inner sections of the attractor while also accurately reconstructing the two orbiting fixed points. Moreover, this leads to improved estimation of latent dynamics, a switching rate that agrees with the true system, and improved multistep inference performance as shown in Table
4.2 Simulated Motor Cortex Data in a Reaching Task
A)
C)
Men
D)
S.DS
P.CLDS
2717 PU o 1999
2267
Beachh Diveckat
Messored EMG
Fartonartos
Time
Figure 3: P-dLDS efficiently captures changes in dynamics. (A) Latent factors computed from empirical EMG of the reaching experiment in [8. Dynamics are characterized by a preparatory and movement phase. (B) Synthetic spikes and LFPs are generated using the wslfp package [17 [6] (C) The trial-averaged coefficients for p-dLDS smoothly vary with reaching angle. DO captures preparatory dynamics while DO 3 captures movement dynamics. (D) Confusion matrix for linear classification of reach directions. p-dLDS predictions closely align to true diagonal.
We now turn an empirically-derive synthetic experiment related to brain comi--c interfaces, where the dynamics and observation functions are nonlinear and derived from analysis of neural data. Our focus is on the reaching task, neuroscience experiment designed to study motor control in non- human primates [99 [22 In this experiment, the subjects are trained to reach towards visually cued targets, while neural activity recorded from motor-related areas such as electromyography (EMG) data from arm muscles. Each trial consists of two distinct phases: preparation and movement. In the preparation phase, the subject plans its movement while keeping their arm still. the movement phase, the subject physically reaches towards the target. The goal of this experiment is to decode reach intention from neural data. We construct dataset by first simulating spiking neural network with known latent factors
[8] trained to reproduce empirical EMG signals from the center-out reach
8

---

Table 2:
Inference performance for the reaching experiment (see Figure 3 on held-out test set. Top-1 and Top-3 accuracies are obtained by predicting reach directions from latent variable features using linear classifiers. State and Dynamics MSE are computed with respect to true latent variables.
State MSE (L)
Dynamics MSE (4)
Model
Top-1 Acc. (f)
Top-3 Acc. (t)
(x10-1)
(x10-2)
SLDS
38.46
57.69
0.5289
0.3942
rSLDS
12.82
32.05
0.5503
292.41
dLDDS
10.25
39.74
0.6742
35.680
pdLDS (ours)
42.31
70.51
0.4061
0.0567
task in [22]] Spikes then converted to 50-channel local field potentials (LFP) recordings via
weighted, delayed sum of synaptic currents (see Figure 3/ and B) [29 [7 [6]. Our dataset contains 150 6-second trials sampled 250 hz, where each trial represents one out of eight reach directions visually cued at a random start time. PCA identifies that three components captures 98% of the variance. Thus, we M 50, N 3, and K DOs discrete states.
Figure 3C shows the trial-averaged DO coefficients from p-dLDS, which change smoothly and cyclically according to the true reach angle. Additionally, the DOs appear to differentiate between the two distinct dynamical regimes, where the activity of fi and fs localize to the preparatory and movement phase respectively. Quantitatively, reach angles using the average state activity over time as features (see Appendix D.3). FigureBD shows that classifiers built from SLDS features fail to capture the full continuum of reach angles. This limitation occurs because the discrete switching states are unable to efficiently capture the diversity of activity present the LFPs, which arise from randomness inherent in the spike sampling process. Consequently, the inferred features from SLDS generalize poorly to held-out data. In contrast, the p-dLDS classifier predictions recover the full spectrum of reach angles since features are naturally continuous and the inferred coefficients can adjust the learned DOs to accurately capture the activity in the held out data. Table shows that p-dLDS outperforms all other models in state and dynamics reconstruction as well as the top-1 and top-3 reach classification accuracy.
compute the linear classification accuracy of the
First Half
Recording Session (Train)
D)
Second Half of Recording Session (Test)
Emotional Content
Watch
Descrbu
Fao
Pooling
Pos liv
Negatied
40
LFP
Spectrogram (H2)
20
C)
Sone
dLDS Coafs.
p-di DS
Cocfs.
200
300
400
500
600
1000
1100
1200
1300
1400
1500
Time
Time
Figure 4: Learned system discovers coherent structure in clinical data. (A, D) LFP data was collected on patients watching videos with different emotional content. (B, E) LFP spectrograms 40-dimensional signals where each channel represents a particular frequency. (C, F) Inferred states and coefficients shows that rSLDS and dL.DS exhibit unpredictable switching behavior. In contrast, p-dLDS captures smooth coefficients and identifies DOs that align with the trial's emotional content. The learned patterns broadly generalize the held-out data.

---

4.3 Clinical Neurophysiology Data
We demonstrate p-dLDS on LFP recordings from the subcallosal cingulate cortex (SCC) in patients with ceatment-resistan depression "ClinicalTriils... identifier NCT01984710), Subjects are asked watch videos with different emotional content (positive, negative, and neutral), describe the videos, and then discuss how the video made them feel. SCC dynamics have been previously shown to provide a quantitative signal for the presence of emotional content [4] and depression recovery []] Thus, we hypothesize that the underlying dynamics may provide information about emotional changes throughout the experiment. We apply p-dLDS to single patient's LFP spectrogram data (Fig. 43. E) within the 0-40 Hz frequency range (M 40). PCA indicates that the first 7 components explain 90% the variance. Therefore, we train a model with = 7 latent dimension and K DOs.
In Figure 4 rSLDS and dLDS produces high degree of state oscillations making it difficult to identify time intervals with consistent emotional content. In contrast, p-dLDS infers coherent structure that corresponds to changes in emotional content in the trial. For example, fi (red) coincides with resting, f2 (orange) with positive videos, and fi and fs (blue and green) to negative videos (Fig. (C). Importantly, this structure persists even on held out data from the second half of the session (Fig. F). We note this preliminary analysis on single subject isn't intended to make claim about specific neurophy responses to emotional content in this brain region, but generally highlights that p-dLDS identifies meaningful dynamical modes where previous models are unable to.
Conclusion
In this work, we present probabilistic decomposed linear dynamical systems model that can be used to discover meaningful representations in neural signals. By marginalizing over uncertainty in latent variable estimates and incorporating an offset into the dynamics, we enhance robustness and improve a variety of performance metrics. Some areas of future work includes exploiting structure in the offsets to automatically identify window size and extending the probabilistic model include more complicated emissions distributions, such as the Poisson likelihood commonly used to model neural spiking data (28).
Acknowledgements and Disclosure of Funding Sources
Y.C. and C.R. were funded by the James S. McDonnell Foundation (grant number 22002039), with Y.C. being further funded by National Institutes of Health (grant number 2T32EB025816) 16), and C.R. being further funded by the Julian T. Hightower Chair. Y.C. and K.J. were part of the Georgia Tech/Emory NIH/NIBIB Training Program Computational Neural Engineering (T32EB0258 16). N.M. was funded by The Kavli Foundation NeuroData Discovery award. A.S.C. were partially supported by the NSF CAREER Award (2340338) and a Johns Hopkins Bridge Grant. S.A. is supported by the National Center for Advancing Translational Sciences of the National Institutes of Health (Award Number ULITR002378 and KL2TR002381), The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.
References
[1] Sankaraleengam Alagapan, Ki Sueng Choi, Stephen Heisig, Patricio Riva-Posse, Andrea Crowell, Vincet Tiruvadi, Mosadoluwa Obatusin, Ashan Veerakumar, Allison C. Waters, Robert E. Gross, Sinead Quinn, Lydia Denison, Matthew O' O'Shaughnessy, Marissa Connor, Gregory Canal, Jungho Cha, Rachel Hershenberg. Tanya Nauvel, Faical Isbaine, Muham- mad Furqan Afzal, Martijn Figee, Brian H. Kopell, Robert Butera, Helen S. Mayberg, and Christopher J. Rozell. Cingulate dynamics track depression recovery with deep brain stim- ulation. Nature, Sep 2023. ISSN 1476-4687. doi: 111003//1111-002-66666363 https://doi.oo htt////oooooro//////// 10. 1038/ 0241188-0223-6641--1
URL
[2] Kyounghwa Bac and Bani K Mallick. Gene selection using two-level hierarchical bayesian model. Bioinformatics, 20(18):3423 34333 2004.
[3] Leonard Berrada, Andrew Zisserman, and M Pawan Kumar. Smooth loss functions for deep top-k classification, arXiv preprint arXiv:1802.07535 2018.
10

---

[4] Carlos M Carvalho, Nicholas G Polson, and James G Scott. The horseshoe estimator for sparse signals. Biometrika, 97(2):465-480, 2010.
[5] George Casella, Malay Ghosh, Jeff Gill, and Minjung Kyung. Penalized regression, standard errors, and bayesian lassos. 2010.
[6] Adam Charles, Aurele Balavoine, and Christopher Rozell. Dynamic filtering of time-varying sparse signals via & minimization. IEEE Transactions on Signal Processing, 64(21):5644-5655 2016.
[7] Gustavo Deco, Morten L Kringelbach, Viktor K Jirsa, and Petra Ritter. The dynamics of resting fluctuations the brain: metastability and its dynamical cortical core. Scientific reports, 7(1): 3095, 2017.
[8] Brian DePasquale, David Sussillo, LF Abbott, and Mark M Churchland. The centrality of population- factors to network computation is demonstrated by versatile approach for training spiking networks. Neuron, 111(5):631 -649, 2023.
[9] John Dormand and Peter J Prince. A family of embedded runge-kutta formulae. Journal of computational and applied mathematics, 6(1):19-26, 1980.
[10] Mário AT Figueiredo. Adaptive sparseness for supervised learning. IEEE transactions on pattern analysis and machine intelligence, 25(9): 150-1159, 2003.
[11] Michael D Fox and Mareus E Raichle. Spontaneous fluctuations in brain activity observed with functional magnetic resonance imaging. Nature reviews neuroscience, 8(9):700-711, 2007.
[12] Zachary T Goodman, Jason S Nomi, Salome Kornfeld, Taylor Bolt, Roger A Saumure, Celia Romero, Sierra A Bainter, and Lucina Q Uddin. Brain signal variability and executive functions across the life span. Network Neuroscience, pages 1-15, 2024.
[13] Trevor Hastie, Robert Tibshirani, Jerome Friedman, and Jerome H Friedman. The elements of statistical learning: data mining, inference, and prediction, volume 2. Springer, 2009.
[14] Julius Huebl, Christof Brücke, Angela Merkl, Malek Bajbouj. Gerd-Helge Schneider, and Andrea A. Kühn. Processing of emotional stimuli is reflected by modulations of beta band activity in the subgenual anterior cingulate cortex in patients with treatment resistant depression. Social Cognitive and Affective Neuroxcience, 1 1290: 03 2016. ISSN 1749-5016. doi:
10. 10.1099//cca//sw0033. URL htthh:://////o
10. http:////....000///aaa//////a//
[15] Kaoru Irie. Bayesian dynamic fused lasso. arXiv preprint arXiv:1905.12277: 2019.
2019,
[16] Kyle Johnsen, and andiiaieeeee Siplab-gt/wslfp: V0.2.1. Zenodo, March 2024. URL
[17] Kyle A Johnsen, Nathanael A Cruzado, Adam A Willats, and Christopher J Rozell. Cleo: A testbed for bridging model and experiment by simulating closed-loop stimulation, electrode recording, and optogenetics. bioRxiv, pages 2023-01, 2023.
[18] Yuko Kakikawa, Kaito Shimamura, and Shuichi Kawano. Bayesian fused lasso modeling via horseshoe prior. Japanese Journal of Statistics and Data Science, 6(2):705-727, 2023.
[19] Matthew T Kaufman, Mark M Churchland, and Krishna V Shenoy. The roles of monkey m1 neuron classes movement preparation and execution. Journal of neurophysiology, 110(4): 817-825, 2013.
[20] Timothy Doycon Kim, Thomas Zhihao Luo, Jonathan W. Pillow, and Carlos D. Brody. Infer- ring latent dynamics underlying neural population activity via neural differential equations. Proceedings of the 38th International Conference on Machine Learning, 2021.
[21] Maksim Lapin, Matthias Hein, and Bernt Schiele. Loss functions for top-k error: Analysis and insights. In Proceedings of the IEFE conference on computer vision and pattern recognition, pages 1468-1477) 2016.
11

---

[22] A. H. Lara, J. P. Cunningham, and M. M. Churchland. Different population dynamics in the supplementary motor area and motor cortex during reaching. Nature Communications, 91):2775, July 2018. ISSN 2041-1723. doi: 10.1038/s41450 1477-01--0514----1 URL. https 7/ //www //www.nateet natura. com/ ctt.a//rrricccae//11400000001110
[23] Hyun Dong Lee, Andrew Warrington, Joshua I Glaser, and Scott W Linderman. Switching autoregressive low-rank tensor models. arXiv preprini arXivv2306.03299, 2023.
[24] Qing Li and Nan Lin. The bayesian elastic net. 2010.
[25] Scott Linderman, Benjamin Antin, David Zoltowski, and Joshua Glaser. SSM: Bayesian Learning and Inference for State Space Models, October 2020. URL htlpttts///tthheie com/ lindermanlabe aaaaanllbbaae lbbmsm
[26] Scott W Linderman, Andrew C Miller, Ryan P Adams, David M Blei, Liam Paninski, and Matthew J Johnson. Recurrent switching linear dynamical systems. arXiv preprint arXiv:1610.08466 2016.
[27] Yuxi Liu, Xian Long, Paul R Martin, Samuel G Solomon, and Pulin Gong. Lévy walk dynamies explain gamma burst patterns in primate cerebral cortex. Communications Biology, 4(1):739, 2021.
[28] Jakob H Macke, Lars Buesing, 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.
[29] Alberto Mazzoni, Henrik Lindén, Hermann Cuntz, Anders Lansner, Stefano Panzeri, and Gaute T. Einevoll. Computing the Local Field Potential (LFP) from Integrate-and-Fing Network Models. PLOS Computational Biology, 11(12)e1004584, 2015. ISSN 1553-7358, doi: SESTIIJOUUNNALPCCL10000
URL htttt:://jjooaaaaalll plos org/ploscompbiol/ article?id-10 1 pcbi 1004584
[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): 44, 2024.
[31] Nandakumar S Narayanan. Ramping activity is cortical mechanism of temporal control of action. Current opinion in behavioral sciences,
2016.
[32] Josue Nassar, Scott W. Linderman, Monica Bugallo, and Il Memming Park. Tree-structured recurrent switching linear dynamical systems for multi-scale modeling. In International Confer- ence Learning Representations (ICLR), 2019. HkZRQQR9YX
[33] Matthew R O' S'haughnessy, Mark A Davenport, and Christopher J Rozell. Sparse bayesian learning with dynamic filtering for inference of time-varying sparse signals. IEEE Transactions on Signal Processing, 68:388-403 2019.
[34] Chethan Pandarinath, Daniel J O' Shea, Jasmine Collins, Rafal Jozefowiez, Sergey 1) Stavisky, Jonathan C Kao, Eric M Trautmann, Matthew Kaufman, Stephen Ryu, Leigh Hochberg, al. Inferring sinale-trial neural population dynamics using sequential auto-encoders. Nature methods, 15(10)=805-88 815, 2018.
[35] Herbert E Rauch, F Tung, and Charlotte T Striebel. Maximum likelihood estimates of linear dynamic systems. AIAA journal, 3(8):1445-1499 1965.
[36] Maria Ribeiro and Miguel Castelo-Branco. Slow fluctuations in ongoing brain activity decrease in amplitude with ageing yet their impact on task-related evoked responses is dissociable from behavior. Elife, 1 1572222 2022.
[37] Michael E Rule, Carlos Vargas-Irwin, John P Donoghue, and Wilson Truccolo. Contribution of Ifp dynamics to single-neuron spiking variability motor cortex during movement execution. Frontiers in systems neuroscience, 9:89, 2015.
12

---

[38] Steffen Schneider, Jin Hwa Lee, and Mackenzie Weygandt Mathis. Learnable latent embeddings for joint behavioural and neural analysis. Nature, 617(7960):360-666 2023.
[39] James M Shine, Patrick G Bissett, Peter Bell, Oluwasanmi Koyejo, Joshua H Balsters, Krzysztof J Gorgolewski, Craig A Moodie, and Russell A Poldrack. The dynamics of functional brain networks: integrated network states during cognitive task performance. Neuron, 92(2): 544-554, 2016.
[40] Chuammeizhi Wang, Bijan Pesaran, and Maryam M Shanechi. Modeling multiscale causal interactions between spiking and field potential signals during behavior. Journal of neural engineering, 19(2):020001, 2022.
[41] Leonhard Waschke, Niels Kloosterman, Jonas Obleser, and Douglas D Garrett, Behavior needs neural variability. Neuron, 109(5):751-766, 2021.
[42] Wei Wu, Michael J Black, David Mumford, Yun Gao, Elie Bienenstock, and John P Donoghue. Modeling and decoding motor cortical activity using a switching kalman filter. IEEE transac- tions on biomedical engineering, 51(6):933-942, 2004.
[43] Andrew Zalesky, Alex Fornio,, Luca Cocchi, Leonardo L Gollo, and Michael Breakspear. Time-resolved resting-state brain networks. Proceedings of the National Academy of Sciences, 111(28):10341-10444 2014.
[44] David M Zoltowski, Jonathan W Pillow, and Scott W Linderman. Unifying and generalizing models of neural dynamics during decision-making. arXiv preprint arXiv:2001.045// 2020.
13

---

A Appendix supplemental material
A.1 p-dLDS Algorithm
Algorithm | describes the proposed inference algorithm. In our experiments, we set n = and " 10-1 and observe that the model converges. Below use the notation hat notation for latent variable estimates or samples and the variable itself to represent the parameters of the variational distributions.
Algorithm 1 Variational EM for Probabilistic dLDS
Require: M observation dimension, N latent state dimension, К number of dynamic operators, S moving average window size, SBL-DF trade-off parameter, n number of samples to estimate expectations, 1/ sparsity threshold, 0 model parameters.
Initialize parameters
C: 1--00 0 Di.j ~N 0²)
fK, ~N(0,0²)
Et f- D+
b Initialize latent state with PCA
while ELBO has not converged do
Update Latent State Posterior
burt MoringAverage.(Garn)
C2 ((())
F. Σk k== feCea
Lr Σ, +- KalmanSmootther(pzr by.r, 0)
// Update Coefficient Posterior Initialize q(c) and ((>) jointly with SBL-DF.
Update g(cr,k) with SGD over equation
10 for densities where |cr, >
Update ((3)) + IG(E + 2: 5.
En (0, k. ox,))*
// Update Parameters
Update with SGD over equation end while
B Latent Variable Inference
B.1 Lemma Derivation
Lemma 1.
Let the transition between any two state vectors Br, F++1 E ER be defined by the linear
dynamics matrix F; E R N and the dynamics offset b, E RV For any X > 0, the objective,
arg min [ T. F,æ, b, + F||
F.
is minimized when F, 0 and b = x1+1
Proof. Let rt F++1 .t. We can rewrite the reconstruction objective the following form,
arg min ||r, Fi byllă + X F
F.
This objective is identical to the standard ridge regression with an unpenalized intercept term [3 The solution is obtained by first centering data, and then solving for the parameters using solution for the standard Tikhonov regression. Below, we define the centered data as á, and 7g for inputs and outputs respectively. Finally, we can use these values obtain the following estimates of the parameters,
14

---

d. p
F. ( XT) -
However, when there is only a single datapoint, we get that ã; = 0, which results in O.
This result arises from having only single observation for any dynamic transition, which leads to a singular design matrix. Although we can improve our estimate of F, by collecting more samples along given trajectory, this is impractical when dealing with naturalistic time-series. For instance, it may be infeasible to collect more data from the exact same initial condition in a naturalistic environment due to noise the experimental setup. chaotic systems, minor deviations can lead to drastically different outcomes over long time horizons. Even if it were possible to precisely control for the initial condition of the signal, the presence of dynamical noise can cause initially aligned time series to quickly drift out of alignment. Consequently, it is not uncommon to observe a single transition between any two time points, as not guaranteed that events across multiple trials will be well-aligned.
B.2 Lemma Derivation
Lemma 2. Let L, 6 E RN be independent random variables such that ~ p(l) and b p(b). Their sum x = 1 b is distributed according to N (pn + ph, Ei + Eb) when 1) p(b) = N Ex) and p(t) (u,, Σi) and when 2) p(b) 8(b ps)) and p(l) N (ms, Ex Ex).
Proof. Case 1. Let 1 ~ N(s, Σd) and ~ (us, Es). The sum of normal random variables follows distribution that results from convolving their individual distributions,
g(x) = q(l + b)
q(l) q(b)
-NW, Σi) NN Es)
=NWm pb: El + Eg)
This is standard result from probability theory.
Case 2. Now let 1 ~ (w,, E, + Σ.) and b ~ 8(b 1a) Similarly, the distribution of the sum of these variables is distributed according their convolution,
q(x) = ql) q(b)
NM Ei + Σn) 8(b 1u5))
N(x T; Mi, E, + ET)5(T - pb)dT
= =N(x + pb: pı, EΙ + Ex)
=N(x; it pb, Σ Ex),
where the fourth line is the result the sifting property delta distributions. Since the final distribution in Case and Case 2 are identical, we complete the proof.
B.3
Optimal g(x) Update
The optimal coordinate ascent variational update is given by the following equation, log q" (æ) x loggp e, y. X|0)]
(12)
(og (l + Σ logp(1. Cr. 0) + Σ log logp(y.l + b., 0)) C.
Conditioned on estimates of b1.T and samples CI:T, the factor graph of equation 12 corresponds exactly to a time-varing Linear Gaussian State Space Model. Thus we can leverage the efficient inference algorithms such as the Kalman filter and RTS smoother when computing marginals of the variational distribution of L.T.
15

---

C Generating Synthetic Examples
C.1 Noisy NASCAR
NASCAR data is generated by partitioning the two-dimensional state space into four regions according
to the rules,
1, 21 > 1
2, 3 < -1
((x)
3, -1 < x1 1,,02
4, -1 -15x1 x E 1, < 0,
where ((x) is the ground truth switching state function that depends on the particular location æ. The ground truth dynamics matrices defined as,
0.1
when Z(x) 1 or
--1
0
A(x) =
To
when Z(x) = or 4,
0
0
and ground truth offsets are defined as,
0 0.005
when Z(x)
[0
-0.005
when Z(x)
b(x)
[0.1 0]
when Z(x)
[-0.1 0] when Z(x) = 4.
Given the current location in state space æt, we can transition to the next point using the continuous time dynamics equation
3t = expm(p (T ((..))) + Tb2(..) + Us.
where each entry of the process noise is sampled from Vt, NoO. 10-4). To modulate the speed of the system, we uniformly sample a speed constant E (0.1, 1]. which applied throughout each segment of the track. We use the continuous time formulation over the discrete-time formulation to ensure that changes the speed do not distort the shape of the original system's state space. To generate noisy observations, we construct linear emissions matrix with random variables such that each entry given by Da.j (O, 1).
C.2 Ramping Lorenz
In order to modulate the speed of the Lorenz system, we adjust the evaluation time points of an ODE integrator, specifically Kunge-Kutta of the order 5(4) (RK54) as implemented in scipy's solve_ine ivp [9. Ramping activity is generated randomly with the following procedure,
1. Uniformly sample an evaluation interval length E [0.25, 1.5).
2. Construct vector 7 that consists evenly spaced numbers over the interval [0, r). In our experiments, we set n to be 100.
3. Perform the transformation exp(7) to obtain vector of ramped evaluation times.
4. Plug in the transformed evaluation times into the RK45 Solver to obtain latent trajectories.
Similar to the NASCAR experiment, we generate noisy observation from a randomly constructed linear emissions matrix such that each entry is given by Dijj N/ 1).
16

---

C.3 Simulated Monkey Reaching Task
Our dataset constructed from publicly available data and code from the center-out reach task in We obtain latent factors from spiking networks that trained to reproduce empirically measured EMG signals, given a 3-dimensional input that specifies the go input and the reach angle. In experiments, these factors considered ground truth. Our trained factor-based spiking network then generates spiking activity for 1200 neurons. Synaptic currents are used as inputs into the Weighted Sum of synaptie currents LFP proxy method (WSLFP) [29], implemented in the wslfp Python package [17 16]. As WSLFP a function of the relative location of neurons and electrodes, place neurons randomly within mm by 10 mm by mm region and electrodes in grid centered in this region. The result is multi-channel LFP dataset with nonlinear dynamics and measurements characteristic of systems neuroscience.
12
A)
C)
Proph
How.
Neurons
Y15
Electrode
LFP Channel
Neuron ID
Yo
Time (s)
Figure 5: Empirícally -Derived Reach Experiment. (A) 1,200 neurons are randomly placed into
5 mm by 10 mm by mm region. Electrodes are placed in grid centered in this region (B) Spiking activity for subset of neurons in example trial produced from factor-based spiking network. (C) First 15 channels in a simulated multi-channel LFP recording. Preparatory and Movement phases marked by the dotted lines.
D Evaluation Metrics
D.1 Multi-step Inference
The multi-step inference performance is computed with the following R-squared metric,
R
1=
Yotk
(13)
Σ5
BL+k
where & is the number of steps from the initial condition, ÿ is the mean estimator for each trajectory and Vc+k the model prediction after applying the inferred dynamics for k steps. When testing, model parameters such as the dynamics and observation matrices are frozen, while specific latent variables are estimated based on the held-out data. Table we show results for 100.
D.2 Inferred Dynamics Error
We measure the accuracy of the latent dynamics with the mean squared error (MSE) of the inferred speed, defined as,
MSE
+ U
(14)
where the true speed à. = 2z+1 xg is computed from the denoised ground truth latent state, and the predicted speed h.- 1 T is computed using the model's 1-step prediction. Since latent trajectories are only identifiable up to a linear transformation, we align the inferred trajectories with the true trajectories using least squares fit before computing this score. More specifically,
17

---

find the optimal linear transformation U E IR trajectories by solving,
NNN
between the estimated and true states across all
1 T
Ú arg min
Σ
22 - URell-
(15)
T
1=1
D.3 Inferred Latent State Space Error
Similarly, we measure the accuracy of the latent state space by computing the MSE after linear aligmment between trajectories from the inferred and true state space. We use this metric only for the reaching example, since the true observation function is a complex nonlinear function,
T
MSE T
Σ c
U2,12
(16)
The linear alignment U N between the estimated and true states across all trajectories is computed by solving the least squares problem in equation (15)
D.4 Inferred switching rate error
Evaluating the accuracy of the switching behavior is more difficult task. In f'ct, developing procedure that matches predicted switch times with true switch times can lead to complicated optimal transport procedure. To simplify the evaluation of switching times, we marginalize over time, and compare only the MSE of the switch rate defined as,
m
MSEawitch
m Σ I|s - T.122
(17)
where m is the number of trials, r; is the true switch rate for the ith trajectory, and Pi 1 Ei= 1f2, * 2 is the predicted switch rate. Intuitively, F; is the number of times that the state or dominant DO changes between consecutive time points normalized by the length of the interval T. In switching models, switch events are defined as time point where the current inferred dynamical state differs from the state the previous time step. Similarly in decomposed models, switch events are defined as time points where the active set of DOs change from the previous time step.
In the NASCAR example, r; is defined with the number of transitions between ground truth segments. In the Lorenz example, r; is defined by the number of times that the trajectory switches between the two lobes in addition to the number of ramping periods.
D.5 Reaching Classification Accuracy
We quantitatively evaluate the reaching experiment with classification task. Here, we want to determine whether the learned systems can be used to distinguish between different reach directions. Recall that switched models infer a switching variable for each time point where 2 {1..... K} while decomposed models infer coefficient vector E R Rather than viewing 2, we can equivalently view as a one-hot encoded vector 2, E (0, 1] which describes whether a particular switching state is active at any given time. This matches the dimensionality of the variables in both switched and decomposed systems.
index,
For simplicity, we focus on linear logistic regression classifiers in our experiment. If we let the inputs be 2, and e, directly, then our classifiers quickly overfits since there are many more input features than trials. Specifically, the number of features scales linearly with the number of time points and systems O(TK). Instead, we marginalize over time and compute features from the estimated latent variables by averaging state activity over time. In switched models, this is the average one-hot encoding value over time. Similarly, this the average coefficient value in decomposed models. However, for each dynamical state, we compute separate features for positive and negative coefficient values to prevent interference between them. In this setup, the input (feature) dimensionality scales according to O(K) while the output dimensionality of the linear classifiers are the reaching directions. For all
18

---

classifiers, we perform a grid search over the values {10'} to identify an appropriate amount of L2 regularization. Top-k accuracies are standard metric in machine learning [ and computed using the estimated class probabilities from the logistic regression classifier.
E
Additional Results
E.1 Synthetic Dynamical Systems
Figure 6A demonstrates that our inference procedure converges to a local optimium while Figure 6 shows full sweep of the multi-step inference metric. Tables [ in the main paper reports the final value. For completeness, we include Tables 3 and 4 which reports the means across seeds of each model, and includes the standard deviations parenthesis.
NASCAR
Lorenz
B)
NASCAR
Lorenz
700
3000
000
2002
500
1002
100
100
150
Spocha
Esodia
Figure 6: (A) ELBO converges in both synthetic dynamical systems. (B) Multi-step inference where k represents the number of steps. Tables and report the final values.
Table 3: Metries for NASCAR. Bold means best performance. (t) indicates higher score is better while (L) indicates that lower is better. X indicates that value diverged towards -00. All MSE values are x x10 while R²2 values are not scaled. We report means across 5 seeds and include standard deviation in parenthesis.
Model
Speed MSE (4)
Switch MSE ())
100-step R² (t)
SL.DS
0.0995 (0.021)
12.89 (1.30)
0.184 (0.024)
iSLDS
0.1065 (0.024)
13.17 (2.84)
0.238 (0.022)
dLDS
123.19 (23.13)
13.28 (5.31)
X
P-LLDS (ours)
0.033 (0.009)
7.34 (3.40)
0.450 (0.027)
Table 4: Metrics for Lorenz. Bold means best performance, (1) indicates higher score is better while (L) indicates that lower is better. X indicates that value diverged towards -oo. We report means across 5 seeds and include standard deviation in parenthesis.
Model
Speed MSE (4)
Switch MSE (L)
100-step It? (1)
SLDS
0.431 (0.233)
0.0204 (0.007)
-3.47 (L.052)
rSLDS
0.304 (0.040)
0.0208 (0.004)
-11.54 (1.353)
dLDS
1.123 (0.089)
0.1529 (0.070)
X
PDLLDS (ours)
0.141 (0.015)
0.0137 (0.014)
0.418 (0.079)
E.2
Reaching Task
For each model, we visualize the trial-averaged dynamic regime activity of each reach direction (Fig. 7. In SLDS, this is visualized by considering the discrete states as a one hot vector over time. When dynamic regime is active, that state will have value of I while the unactive states will have value of O. Thus the trial averaged value of each state must have a value the interval [o, 1]. In dLDS, we plot the inferred coefficient value without any modification.
Although SLDS correctly identifies preparatory and movement phases using states and respectively, it fails to differentiate dynamics occurring outside of these expected phases, incorrectly grouping
19

---

unrelated regions together. Furthermore. the discrete formulation produces very similar activity patterns across all reach angles, obseuring any differences that are present. In dLDS, we observe that the features change smoothly and cyclically with the reach angle. However, the dynamic operator activity do not localize to the preparatory and movement phases due to limited inference procedure.
Table 5: Inference performance for the reaching experiment (see Figure 3 on held-out test set. Top-1 and Top-3 accuracies are obtained by predicting reach directions from latent variable features using linear classifiers. State and Dynamics MSE are computed with respect to true latent variables. We report standard deviations in parenthesis across 5 seeds.
Model
Top-1 Acc.
Top-3 Acc.
State MSE (x10-1)
Dynamics MSE (x (x10-2)
SLDS
38.46 (2.84)
57.69 (7.53)
0.5289 (0.13)
0.3942 (0.23)
rSLDS
12.82 (3.05)
32.05 (8.31)
0.5503 (0.23)
292.41 (13.96)
dLDS
10.25 (5.97)
39.74 (10.29)
0.6742 (0.52)
35.680 (5.76)
pdLDD (ours)
42.31 (3.50)
70.51 (6.45)
0.4061 (0.38)
0.0567 (0.04)
A)
SLDS Discrete States
B)
dLDS DO Coefficients
State
State 2
D0
DO
Reach Angle
45
90'
State
State
DO
DO
2
2
2
Time 15)
Time (s)
Time (s)
Time [s]
Figure 7: Trial-averaged activity for (A) SLDS discrete states and (B) dLDS DO coefficients for each reach angle. The preparatory and movement phases occur between the dashed lines similar to Figure 3
Time 0 represents the onset of the stimulus.
F Experimental Setup
E1
Hyperparameter Settings
For switching models, we rely on the ssm package which allows for efficient Bayesian inference for variety of state space models (25]. We set the variational posterior to structured p and the fitting procedure to laplace_eed _e as recommended by the developers. Additionally, we set the distributional form of the dynamics and emissions matrices Gaussian.
The hyperparameters of dLDS primarily consists of the lagrange multipliers in the BPDN-DF objective including λς, A1 X2. We find the optimal value of these hyperparameters using a random search with fixed budget of 1000 evaluations. For each hyperparameter, we uniformly sample over the log of the interval [10-3, 103] and evaluate it against the BPDN-DF objective. For the NASCAR experiment, we found that Ào 1.044, a = 0.254, and À = 0.023 resulted in the best performance. For the Lorenz experiment, found that Ao 0.628, λ₁ 2.010, and X2 0.0124 yielded the best performance.
For P-dLDS, the relevant hyperparameters consists of the SBL-DF dynamics tradeoff &. and the offset window size S. We use random search with a budget of 1000 samples to determine the values of S and & and fit a separate model for each set of hyperparameters. In the NASCAR experiment, we isolate the effect of the probabilistic inference procedure by setting = T, removing the influence of the time-varying offset term. For &. we perform a random search by uniformly sample over the log of the interval [10-3 3, 10 and found that E = 0.945 was optimal. For the Lorenz experiment, we also optimize for the window size S by uniformly sample discrete index on the interval [?).
20

---

For the Lorenz experiment, the optimized hyperparameters are S = 85 and & 8.928. For the real dataset, the optimal offset 76 which smaller than the timescale p-dLDS coefficient switching (around 150 time points), suggesting that the same DO dynamics may persist even as the fixed points of the system fluctuates throughout the experiment.
F.2 Hardware Specification
We perform hyperparameter sweep on our institution's HPC cluster using small-scale CPU resources which consists of Dual Intel Xeon Gold 6226 CPUs. Once hyperparameters have been optimized. it is possible to run each experiment within approximately 2 hours on the 2020 edition of the M1 Macbook Pro.
G Description of Clinical Neurophysiology Data
Data was collected as part of a study investigating deep brain stimulation for treatment-resistant depression (TRD). The study is pre-registered in TrinccaTriialsggv (identifier NCT04106466) The study protocol was approved by the IRB (identifier IRB00066843). Informed consent was obtained from participants before participation instead have their DBS electrodes and Summit RC+S IPG device provided free of charge. The analysis focused on LFP signals from a single participant with all personally identifiable information removed.
the trial. Patients receive no monetary compensation, but
H Limitations
While our proposed method demonstrates strong performance in our experiments, there are many limitations. For instance, our approach does not have a strong mechanism for generating future unseen coefficients. Our assumed coefficient transition model is primarily motivated by our desire to obtain smooth coefficients over time. However, we believe that they may be more complex transition models that can both capture persistent activity in challenging systems while also being accurate forecaster, such as deep learning based transition model. Another limitation of our approach is that our method assumes smoothness in the latent space. However, we do not explore the possibility of having sparse structure the latent space which can be easily accomplished BPDN-DF by adding an Ll penalty over
21

---

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: We propose a probabilistic treatment and an extended dynamics formulation in decomposed models (Section 3). We demonstrate that these changes reduce estimation errors and finds coherent structure where previous models fail in many challenging synthetic examples, and noisy real-world example (Section 4).
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in 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 or 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.
It is fine to include aspirational goals as motivation as long as is clear that these goals attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification:
A limitations section is provided in Appendix H due to space constraints.
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 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 (e.g., independence assumptions, noiseless settings, model well-specification, 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 the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect 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-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency 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 fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, 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?
22

---

Answer: [Yes]
Justification:
We provide proofs of our lemmas in B. and B.2
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 proofs 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 be 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 of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
metric definitions in Appendix D and hyperparameter details in Appendix F.I Moreover, we release our code with our submission.
Justification: We provide synthetic data details in Appendix C
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 novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, other means that appropriate to the research performed.
also be provided via detailed
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a 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 reprodueibility may be tricky in some cases, in which authors are welcome to 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 be possible for other researchers have some path reproducing or verifying the results.
5. Open access to data and code
23

---

Question: 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?
Answer: [Yes]
Justification: We provide functions to generate the synthetic datasets from experiments. Unfortunately, we are unable to release the neurophysiological dataset due to request from collaborators.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark).
The instructions should contain exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https for more details.
y /nips. /nips.cc////p https//ccccciigggeeeeeee///////e wwaccc//uudss/eee haabll/c/guideesodeeuuubiiii
The authors should provide instructions on data access and preparation, including how 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, to 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/ettails
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 understand the results?
Answer: [Yes]
Justification: We specify that all datasets are split 50:50 for train and test in the main paper. Moreover, we provide optimizer and hyperparameter settings in Appendix F.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to 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:
We provide standard deviations in Appendix
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
24

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run 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.
OK to report 1-sigma error bars, but one should state 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 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 details about hardware in Appendix
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 of 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 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 httts://neurips.c
Answer: [Yes]
Justification: Most of our datasets are synthetically generated or derived from publicly available sources, which we refer to throughout the paper. For our real-world experiment, we include the IRB identifier.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (c.g.. if there special consid- eration due to laws or regulations in 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 proposed model aims to reveal coherent patterns time-series data, serving as a scientific tool similar to PCA. As probabilistic technique not designed for content generation or automated decision-making, we not anticipate direct societal impacts.
25

---

Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or 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 partieular applications, let alone deployments. However, if there is direct path any negative applications, the authors should point it out, For example, it is legitimate to point that 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 out that a generic algorithm for optimizing neural networks could enable people 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 are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI..
11. Safeguards
Does the paper deseribe 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: Our method does generate synthetic content that risk of abuse. Guidelines:
The answer NA means that the paper poses such risks.
Released models that have 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
Question: 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?
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or 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.
26

---

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 thcode. com/ com/dataseta has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
package should be provided. For popular datasets, paperswi
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 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 include details where obtained our data (either synthetically gen- erated, or part of an existing dataset) in the main text and appendix. We also point to the commonly package used for SLDS and rSLDS in the appendix.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset 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 to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [Yes]
Justification: We provide the IRB identifier in Appendix G which includes the informed consent form with experiment details.
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 the main paper.
According the NeurIPS Code of Ethics, workers involved 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: [Yes]
Justification: We include the IRB identifier number Appendix G but removed any information related to the institution conducting the experiments.
Guidelines:
27

---

The answer NA means that the paper does not involve crowdsoureing nor research with
human subjects.
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 expect authors adhere 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.
28


Reviews:
[{"review_id": "niTnavp9qV", "reviewer": "Reviewer_jjeY", "summary": "The paper presents a probabilistic version of the dLDS model, first presented by Mudrik et al (2024). The primary impetus for this model was to present a version of the dLDS that was more robust to noise, although the authors here also included a slowly-evolving offset term meant to capture evolving fixed points. Inference is conducted by a variational objective which is optimized with a combination of ad-hoc but reasonable methods. The paper paves the way for many other probabilistic state-space models including those with non-Gaussian measurements.", "questions": "I am a bit confused about the decomposition of the coefficient prior defined on lines 168-169 following equation (5). Here the authors claim p(c_t|c_t-1, \\gamma_t) is defined by a factorization N(c_t-1,\\sigma_t-1)N(0,\\gamma_t). The authors state that the first of these factors encourages smoothness while the second encourages sparsity. However, it is not at all clear to me that this makes sense. First, the prior is a conditional distribution over c_t for a single coefficient, so what is the factorization over exactly? Did the authors mean to present a joint distribution? Over which variables? I can see that the authors would like c_t to have both properties (sparsity, smoothness) but this particular specification doesn’t make sense. Second, \\sigma_t and its inference is undefined anywhere in the paper. Third, on lines 215-216 the authors state that parameters can be learned by SGD “which is possible when we assume that the covariance matrices have diagonal structure.” Do the authors specifically mean the posterior covariances (eq 8), the prior covariances over the latent variables (line 158)? Either way, what is the limitation on SGD learning the parameters when these covariance matrices are not diagonal? Lastly, the authors state on line 187 “the approximate posterior becomes…[some Gaussians]”. It is not clear from the text if the Gaussian form of this variational posterior is a choice of the authors or if this is a variational optimum that drops out of the model structure and factorized structure of q. Can the authors please clarify?", "limitations": "No major conceptual limitations but technical details should be clarified", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "Originality - The model is original in the sense that it has not been presented previously.\nQuality - The paper is mostly well written and well reasoned. I see virtually no major conceptual flaws and only minor problems. Both the simulation experiments and analysis demonstration are well conceived and of interest to the computational neuroscience community. The model is well conceived and represents a natural extension of the dLDS model.\nClarity - Very well organized, easy to follow, and a pleasure to read.\nSignificance - I wouldn’t exactly call this a weakness but the extension to dLDS model is so natural as to be virtually obvious. The contribution itself, while seemingly inevitable, was nonetheless executed by the authors first (props to them) and with high quality. It serves as a strong rung in what is sure to be a dLDS ladder.", "weaknesses": "The primary weakness is in some technical confusion that i have the with paper. I detail below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for taking the time to review our paper thoroughly. We are encouraged that you find our work to be \"original\", \"a pleasure to read\", \"well-conceived\", and of \"interest to the computational neuroscience community\".\n\n**Q1 and Q2: Factorization of the coefficient transition and Definition of $\\sigma_t$**\n\nWe clarify that the transition density $p(c_{t,k} | c_{t-1,k}, \\gamma_{t,k} )$ over the coefficients $c_t$ has the following functional form:\n\n$$\n\\begin{align}\n    p(c_{t,k} | c_{t-1, k} , \\gamma_t) &\\propto \\exp \\left(- \\frac{c_{t,k}^2 }{2\\gamma_{t,k}} - \\frac{(c_{t,k} - c_{t-1,k})^2} {2 \\sigma_{t-1,k}^2}\\right) \\\\\\\\\n    &= \\exp \\left(- \\frac{c_{t,k}^2 }{2\\gamma_{t,k}} \\right) \\exp \\left( - \\frac{(c_{t,k}- c_{t-1,k})^2} {2 \\sigma_{t-1,k}^2}\\right) \\\\\\\\\n    &\\propto N(c_{t,k} ; 0, \\gamma_{t,k}) N(c_{t,k} ; \\hat{c}\\_{t-1,k}, \\sigma_{t-1,k}^2)\n\\end{align}\n$$\n\nIn the first line, we propose a density that captures the constraints of sparsity and smoothness for the inferred coefficients $c_{t,k}$. A small  variance around 0, $\\gamma_{t,k}$, promotes sparsity while a small variance around the previous coefficient, $\\sigma^2_{t-1,k}$, promotes smoothness. The second line factorizes this density into two separate terms and the third line shows that each term is proportional to a normal distribution with their respective parameters. \n\nWhile the idea of combining two shrinkage effects in a single density has been explored in previous works [1,2,3,4], those approaches generally require manual balancing of the two penalties. In contrast, we estimate these parameters using SBL-DF results. Specifically, SBL-DF produces estimates for the sparsity variance $\\gamma_{t,k}$ through the estimated hyperparameter posterior $q(\\gamma_{t,k} )$ and the smoothness variance $\\sigma^2_{t-1,k}$ through the approximate coefficient posterior $q(c_{t,k} ) = N(c_{t,k} ; \\hat{c}\\_{t,k} ,\\sigma_{t,k}^2 )$.\n\nWe thank the reviewer for pointing out that our discussion on the coefficient transition densities could be improved and will incorporate this discussion into the camera-ready version if accepted. \n\n[1] Irie, Kaoru. \"Bayesian dynamic fused LASSO.\" 2019.\n\n[2] Casella et al. \"Penalized regression, standard errors, and Bayesian lassos.\" 2010.\n\n[3] Li and Lin. \"The Bayesian elastic net.\" 2010.\n\n[4] Kakikawa et al. \"Bayesian fused lasso modeling via horseshoe prior.\" 2023.\n\n**Q3: SGD over covariance matrices.**\n\nWe clarify that we assume a diagonal structure for the posterior covariances. Estimating general covariance matrices $\\Sigma$ is challenging both mathematically and computationally due to the requirements for a valid $\\Sigma$ (e.g., symmetry, positive semidefiniteness). Vanilla SGD over individual matrix entries does not inherently respect these constraints and does not guarantee that the solution is a valid covariance matrix. Typically, solving for $\\Sigma$ involves a semidefinite program (SDP) requiring specialized methods like ellipsoid or interior point solvers [1]. However, we simplify this problem by considering only when $\\Sigma$ is diagonal (i.e., $\\Sigma = {\\rm diag} (\\sigma_1 ,\\dots, \\sigma_n)$ ). This allows us to simplify general SDP optimization to an unconstrained problem by operating in the log space of the diagonal variance terms. We find that this approach works well in our experiments and is computationally efficient (see additional pdf for runtimes).\n\n [1] Vandenberghe and Boyd. ``Semidefinite programming.'' SIAM review 1996.\n\n**Q4: Derivation of the approximate posterior.**\n\nThank you for pointing out that equation 8 (following line 187) could be improved in clarity. It is the variational optimum that results from the model structure assumed in line 157. We derive this result below. First, we acknowledge two typos:\n\n1. The joint distribution of pdLDS (equation 6) is missing emission densities and will be updated to,\n\n$$ \\begin{align}\np(x, y,c, \\gamma | \\theta ) = p(x_1) \\left[ \\prod_{t=1}^ T p(y_t | x_t) \\right] \\left[ \\prod_{t=1} ^ {T-1} p(x_{t+1} | x_{t}, c_{t})\\left[ \\prod_{k=1}^{T} p(c_{t+1,k} | c_{t,k}, \\gamma_{t+1,k}) p (\\gamma_{t+1, k} | c_{t,k}) \\right] \\right] \\tag{1}\n\\end{align}$$\n\n2. The formula for the optimal coordinate ascent variational update (equation 7) should be the exponentiated log of the joint distribution [1] and will be update to,\n\n$$\n\\begin{align}\nq(x) \\propto \\exp \\left( \\mathbb{E}_{q(c, \\gamma)} \\left[ \\log p(x, y, c, \\gamma | \\theta)  \\right] \\right). \\tag{2}\n\\end{align}$$\n\nWe proceed by substituting equation 1 above into equation 2, dropping out constant terms,\n\n$$\\\\mathbb{E}\\_{q(c, \\\\gamma)} \\\\left[ \\\\log p(x, y, c, \\\\gamma | \\\\theta)  \\\\right] = \\\\mathbb{E}\\_{q(c, \\\\gamma)} \\\\left[   \\\\log p(x_1) + \\sum_{t=1}^ T \\\\log p(y_t | x_t)  + \\\\sum_{t=1} ^ {T-1} \\\\log p(x_{t+1} | x_{t}, c_{t}) \\\\right] + {\\\\rm const.} \\tag{3}$$\n\n\nNext, our assumed decomposition in line 157 gives us that $p(x_{t+1} | x_{t}, c_t) = p(x_{t+1} = l_{t+1} + b_{t+1})$. The distribution over $x_{t+1}$ is computed as,\n\n$$\n\\begin{align}\n    p(x_{t+1}=l_{t+1} + b_{t+1}) &= p(l_{t+1} | l_t, c_t) \\star p(b_{t+1} | b_{t})  \\\\\\\\ \n    &= N(x_{t+1} ; l_t + F_t l_t + b_t, \\Sigma_x) \\tag{5}\n\\end{align}\n$$\n\nwhere we have defined the convolution operator $\\star$ and $\\Sigma_x =\\Sigma_b + \\Sigma_l$. Substituting the above equations 5 and 3 into the above equation 2 gives us the optimal variational update,\n\n$$q(x) \\propto N(\\mu_1, \\Sigma_1) \\prod_{t=1}^T N(y_t; D x_t, \\Sigma_y)  \\prod_{t=1}^{T-1} N(x_{t+1} ; l_t + \\hat{F}_t l_t + b_t, \\Sigma_x)\n$$\n\nwhere $\\hat{F}\\_t$ is estimated from samples from $\\hat{c}\\_{t} \\sim q(c,\\gamma)$. We note that the additional emissions term in the above equation 1 leads to an additional term in our above equation 8. \n\nAgain, we thank the reviewer for bringing this to our attention. We plan on correcting the typos and will include this discussion into the camera-ready version if accepted.\n\n[1] Blei et al. \"Variational inference: A review for statisticians.\" 2017.", "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": "mEIU4AljqH", "reviewer": "Reviewer_M5AA", "summary": "The paper proposes the probabilistic decomposed linear dynamical systems (p-dLDS) model, which extends the existing dLDS model. With the time-varying offset term and probabilistic formulation, p-dLDS improves upon the dLDS model in terms of robustness to noise and the ability to capture systems that orbit multiple fixed points. The authors show the advantage of the p-dLDS model over SLDS and dLDS through simulated and real datasets.", "questions": "- When does the model break due to the inclusion of time-varying offset terms? In other words, are there cases or hyperparameter settings when the model fails? How should the users avoid this? Are the covariances in line 158 learned or user-specified?\n\n- The paper resorts to using PCA to determine the dimensionality of the model. Is there a reason why the dimensionality is not chosen via cross-validation?\n\n- Minor: In Section E.1 Figure 6 panel A, it seems like the ELBO peaks at the first few iterations, then drops, and finally converges to a value smaller than the peak. I wonder if there's an explanation for this.\n\n- How many hyperparameters are there for p-dLDS? How are these chosen, and how simple or complex is the model selection process for p-dLDS compared to the hyperparameter sweep for rSLDS?", "limitations": "- As the authors noted, one limitation of p-dLDS is that it assumes Gaussian observation noise.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- Clarity: The paper is presented clearly, with easy-to-understand figures and descriptions for the model, inference procedure, and experimental results. The background and related work sections are also well-written, making the paper easy to understand for readers who are not experts in the field.\n\n- Reproducibility: The paper includes code and hyperparameter settings to reproduce its experiments, which is crucial for the ML community.\n\n- Significance: The probabilistic extension, in addition to the time-varying offset term, adds value to the recently developed dLDS, which suffers from dynamics noise.", "weaknesses": "- Experiments\nAlthough existing experimental results are exciting, I have a couple of suggestions. \n\nOne is that, while the paper claims the robustness of p-dLDS to dynamics noise, it does not have experiments that show how robust p-dLDS is to noise. In other words, as you sweep the experiments with different dynamics noises, when do SLDS, dLDS, and p-dLDS show similar/different results? When does p-dLDS fail? Showing the advantage of p-dLDS over other models in a realistic setting of dynamics noise would be important.\n\nIn addition, the experiments set the number of DOs and discrete states the same for p-dLDS and rSLDS (e.g., Section 4.2). For a fair comparison, I think that the number of discrete states for rSLDS should be set separately on its own via cross-validation. It could be that rSLDS needs more discrete states to perform as well as p-dLDS on e.g., linear classification of reach directions.\n\n- Interpretation of discrete states of p-dLDS\nTo my understanding, the paper does not explain how the learned p-dLDS parameters are used to segment the data to infer \"discrete\" states, as in Figures 2 and 4. In addition, for the NASCAR experiment, p-dLDS leads to fewer segmentations than four true segments. I am unsure whether we could say that this representation is more \"parsimonious\" while it seems like it leads to incorrect segmentations. In the Lorenz attractor experiment, it also seems like the rSLDS segmentation is more interpretable than p-dLDS segmentation.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful reading of our submission and are encouraged that you find our work to be \"well-written\" and \"adds value to the recently developed dLDS\".\n\n**W1. Sweeping dynamics noise. Advantage of p-dLDS in a realistic setting.**\n\nWe agree that it's important to understand our model's behavior under different noise conditions. In our experiments, we focus on covering a range of neurally plausible sources of noise and demonstrate that p-dLDS improves significantly over existing models on a variety of metrics. Moreover, we emphasize that our clinical neurophysiological experiment showcases exactly “the advantage of p-dLDS over other models in a realistic setting”. We use real data with real noise (e.g. noisy trial structure, unexpected head movements, etc.), and demonstrate that p-dLDS can capture meaningful structure where other models cannot.\n\nAs discussed in our introduction, dynamical noise can arise from numerous sources, making it challenging to thoroughly explore all scenarios within the scope of this paper.  Without a specific end application in mind, it is difficult to construct a experiment that fully addresses this point. \n\n\n**W2. Selecting number of discrete states for rSLDS via cross validation.**\n\nGreat point! In our original experiments, we focus on approximately matching the model parameters to control for the model complexity. However, we include an additional experiment where we sweep the number of discrete states for rSLDS (see additional pdf). While more discrete states did improve rSLDS performance, it still underperforms with respect to p-dLDS at a lower parameter count, highlighting the inherent limitations of using a switched formulation to model a continuum of signals. \n\n**W3. How are \"discrete\" states defined for pdLDS? Parsimonious and interpretable segmentations.**\n\nAs noted in the captions of Figures 2B and 2D, “discrete” states for dLDS models are colored according to the “dominant coefficient” which we define as the dynamic operator (DO) state with the largest coefficient magnitude. We plan on making our definition more clear in the main text in the camera-ready version of the paper, if accepted.\n\nFor NASCAR, we respectfully disagree that pdLDS \"leads to incorrect segmentations\". Due to the pdLDS model formulation, the learned model components are capturing the notion that different track segments can have their movement captured by the same DO (i.e., using the same movement but backward). While this results in fewer segments than the ground truth switched model the data was generated with, it correctly identifies redundancies in the switched segments that may be important to capture in many applications. The ability to identify these redundancies is a key advantage of pdLDS over switched models. While you certainly may want to distinguish these cases in a specific application, this is easily done by looking at the sign of the coefficients. We emphasize here again that pdLDS recovers the track segments despite noise (unlike other methods), leading to a variety of quantitative performance improvements.\n\n\nIn the Lorenz experiment, while rSLDS segmentations may seem simpler and more interpretable, it inappropriately under-segments the true system by assigning the same state to both fast and slow segments, which loses important dynamics information about speed. This behavior can be problematic for signals that vary across a continuum, such as those from the reaching experiment. Moreover, our quantitative assessment (Table 1) shows that p-dLDS's representation leads to improved performance across various metrics, supporting its value despite apparent complexity.\n\n**Q1:** While we agree that understanding the limitations of the time-varying offset term is important, fully characterizing the model's failure modes without a specific application in mind is challenging. There are numerous ways to study when the model breaks due to the offset terms, and the most relevant ones depend on the specific use case. \n\nRegarding the covariances, lemma 2 allows users to implicitly specify the sum of both covariances $\\Sigma_l + \\Sigma_b$ through the selection of the smoothing window hyperparameter $S$. Given a particular $S$, our fitting procedure estimates the value of the covariances $\\Sigma_l + \\Sigma_b$. \n\n**Q2:** We use PCA to determine the size of the latent space due to its simplicity, computational efficiency, and ability to control for parameter count across different models. In computational neuroscience, PCA is a widely used and well established method that has provided numerous scientific insights [citations available upon request]. In contrast, cross-validation can be costly and may select different latent space sizes for different models, which can confound differences due to parameter count with differences in the latent dynamics model. By selecting the latent dimension across all models using PCA, we control for model complexity and can more directly compare their learned representations.\n\n**Q3:**  We believe this results from the training dynamics introduced by using SGD and Momentum, combined with the challenging non-convex optimization landscape defined by our ELBO objective. Generally speaking, variational inference is only guaranteed to converge to a local optimum.\n\n**Q4:** Hyperparameters (HPs) include $M$, $N$, and $K$ for model dimensionality; $\\xi$ and $S$ for coefficient inference; and any SGD HPs (e.g., learning rate, momentum). In general, HPs can be estimated using either domain knowledge or cross-validation. While p-dLDS has a more complex parameter sweep than rSLDS, it is easier to fit than dLDS since it does not require the manual balancing of the BPDN-DF Lagrange multipliers. This greatly simplifies the fitting process for decomposed models.", "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": "4abE2lZI2q", "reviewer": "Reviewer_3hfN", "summary": "The Authors propose a probabilistic extension to the \"Decomposed Linear Dynamical Systems (dDLS)\" method by Mudrik et al. (2024). The proposed p-dLDS belongs to a family of models specifically designed to describe neural activity from high-dimensional dynamical data. The effectiveness of p-dLDS is evaluated on a set of synthetic benchmarks and a real dataset. \n\nAs the primary use-case of p-dLDS -- modeling neural data -- is not my expertise, my evaluation capability for this paper is quite limited, as reflected in my confidence assessment.", "questions": "- What is the running time of the different baselines? What is the computational complexity of p-dLDS?\n- From the experimental evaluation, I can see that dLDS and p-dLDS behave quite differently. Have you observed any relationship between the posterior mean of the coefficients $\\mathbf{c}_{t}$ in p-dLDS and the deterministic analog in dLDS?", "limitations": "N / A", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "On top of giving accurate predictions of future states, the proposed method recovers a sensible switching behavior for every one of the examples reported.", "weaknesses": "- The experimental results do not report the variance across independent runs. As most of the datasets are synthetic, these should be fairly straightforward to produce. \n- Apart from adding a probabilistic structure to the sparse coefficients of dLDS, p-dLDS adds a slow-fast decomposition of the latent state (Eq. 4). Yet, no ablation study assessing the importance of this modeling choice is reported.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer 3hfN for taking their time to review our work.\n\n**W1:** We report the standard deviation across independent runs in Tables 3, 4 and 5 in the supplementary materials. We omit these values in the main text due to limited space.\n\n**W2:** While we do not perform an ablation study for the offset term, we discuss theoretically in section 3.1 why the lack of an offset term in dLDS leads to a representation that does not generalize in systems that contain multiple fixed points. Namely, that any parameter setting in the original dLDS dynamics model reduces to a linear dynamical system (LDS) which is characterized by a single fixed point centered around the origin. The Lorenz experiment (section 4.1) illustrates how this limitation leads to a representation that does not align with the attractor lobes, inappropriately segments the system radially with respect to the origin (Fig. 2D), and leads to decreased performance in our quantitative metrics (Table 1). In contrast, we show that the offset term in p-dLDS alleviates these problems and enables the learning of a representation that aligns with the multiple fixed points of the Lorenz attractor.\n\n\n**Q1:** We include additional runtime experiments that sweep across a range of dictionary sizes and latent dimensions (see attached pdf). We observe that in general, p-dLDS is significantly faster than dLDS, and even matches rSLDS and SLDS in certain parameter settings. Our approach adopts the worst-case computational complexity of the SBL inference procedure which is $\\mathcal{O}(TK^3)$ in time and $\\mathcal{O}(TK^2)$ in memory due to a matrix inversion required to compute the posterior coefficient covariance [1]. However, as noted by [2], the introduction of an dynamically informed hyperprior reduces the number of overall iterations required for convergence and typically reduces the actual runtime by an order of magnitude. \n\n\n[1] Tipping. \"Sparse Bayesian learning and the relevance vector machine.\" 2001.\n\n[2] O’Shaughnessy et al. \"Sparse Bayesian learning with dynamic filtering for inference of time-varying sparse signals.\" 2019.\n\n\n\n**Q2:** In general, we observe that the posterior mean of the coefficients in p-dLDS will be different from the deterministic analog in dLDS. In Figures 2B, 2D, 4C and 4F, we illustrate these differences by plotting the inferred coefficients from both models side by side. We find that without the probabilistic inference procedure and the time-varying offset term, the inferred coefficients of dLDS can produce switching patterns that may not align with the ground truth (Figures 2B and 2D) and can oscillate unpredictably (Figures 4C and 4F).", "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": "Va05mOGAvU", "reviewer": "Reviewer_mAN7", "summary": "This paper develops a novel model of latent neural dynamics that builds on existing work in switching and time-varying linear state space models. In particular, the authors propose a full probabilistic formulation of decomposed LDS models and provide a variational EM algorithm to do inference in these models. They apply their method to various synthetic and real neural datasets, showing improved performance over the primary time-varying linear alternative methods.", "questions": "- The approach for determining `b_t` via a moving average accomplishes the goals of learning a changing fixed point. However, it appears to be suboptimal for cases where the fixed point can change both rapidly or slowly (i.e. across different timescales) because of the fixed window size. Additionally, the requirement to search over a wide number of window sizes pose schallenges for optimization. I'd encourage the authors to consider other approaches. Could the bias term also be determined via a sum over a dictionary with shared coefficients `c`? \n\n- Relatedly, can the authors provide more details of the hyperparameter selection scheme for the window size? The paper states that for the Lorenz experiment, the window size was sampled uniformly from 2 to the length of the time series. Does the procedure follow the dLDS hyperparameter selection, in which 1000 random choices of `S` and `\\xi` were generated and a separate model was fit given each setting of hyperparameters? What was the optimal offset values on the real datasets, and does this relate to the timescales of the p-dLDS coefficient switching? \n\n- The parameterization and update step for the sparsity coefficients appears to work well on the simulated and real datasets. However, it does require a multistep approximate procedure for the update. I'm wondering if the authors have considered simplifications to the model that may allow for using simpler inference in this step while still achieving the goals of learning sparse coefficients? For example, does it work to model each $\\gamma_t$ independently? This appears to be very similar to the original SBL-DF algorithm and the first step of the sparsity coefficient update which initializes $q(c) q(\\gamma)$ using SBL-DF. How much improved performance does the proposed model & methods have relative to this baseline?\n\n- In the reaching experiment, the reach angles are classified using the discrete states or p-dLDS coefficients and not the continuous states. However, I would expect classification using the continuous states to outperform this in both cases. I would suggest the authors further justify that comparison in the text and compare to classification based on the continuous states.\n\nMinor comments\n\n- It appears that eq 5 describes the joint distribution $p(c_t, \\gamma_t | c_{t-1})$, not just $p(c_t \\mid c_{t-1})$.", "limitations": "The authors address the limitations in the text.", "rating": 7, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Overall, the proposed model and method are original and are a positive contribution to the modeling of neural dynamics. The authors correctly point out shortcomings of both the SLDS and dLDS approaches, and develop a new model that addresses shortcomings. In particular, the probabilistic formulation of dLDS dramatically improves the models performance, and the flexibility given by the time-varying linear combination of DOs helps p-dLDS perform well in situations that the SLDS is not well suited for. The evaluation is high-quality, as the paper compares relevant methods in four different tasks. Finally, the approach is well-motivated, and more often than not the methods and model are clearly described.", "weaknesses": "The model and fitting method have two added complexities that may make it harder than needed to reliably fit the model parameters and make inferences about the data. In particular, the update step sparsity coefficients requires a multi-step approximate procedure and the time-varying offset term depends on a hyperparameter that can take on a wide number of values. These points are discussed in more detail in the questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer mAN7 for their thoughtful feedback and suggested improvements to our experiments. We are encouraged that they found our work to be \"high-quality\" and a \"positive contribution to the modeling of neural dynamics\".\n\n**W1:** While p-dLDS does require more complex mathematical concepts, the main goal of this work is to develop a procedure that is easier and more reliable to work with than dLDS. Specifically, p-dLDS actually simplifies the fitting process by reducing the total number of hyperparameters for inference to two: the SBL-DF tradeoff $\\xi$ and the offset window size $S$. In contrast, dLDS needs three hyperparameters ($\\lambda_0$, $\\lambda_1$, and $\\lambda_2$) for determining the tradeoff between dynamic reconstruction, coefficient sparsity, and coefficient smoothness. In practice, manually balancing these tradeoffs between each of the Lagrange multipliers can be difficult and time-consuming. p-dLDS eliminates the need for manual tuning of these Lagrange multipliers through probabilistic inference. We estimate $\\lambda_0$ with the inferred covariance $\\Sigma_x$, $\\lambda_1$ with $\\gamma_{t,k}$'s from SBL-DF, and $\\lambda_2$ with $\\sigma^2_{t-1,k}$ also from SBL-DF. In practice, we find that our model is much easier to tune compared to the original dLDS model. As demonstrated in the results presented here, the outcomes are more accurate and robust for p-dLDS than dLDS. \n\n**Q1:** We agree that our approach for the offset may be limited when the fixed point changes across different timescales. We did consider a cost-based dictionary learning approach for the offset term early on in our work, but found it challenging to prevent the learned dynamics from collapsing to the degenerate scenario described in Lemma 1. Often, our offset dictionary would span the same subspace as the direction of the dynamics, preventing the learning of meaningful structure in the DOs. \nMoreover, this approach introduces additional complexity with multiple hyperparameters (e.g. for dictionary size, coefficient structure, and more) and requires an expensive iterative solver to infer offset coefficients, further slowing down training. In contrast, our moving average approach is simple, requiring only a single hyperparameter (the window size), and can be efficiently computed in a single pass. \n\n**Q2:** Yes, the p-dLDS hyperparameter selection procedure is very similar to the dLDS one. We use random search with a budget of 1000 samples to determine the values of $S$ and $\\xi$ and fit a separate model for each set of hyperparameters. In our search, we uniformly sample an integer from 2 to the length of the time series $T$. For the real dataset, the optimal offset is $S=76$ which is smaller than the timescale of the p-dLDS coefficient switching (around 150 time points). This suggests that the same DO dynamics may persist even as the fixed point of the system fluctuates throughout the experiment. We appreciate the reviewer's suggestion for clarification and will include these additional details about the hyperparameters in the supplementary materials and the camera-ready version if accepted.\n\n**Q3:** Great question! While we did not directly explore modeling $\\gamma_t$'s independently, previous works [1,2,3] have shown that the introduction of a dynamically informed prior leads to lower error and faster convergence when compared to their independent counterparts in both probabilistic and cost-based inference approaches. In particular, Figures 8 and 10 in [1] demonstrates that SBL-DF (with dynamically informed $\\gamma_t$'s) shows an improvement in error and convergence time over the static SBL (with independent $\\gamma_t$'s). Additionally, we highlight that our inference procedure is computationally efficient, and significantly reduces the training time required for decomposed models (see additional pdf for runtime experiments).\n\n[1] O’Shaughnessy et al. \"Sparse Bayesian learning with dynamic filtering for inference of time-varying sparse signals.\" 2019.\n\n[2] Charles et al. \"Dynamic Filtering of Time-Varying Sparse Signals via $\\ell _1 $ Minimization.\" 2016.\n\n[3] Charles et al. \"Convergence of basis pursuit de-noising with dynamic filtering.\" 2014.\n\n**Q4:**  Thank you for your suggestions! We originally did not consider using the continuous latent state to classify because we wanted to show that the learned dynamics contained information about the reach angle. However, following your suggestion, we found that the classification accuracy based on the continuous latent state $x_t$ outperformed the classification accuracy from previous experiments in all models. Similar to our previous classification set up, features are generated from the continuous latent states by averaging $x_t$'s over all points in time. We highlight that classification based on the continuous states improved p-dLDS's top-1 accuracy to 57.60\\% and top-3 accuracy to 94.87\\% (see additional pdf). \n\n**Minor.** Thank you for your careful review! We will update these typos.", "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:
Thank you to all four reviewers for their insightful feedback. We respond to reviewers individually below, but provide a description of our additional experiments here. Please see our attached pdf for figures relevant to the experiments below. 

**Experiments:**

1. **Increasing the number of rSLDS discrete states in the Reach Experiment:** Reviewer M5AA suggests that the optimal number of rSLDS discrete states may differ from the one used. We sweep the number of discrete states $K$ from 4 to 20 with a step size of 2 and compare against the p-dLDS results at $K=4$. We find that although increasing the number of discrete states leads to improvements in rSLDS performance, it does not outperform p-dLDS performance at $K=4$.  
    

2. **Training Runtime:** Reviewer 3hfN recommends that we compare the runtimes across different models. We report the time per training iteration when sweeping across the dictionary size $K$ and the latent dimensions $N$ separately from 5 to 50 with a step size of 5. For all experiments, we report training times on a single time series with length $T=1000$ and observation dimension $M=100$. When sweeping $K$, we fix $N=10$. Similarly when sweeping $N$, we fix $K=10$. We see that the p-dLDS inference procedure significantly reduces the training time required for decomposed models.
    

3. **Reach Classification using Continuous Latent States:** Following Reviewer mAN7's suggestion, we compared classification accuracy based on the dynamics ($z$ or $c$) to that based on the continuous states ($x$). We computed classifier features from $x_{1:T}$ by taking the average over time for each trajectory, similar to section D.5. This improved the Top-1 and Top-3 accuracies for all models, indicating that location in latent state space can be highly informative about specific reach directions. In this experiment, p-dLDS achieved Top-1 and Top-3 accuracies of 57.69\% and 94.87\% respectively.

Meta-review: This is a nice extension of recent work on decomposed LDS models. The novel contribution is a probabilistic formulation of the model and a variational EM framework for state and parameter estimation. I was particularly impressed with the improvements in state segmentation on real datasets. Please incorporate the changes and clarifications described during the discussion period. I think this will be a solid contribution to NeurIPS.
Decision: Accept (poster)