Learning Mixtures of Linear Dynamical Systems via Hybrid Tensor-EM Method
Authors: Lulu Gong
Shreya Saxena
Venue: ICLR 2026 Poster (2026)
Keywords: Mixture of linear dynamical systems
Tensor-based moment method
Expectation-Maximization
Latent dynamical systems
Neural data analysis
Primary area: learning on time series and dynamical systems
TL;DR: We propose a novel approach combining tensor-based moments and EM refinement for learning mixtures of linear dynamical systems, which enables reliable and improved recovery of latent systems and is validated on synthetic and neural data.

Abstract: Linear dynamical systems (LDSs) have been powerful tools for modeling high-dimensional time-series data across many domains, including neuroscience. However, a single LDS often struggles to capture the heterogeneity of neural data, where trajectories recorded under different conditions can have variations in their dynamics. Mixtures of linear dynamical systems (MoLDS) provide a path to model these variations in temporal dynamics for different observed trajectories. 
However, MoLDS remains difficult to apply in complex and noisy settings, limiting its practical use in neural data analysis. Tensor-based moment methods can provide global identifiability guarantees for MoLDS, but their practical performance degrades in high-noise or complex scenarios. Commonly used expectation-maximization (EM) methods offer flexibility in fitting latent models but are highly sensitive to initialization and prone to poor local minima. Here, we propose a tensor-based moment method that provides identifiability guarantees for learning MoLDS, which can be followed by EM updates to combine the strengths of both approaches. The novelty in our approach lies in the construction of moment tensors using the input-output data, on which we then apply Simultaneous Matrix Diagonalization (SMD) to recover globally consistent estimates of mixture weights and system parameters. These estimates can then be refined through a full Kalman EM algorithm, with closed-form updates for all LDS parameters. We validate our framework on synthetic benchmarks and real-world datasets. On synthetic data, the proposed Tensor-EM method achieves more reliable recovery and improved robustness compared to either pure tensor or randomly initialized EM methods.  
We then apply this method to two neural datasets from non-human primates doing reaching tasks. For both datasets, our method successfully
models and clusters different conditions as separate subsystems.
These results demonstrate that MoLDS provides an effective framework for modeling complex neural data in different brain regions, and that Tensor-EM is a principled and reliable approach to MoLDS learning for these applications.

Published as a conference paper at ICLR 2026
LEARNING MIXTURES OF LINEAR DYNAMICAL SYS- TEMS VIA HYBRID TENSOR-EM METHOD
Lulu Gong and Shreya Saxena
Department of Biomedical Engineering Wu Tsai Institute, Center for Neurocomputation and Machine Intelligence Yale University, New Haven, CT, 06510, USA {lulu. gong) gongeoyalee......o..ooo. @yale.
ABSTRACT
Linear dynamical systems (LDSs) have been powerful tools for modeling high- dimensional time-series data across many domains, including neuroscience. How- ever, single LDS often struggles to capture the heterogeneity of neural data, where trajectories recorded under different conditions can have variations in their dynamics. Mixtures of linear dynamical systems (MoLDS) provide a path to model these variations in temporal dynamics for different observed trajectories. However, MoLDS remains difficult to apply in complex and noisy settings, limiting its practi- cal use in neural data analysis. Tensor-based moment methods can provide global identifiability guarantees for MoLDS, but their practical performance degrades in high-noise or complex scenarios. Commonly used expectation masiiinzaaiiiiiiiinn (EM) methods offer flexibility in fitting latent models but are highly sensitive initialization and prone to poor local minima. Here, we propose tensor-based moment method that provides identifiability guarantees for learning MoL.DS, which can be followed by EM updates to combine the strengths of both approaches. The novelty in our approach lies the construction of moment tensors using the input- output data, which we then apply Simultaneous Matrix Diagonalization (SMD) to recover globally consistent estimates of mixture weights and system parame- ters. These estimates can then be refined through full Kalman EM algorithm, with closed-form updates for all LDS parameters. We validate our framework on synthetic benchmarks and real-world datasets. On synthetic data, the proposed Tensor-FM method achieves more reliable recovery and improved robustness com- pared to either pure tensor or randomly initialized EM methods. We then apply this method to two neural datasets from non-human primates doing reaching tasks. For both datasets, our method successfully models and clusters different conditions as separate subsystems. These results demonstrate that MoLDS provides an effective framework for modeling complex neural data different brain regions, and that Tensor-EM is principled and reliable approach to MoLDS learning for these applications.
INTRODUCTION
Neuroscience experiments now produce large volumes of high-dimensional time-series datasets, calling for new computational tools to uncover the dynamical principles underlying brain function (Paninski & Cunningham] 2018 Stringer Pachitariu 2024 Urai et al. 2022 Saxena & Cunningham 2019). These recorded data often originate from multiple, distinct underlying dynamical processes, yet the identity of the generating systems at any given time unknown. Estimating and recovering the parameters and dynamics of such latent systems from these mixed trajectories is central challenge in system identification and machine learning (Ljung 1998 Durbin & Koopman 2012 Bishop 2006).
Classical mixture models, such as mixtures of Gaussians (MoG) (Dempster et al. 1977 and mixtures of linear regressions (MLR) (De Veaux 1989). provide valuable tools for modeling these heteroge- neous data. However, they primarily designed for static settings and do not explicitly capture temporal dependencies. limiting their applicability to sequential data where temporal dynamics are central. In contrast, dynamical models such as linear dynamical systems (LDS) and their modern

---

Published as a conference paper at ICLR 2026
extensions, such as the (recurrent) switching LDS (SLDS) and the decomposed LDS (dLDS) (Gahraa
mani & Hinton 2000 Fox 2009: Linderman et al. 2017: 7 Mudrik et al. 2024: Zhang & Saxena 2024), are suitable for tumesceres data modeling with latent states and potential regime switches.
Switching models typically target a single long trajectory and are prone to solutions that contain frequent switching; they not suitable when the goal is to learn parameters of multiple distinct
LDSs from independent trials, which is precisely the common structure in neural experiments.
Mixtures of linear dynamical systems (MoLDS) (Chen & Poor 2022 Bakshi al., [023 Rui
& Dahleh 2025) effectively address this setting by treating each trajectory as originating from a single latent LDS. Then the learning of MoLDS aims to uncover the number of subsystems. their parameters, and the mixture weights from collections of input-output trajectories. This formulation enables direct parameter recovery and interpretability, making it appealing for many applications, including large-scale neural data analysis. Neural datasets usually have large numbers of recordings across behavioral or task conditions, with many trials under each condition. An important question is whether neural populations reuse common latent dynamics across conditions or exhibit distinct dynamics for different behaviors (Vyas et al. 2020 Athalye et al. 2023). MoLDS can provide a latent dynamics, it can reveal the mixture structure of all trials in the dataset and offer a comprehensive view of the underlying dynamical structures and motifs.
principled way explore this question: by identifying shared and
For the inference of MoLDS, tensor-based moment methods are commonly employed, where high- order statistical moments of the input-output data reorganized into structured tensors, and their decomposition provides estimates of mixture weights and LDS parameters. The appeal of these algebraic approaches lies in their global identifiability: unlike iterative optimization methods that navigate non-convex landscapes, tensor decomposition exploits the algebraic structure of moments to directly recover parameters through polynomial equation systems that admit unique solutions under ideal conditions (Anandkumar et al.) 2014). However, their practical performance is often limited because moment estimates become imperfect in realistic, noisy datasets, leading to de- graded parameter recovery (Kuleshov et al. 2015). In paraliel, likeiihood-base approaches such as expectation-maxtt mantttionmmaiiizattionn (EM) have long been widely used for fitting classical mixture models and LDSs. In this context, EM provides a flexible iterative procedure for jointly estimating latent states, mixture weights, and system parameters through local likelihood optimization. While powerful and widely adopted, EM suffers from well-known sensitivity to initialization and susceptibility to poor local minima (Xu & Jordan [996] Bishopp 2006). These limitations become particularly problematic in the MoLDS setting where the parameter space is high-dimensional and the likelihood surface is highly multimodal.
Here, we propose hybrid Tensor-EM framework that strategically combines global initialization with local refinement of these two methods. We first apply tensor decomposition based on Simultancous Matrix Diagonalization (SMD) (Kuleshov et al. 2015) to obtain stable and accurate initial estimates. We then use these estimates to initialize a full Kalman filter-smoother EM procedure, which refines all parameters through closed-fomm updates over all trajectories. This hybrid approach harnesses the global identifiability of tensor methods for robust initialization, while leveraging EM's superior local optimization, to achieve both reliability and efficiency.
We validate this framework on both synthetic and real-world neural datasets. On synthetic bench- marks, the proposed Tensor-EM method achieves more reliable recovery and improved robustness compared to (i) pure tensor methods and (ii) EM with random initialization. Next, we analyze neural recordings from two different experiments: (1) Recordings from monkey somatosensory cortex during center-out reaches in different directions, where Tensor-EM identifies distinct dynamical clusters corresponding to the reaching directions, matching supervised LDS fits per direction but achieved in fully unsupervised manner; (2) Recordings from the dorsal premotor cortex while monkey performs sequential reaches in more distributed directions, where Tensor-EM succeeds in parsing the different trials into direction-specifi dynamical models. These results establish MoLDS as an effective framework for modeling heterogeneous neural systems, and demonstrate that Tensor-EM provides principled and reliable solution for learning MoLDS in both synthetic and challenging real-world settings.
2

---

Published as a conference paper at ICLR 2026
2
RELATED WORK
Mixtures models. Classical mixture models (c.g. MoG and MLR) capture heterogeneity but not explicit temporal structure (Dempster et al. 1977 De Veaux 1989 Li Liang 2018). Importantly, MoLDS is related to MLR through lagged- input representations: by augmenting inputs with their past values, an MLR model can approximate certain temporal dependencies (Rui Dahleh (025). Through this connection, MoLDS inherits useful algorithmic tools, including spectral tensor methods
and optimization approaches (Anandkumar et al. 2014 Yi et al. 2016 Pal et al. 2022 Li & Liang 2018). while maintaining its superior modeling capacity for dynamical systems.
LDS models and their variants. LDS models have been widely used to model time-series data, and several extensions have been developed to better handle nonstationarity and nonlinear structures. SLDS methods (Ghahramani & Hinton 2000; Linderman et al.] (017) model long trajectories that switch between different dynamical regimes over time, requiring joint inference of both continuous latent states and discrete mode sequences. The dLDS method (Mudrik et al. 2024] Chen et al. 2024) captures more gradual or overlapping changes by expressing dynamics sparse, vire-nryiing combinations of basis operators. These frameworks are designed for long, nonstationary sequences where the dynamics themselves evolve over time. In contrast, the MoLDS setting assumes that each short trajectory is well described by single LDS drawn from a collection of LDS components. The goal then is to identify the set of latent dynamical systems and assign trials to those components, rather than modeling intra-trial dynamics changes.
Tensor methods and EM. Tensor decomposition methods offer a principled algebraic approach to parameter estimation in latent variable models, with polynomial-time algorithms and theoretical identifiability guarantees (Anandkumar et al.] 2014 Panagakis et al. 2021). Recent work in tensor- based MoLDS learning has explored different moment construction strategies and decomposition algorithms. Early approaches applied Jennrich's 's algorithm directly to input-output moments (Bakshi al. 2023), while more recent work incorporates temporal lag structure using the MLR reformulation and the robust tensor power method (Rui & Dahleh, (025). Our approach adopts the Simultaneous Matrix Diagonalization (SMD) method (Kuleshov et al.) 2015), which also operates on whitened tensor slices and offers improved numerical stability and robustness to noise, which are critical advantages in the challenging MoLDS setting and neural data.
The combination of tensor initialization followed by EM-style refinement has been proven effective in mixture model settings. the MLR literature, tensor methods have been shown to provide globally consistent initial estimates that lie within the basin of attraction of the maximum likelihood estimator, which are then refined using alternating minimization (Yi et al. 2016 Zhong et al. 2016 Chen et al. 202T). However, alternating minimization represents a simplified version of EM that uses hard assignments rather than the probabilistic responsibilities essential for handling uncertainty in noisy settings. Our work extends this paradigm to the more complex MoLDS setting by combining tensor initialization with more powerful Kalman smoothh EM, including proper handling of latent state inference and closed-form parameter updates.
Contributions. Our work makes both methodological and empirical contributions to MoLDS learning and shows its practical utility in neuroscience settings. Methodologically, relative existing MoL.DS tensor methods, our approach makes several key advances: (i) we employ SMD for more stable decomposition of whitened moment tensors, (ii) we provide principled initialization strategy for noise parameters (Q, R) based on residual covariances, which were missing from prior tensor-based sufficient
approaches, and (iii) we integrate a complete EM procedure with peponnibilitywweighee statistics for all LDS parameter updates. Compared to existing tensor-alet alternating minimization pipelines for MLR, our method leverages the full flexibility of Kalman filtering and smoothing, which is essential for addressing the temporal dependencies and uncertainty quantification required in real MoI.DS applications.
Empirically, we demonstrate the successful applications of Tensor-EM MoLDS methods to com- plex and real-world data analysis. While prior MoLDS tensor work has been limited to synthetic evaluations, we show that our Te framework can effectively analyze neural recordings from different brain regions during distinct reaching tasks and successfully identify distinct dynamical regimes corresponding to different movement directions in fully unsupervised manner. This repre- sents an important step toward making MoLDS a practical tool for important applications, particularly in neuroscience, where capturing heterogeneous dynamics across experimental conditions is central
3

---

Published as a conference paper at ICLR 2026
challenge. Together, these methodological and empirical advances demonstrate that our Tensor- MoLDS framework provides improved robustness and accuracy in both controlled and challenging real-world settings.
MOLDS: MODEL AND TENSOR-EM METHOD
3.1 MIXTURE OF LINEAR DYNAMICAL SYSTEMS (MoLDS)
In the MoLDS setting (Figure D. we observe N input-output trajectories {(u,,0:T)-1, V5,o.T,-)] each generated by of K latent LDS components. Let 2; [K] denote the (unknown) component for trajectory i, drawn i.i.d. as 2; ~ Multinomial(pt.x) where pk E (0, 1) are the mixture weights
with Σ Ph = 1, indicating the probability of a trajectory being generated by component k. Conditional on 2 = k, the data is generated from the following LDS:
1++1 = AATT + Bg + Wt,
N(0, Qk),
(1)
M++1 CL.. + Dr + Des
D4 N0, Rx),
(2)
with A& ERxn E C D2 ER and Qx O, RE > O. The goal of MoLDS
learning is to recover the mixture weights and LDS parameters {Pk, (As, BR. Cs. DR. Qs. R&) H==
These parameters are identifiable only up to two natural ambiguities: the ordering of the components (permutation) and similarity transformations of the latent state realization (which leave the input- output behavior unchanged).
To make this recovery possible, we adopt several standard conditions: (i) inputs are persistently
exciting, (ii) each LDS component is controllable and observable, and (iii) the components are sufficiently separated to ensure identifiability (Bakshi et al. 2023 Rui & Dahleh 2025)
3.2 TENSOR-EM APPROACH OVERVIEW
Our complete approach (Algorithm]) T consists of two stages: a tensor initialization stage (see Algorithm in the Appendix). which provides globally consistent estimates of the mixture weights and system param- eters, and an EM refinement stage (see Algorithm + which further improves these estimates and achieves statistical efficiency. Between these two stages, a key step is the initialization of the noise parameters (Qx, Rx), since these are not identifiable from the censor-based estimates. We address this gap by esti- mating them from the residual covariances computed using the tensor-based parameter estimates (detailed in Appendix D
Maxture of
Linear Dynamical
Systems (MoLDS)
LDS
LDS
Inputs
DD
B
B
Latent
States
A
Az
Neural
Observations
fiving rate
This Tensor-EM approach combines the global iden- tifiability guarantees of algebraic methods with the
statistical optimality of likelihood- This hybrid approach is particularly effective in chal- lenging settings with limited data, high noise, or poor component separation scenarios where neither pure tensor methods nor randomly initialized EM perform reliably.
inference.
Figure 1: Overview of MoLDS and the appli- cation to neural data analysis.
3.3
TENSOR INITIALIZATION FOR MOLDS
The tensor initialization leverages the key insight that MoLDS can be reformulated as MLR through lagged input representations (Rui & Dahleh 2025). This transformation exposes the mixture structure in high-order moments, enabling algebraic recovery of component parameters via tensor decomposi-
tion (see Appendix B for details). The method works by exploiting the impulse-response (Markov
Formally, any invertible matrix M inducing equivalent realization via A - M- M-1A,M Bg +
M Br, Ck + M, and Da unchanged, will yield the same input-output mapping.
4

---

Published as a conference paper at ICLR 2026
Algorithm
Tensor-EM Pipeline for MoLDS
Require:
Trajectories
truncation length L,, L.S order n, #components
K
Ensure: Mixture weights {px and LDS parameters {(Ax, Be, Ckz De, Qu. Re)
Stage 1: Tensor Initialization
(Appendix Alg.
1: Transform MoLDS to MLR via lagged inputs; construct moment tensors M2, Mg
2: Apply whitening and SMD to recover mixture weights {px} and Markov parameters
3: Realize LDS matrices {(Ak, Bg, C,, Dk)3 via Ho- Ho-Kalman algorithm
4: Initialize noise parameters {(QR) R }] from residual covariances
(Appendix (Appendix Alg.
Stage EM Refinement
5: repeat
6: E-step: Compute trajectory responsibiliities "H, via Kalman filter likelihoods
7: Run Kalman smoother to obtain responsibility sepponibblity-weiihteed sufficient statistics
8: M-step: Update mixture weights and all LDS parameters via dlsed-form MLE
9: until convergence in log-likelihood
10: return Refined parameters {px, (A,, BRs Cz2 DR, Qs. R
parameter) representation LDSs. We first construct second- and third-order cross moments of the lagged snput-output data, denoted as M2 and M3
M2 = 2
Σ Ov, Ia), JENS
M3 6|NV
Σ -x()))))()))) JEN
(3)
where No and Ng are two disjoint subsets of samples, normalized lagged input, Ia and are correction terms (see Appendix [B.2] for notations). Then apply the whitening transformation W and decompose the resulting tensor into weights and Markov parameter estimates. At last, we recover LDS state-space parameters through the Ho- Kalman realization (Oymak & Ozay] (219)). The core mathematical insight is that if M2 and M3 appropriately constructed from sub-sampled lagged covariates, the whitened tensor
Wi, is the observed output, vj is the
T Ma(W, W) = Σ k=1
63
Pk a
symmetric and orthogonally decomposable, where ag R K are the whitened representations of the regression vectors (scaled Markov parameters) for each LDS component, and pk the corresponding mixture weights. This decomposition is unique and recovers the mixture components (oN,P&) up to permutation and sign. Importantly, we employ SMD (Kuleshov et al. 2015) for this decomposition step, which is shown to be more stable and accurate empirically in Section 4.1 We then apply the Ho-Kalman realization (Oymak & Ozay 2019) recover the state-space parameters for each component, which are then used as principled initializations for the subsequent refinement stage. The complete procedure is provided in Algorithm 3 in Appendix B.4 together with the
MoLDD-tooMMLR reformulation and tensor construction details in Appendix B
3.4
EM REFINEMENT FOR MOLDS
The tensor initialization provides globally consistent estimates of mixture weights and system matrices, but does not recover the noise parameters (Qk, R2) nor achieve optimal statistical accuracy. We therefore need to refine these estimates using a full Kalman filter-smoother EM algorithm that maximizes the observed-data likelihood.
Our EM formulation extends classical mixture EM to the MoLDS setting by computing trajectory- wise responsibilities via Kalman filter likelihoods, then updating parameters from responsibility- weighted sufficient statistics (see Appendix C and Algorithm 4 for details). In brief, at iteration t,
given current parameters 0(0) responsibilities
((p( A(t) B (t) 6(t (() ((0 f( Hem the E-step computes
(t)
exp( log Pk + log 1ogp(y,,0,,11 / ((
Pe.k
Σk exp log 1/ (logr()
+ log | (()))) u,,0:7;;11
(5)
5

---

multiple met- Published as a conference paper at ICLR 2026
Next, we a Kalman smoother to compute responsibiliy-weighhh sufficient statisties s for each
component. The M-step then updates all parameters via closed-form maximum likelihood estimates (MLE)
B(tt)) O (A B, O
o o R
1) MLE-LDS(S 79.6
(6)
Pk (1+1)
1 N
Σ
E
(7)
4 TENSOR-EM PERFORMANCE ON SYNTHETIC DATA
4.1 SMD-TENSOR METHOD PROVIDES MORE RELIABLE RECOVERY ON SIMULATED MOLDS
We demonstrate that the employed
method reliably recovers parameters in synthetic MoLDS with small number of components, low latent dimensionality, and well-separated dynamics. We empirically compare SMD-Tensor against the RTPM-Tensor baseline across a wide range of such settings. Our results show that SMD-Tensoo consistently outperforms RTPM-Tenso rics. In particular, when the number of mixture components is small and sufficient data are available, SMD achieves near-optimal recovery of both the mixture weights and the Markov parameters.
(a)
Error Number of Trajectories
(b)
Error
Trajectory Langth
(c)
Beet Error Heatmmp
2.5
1200
1000
4.5
1-1249
1.5
8:43
090
Mean
559 1000
200
1200
Number oftnjector
lenth
(d)
Ernor vs Number of Trajecterles
(e)
Error vs Trajectory Longth
(9)
Best Error Heatmap
2.5
Weem Manto eree
1000
0.5
000
10
200
400
3248
20c
Number Urapeton
tagettary kmn
pifference
(9)
Mathod Comparison
(h)
(RTPM - SHD)
SMD: Paramoter Recovery
5ND Mean Dme
KCD
Parameter Recovery
400
X00
0.0
XTTM Moan Ertrr
Conficiint Inde:
Figure 2: Comparison of MMDDTTnsor and CTPM-TenssOn methods for MoLDS: (a,b) Mean Markov parameter estimation errors of the SMD method decrease as the number of trajectories as N and/or T increase; (c) Heatmap of the best trial result across all (N, T) configurations. (d-f) Corresponding results for the RTPM method. (g) Scatter plot comparing mean Markov errors of RTPM vs. SMD across configurations, with SMD outperforming in 91% of cases. (h) Difference in mean Markov errors between RTPM and SMD (positive values indicate SMD performs better). (i) Example recovery for N = T 1280, where SMD recovers both mixture weights and Markov parameters more accurately.
The first row of Figure[]reorrts results of the SMD-Tensor method for a K = 3 mixture model with LDS dimensions n = 2, m = p = 1. The LDS parameters are randomly generated, with eigenvalues of A constrained to lie inside the unit circle. We vary the trajectory length T' and the number of trajectories N for each configuration. We run multiple independent trials, calculating the discrepancy between estimated and true Markov parameters. The second row of Figure shows the result for the -TPM-Tensor method. It is noticed that the mean Markov parameter errors decrease with increasing trajectory length and number of trajectories for both methods, but the trend is more pronounced and

---

Published as a conference paper at ICLR 2026
(a) Relative Accuracy Comparison
(b) Markov Parameter
(c) Weight Errors
belier)
Acc
2.02
Merlov Panemeler Crror
0.3
Weight Eroor
Asprozate
Ace
0.2
"enste
Randam EM
0.8
0.5
ThnsoneM
Randum
Randum
(d) Aggregate Paramater Errors
(a) Iteration Counts
(1) Error vs Efficiency
lacom. mean of A.B.C)
Aggnegatte Param Bror
EM tareioes
Marhov Parameter
4.0
20
0.5
TensonEM
Randem-D4 -34
tervol
TM
Random-DD
0
10
EM Iterations
20
Figure 3: Performance comparison of pure tensor, Tensor-EM, and random-initialized EM on simulated MoLDS: (a) Relative accuracy radar plot across metrics of Markov parameter accuracy, weight accuracy, and aggregate parameter accuracy. (b-c) Tensor-EM achieves the lowest Markov parameter and weight errors, while random EM performs the worst and shows high variability. (d) reports aggregate parameter emrors (geometric mcan of A, B, C' errors). (c) Tensor-EM converges in far fewer iterations than random EM, highlighting its efficiency. (f) The error-efficiency plot shows that Tensor-EM combines low error with moderate iteration cost, yielding robust and accurate recovery.
stable for SMD, as reflected in the heat maps (Figure 2..). Across nearly all (N,T) cases, SMD yields consistently lower errors (Figures (FiguresPg.h). hg. In cases with larger and N, SMD achieves highly accurate recovery of both mixture weights and Markov parameters (Figure 2 underscoring its robustness and reliability for the MoLDS setting.
4.2 TENSOR-EM IMPROVES ROBUSTNESS AND ACCURACY FOR COMPLEX SYNTHETIC MoLDS
In complex MoLDS settings with many components, purely tensor-based methods and randomly initialized EM often fail to achieve accurate recovery, either due to noisy parameter estimates or convergence to poor local optima. Our proposed Tensor-EM approach overcomes these limitations by combining globally consistent tensor initialization with EM-based refinement, resulting in more robust and accurate learning. We demonstrate these advantages in this section.
p = 2) with zero direct feedthrough (D 0). Across metrics including Markov parameter errors, mixture weight errors, aggregate LDS parameter errors, and iteration counts, Tensor-E. consistently outperforms all base- lines. It achieves substantially lower Markov and weight errors (panels b and c), while maintaining an aggregate parameter еигог comparable the pure tensor solution (panel d) At the same time, Tensor-EM converges in much fewer iterations than random EM (panel e), leading to a better error- efficiency tradeoff (panel f). The radar plot Figure ()) summarizes these improvements, showing that Tensor-EM achieves the most balanced performance across all evaluation metrics. These results highlight that Tensor-EM effectively combines the strengths of tensor methods and EM and enables more reliable parameter recovery in synthetic MoLDS settings where standalone tensor and randomly initialized EM approaches often underperform.
Figure 3 presents results simulated K MoLDS (n
"The aggregate parameter error. i.e., the geometric mean of A, B, C' errors, is coarse summary: unlike Markov parameter errors, it may be less precise since A, B, C can differ by similarity transformations without altering the underlying dynamies.
7

---

Published as a conference paper at ICLR 2026
5 TENSOR-EM MoLDS PROVIDES RELIABLE AND INTERPRETABLE RECOVERY ON REAL-WORLD APPLICATIONS
We next apply the proposed MoLDS method to two real-world neural datasets.
5.1 ARFA2 DATASET
We first analyze a neural dataset from a sin- gle recording session of macaque performing center-out reaches, with neural activity recorded from the somatosensory cortex (FigureHaa. Dur- ing the experiment, the monkey used a ma- nipulandum to direct a cursor from the cen- ter of a display to one of eight targets, and neural activity was recorded from Brodmann's area2 of the somatosensory cortex (see (Chowd- hury et al. (020) for experimental details, and (illlrr 2022) for dataset location). The neural recordings consist of 65 single-unit spike times, which converted spike rates. In addition to neural data, the position of the monkey's 's hand, cursor position, force applied to the manipulan- dum, and hand velocity were recorded during the experiment. There are directions in the task, and each direction has multiple trials of trajectories. For the MoLDS fitting, we extract the movement-relas relatte segment of each trial, de- fined as the window from 100 ms before to 500 ms after movement onset. The 65-dimensional observations are reduced to 6-dimension by using the standard PCA, which explains more than 90 percent variance of neural activitics. These PCs are then taken outputs, while the hand velocity variables are taken as inputs. In addition, we also evaluate the dataset with 20 PCs using MoLDS, and consistent results are found (see Figure [T]. which confirms the effectiveness of the Tensor- algorithm to a MoLDS model on this dataset.
Poekay feech Ingecteries
Searon 2000 Rar
Peu Populattion Authity
Figure 4: Neural datasets overview: (a) The Area2 Dataset contains neural trajectories from monkey primary somatosensory cortex; reach tra- jectories and firing rate from multiple neurons. (b) The PMd Dataset contains recordings from mon- key dorsal premotor cortex; hand speed and neu- rons' rasters.
We evaluate the MoLDS with Tensor-EM, Random- -M,, and pure tensor methods using standard train/ seinn/validation/ees split of the dataset (see full pipeline in App. E). For each hyperparameter setting (including K), models are trained on the training set and scored the validation set using negative log- ligliiildhood (NLL), one-step-ahead RMSE, and BIC. The model used for test-time analysis is the one minimizing validation BIC (we also report NLL/RMSE). For trial and component k. we compute a responsibility Fak X exp((()) where e(() is the one-step Kalman log-likelihood under component ks For each movement direction, the dominant component is the arg max; of the mean responsibility across its trials.
Figure 5 summarizes the results. As shown in Figure 5p. the validation criteria consistently favor 3-component MoLDS trained with Tensor- EM. Moreover, the one-step predictions ÿe closely track observations y on the test data as in Figure 5p. In Figure the dominant-component maps reveal three cross-trial trial clusters aligned with directions, and the usage fractions in Figure 5 quantify the prevalence of MoLDS components on the test set.
In this setup, we also compare the Tensor-EM MoLDS method with the supervised learning results of LDS, where we train (1) one LDS on all trials where the monkey reaches specific direction (per-dir LDS), and separately (2) separate LDS fit on each trial regardless of reaching direction (single-trial LDS). Finally, we cluster the parameters of the different LDS's (per-dir or single-trial). See Figure 8 in the Appendix for representation of the clustered parameters and more details. The per-direction LDS baseline closely matches the result of MoLDS with the Tensor- EM method rather than Random-EM in trial groupings (see Figure 5L.), and the impulse responses are highly similar (see Figure 10 in the Appendix). We also train an SLDS in similar unsupervised way as MoLDS (no direction labels) and find that it does not yield meaningful cross-trial clusters here, while it is effective
8

---

Published as a conference paper at ICLR 2026
(A] Validation metrics
(b) Test 1-stop predictions
MIC
127
Bha mentoee The e13 gest
187
E5
Xx]
E4-
181
12
Metrod
Method
Mathee
(c) Per vo MoLDS
SLDS va MOLDS
30
Single-trial LOS
20'
Random TensoNEN
McLDS component usage
335
335
135T
3351
20
282
0.11
Figure 5: MoLDS application on Area2 Dataset: (a) Validation metrics of MoLDS for different K. (b) One-step predictions for two example trials using corresponding LDS components from the MoLDS selected by the lowest validation metrics. (c) Agreement between Tensor-EN MoLDS trial assignments (outer ring) and per-direction LDS, SLDS, single-trial LDS, Random-EM MoLDS clusters (inner rings in each polar plot); the SLDS-based method cannot provide meaningful cross-trial clusters. (d) Tensor-EM MoLDS component usage fractions on held-out test trials.
for possible within-trial regime switches (see Figure 9 in Appendix). These highlight MoLDS's strength in capturing between-trial heterogeneity. In addition, when compared with the Random-EM method, Tensor-EM offers key advantage: the tensor initialization step yields stable starting point that reduces variability across runs, leading to more consistent parameter recovery and recovered model structure as seen Figures 5l-c.
5.2
PMD DATASET
We apply our method recordings from monkey dorsal premotor cortex (PMd) during se- quential reaches (Figure A.). in which trial-wise movement directions are contimuously distributed over the circle. Full experimental details are in (Lawlor et al.. 2018) and the dataset is pro- vided in (Perich et al. (218). The trials' reach direction spans the full polar range rather than clearly separated directions in the Area2 dataset. For the PMd dataset, we similarly extract movement-related activity by taking fixed window (-100 to +500 ms) around the movement onset for each trial. We analyze PMd activity with the 5-dimensional kinematic as inputs, i.e., (x-velocity, y-velocity, x-acceleration, y-acceleration, speed), and the first 16 PCs as outputs. Fol- lowing the Area2 protocol, we train MoLDS across hyperparameters on the PMd training split and select the model by validation criteria (NLL/RMSE/BIC), Figure 64 reports validation results where the 4-component MoLDS trained with the Tensor- EM method has better performance. Figure 6p shows test one-step- ahead reconstructions of the corresponding component from this MoLDS fit.
To accommodate continuously distributed reach angles, we uniformly binned trials into 12 angular bins based on their reach angles. For each bin, assigned dominant component based on the highest mean responsibilityy Figure shows these bins colored by their dominant LDS component and overlays responsibility-weighh we also compare Tensor-EM and Random-EN, MoLDS results. Figure 6 plots impulse-response magnitude across lags of each component, revealing component-speciffe gains and temporal decay; Appendix E.2 (Figure 14) further shows the decomposed responses across input channels. Figure 6 reports Tensor-Ei MoLDS component usage fractions on held-out test trials. Overall, the mixture components of the PMd dataset specialize in different movement directions and exhibit distinct dynamical response profiles.
wreferred-direction (PD) arrows for each component, where
CONCLUSION
MoLDS provides an interpretable framework for modeling repeated, trajectory-level dynamics in heterogeneous data. In this work, proposed hybrid Tensor-EM approach for efficient and

---

Published as a conference paper at ICLR 2026
BMSE
2-114
22"
Componerti
Sast will
(2) Angle bia domimant composentt
4
mg
you
Ireputionn
Tinson-EM
Random-EM
2:22
0214
Comab
Cara
Come)
Figure 6: MoLDS application on PMd Dataset: (a) Validation metrics of MoLDS for different K (NLL, RMSE, BIC). (b) One-step prediction on an example neural trial using the corresponding component from the validation-chosen MoLDS. (c) Angle bins: dominant component; PD arrows sisponsibilityyweiihhee (d) Impulse response magnitude of fensor-EM MoLDS components. (e) Tensor-EM MoLDS component usage fractions on held-out test trials.
reliable MoLDS inference. Specifically, we used SMD for censor-based initialization, followed by Kalman filter-smoother EM refinement. This combination enables accurate parameter recovery and interpretable trajectory clustering in both synthetic benchmarks and real neural datasets. This hybrid framework makes MoLDS practically usable for large-scale, noisy datasets.
Our empirical evaluation focused scenarios with linear dynamics: synthetic data generated from linear systems and neural data consisting of smoothed firing rates over selected movement windows. While these experiments validate the effectiveness of the proposed inference procedure, we have not explored cases with a high likelihood of model mismatch or strong nonlinearities. We envision few extensions that help broaden the applicability of MoLDS. First, nonlinear dynamics may be approximated by segmenting trajectories into shorter sessions that are locally -eel-aaprox by linear dynamical systems. Second, linear components could be augmented with nonlinear features to capture moderate deviations from linear structure. In addition, future directions can include extending the framework to autonomous settings and developing theoretical guarantees under simplified assumptions.
More broadly, MoLDS is well-suited various neuroscience tasks in which multiple experimental conditions share latent dynamical structure but differ in observable trajectories. For example, in sensory decision-making with different stimulus strengths, distinct conditions may reuse a small set of dynamics (Hanks et al. 2015). In context-switching tasks, neural populations in the prefrontal cortex may transition among a number of rule-dependent dynamical modes (Mante et al.) 2013). In motor learning experiments, early, intermediate, and late learning stages may correspond to transitions among a small set of dynamical regimes (Shenoy et al. 2013). By developing a principled Tensor-MM learning and inference framework, this work establishes practical and extensible foundation for applying MoLDS to complex neural datasets. Beyond neuroscience, the framework may also support broader applications. such as identifying clinically relevant behavioral dynamics (Bulteel et al. 2016) or uncovering cell-type-specfii dynamical structure (Luecke et al. 2021).
7 ETHICS STATEMENT
Here, we aim to make methodological and neuroscientific insights, and do not note any negative societal ethical implications.
10

---

Published as a conference paper at ICLR 2026
8 REPRODUCIBILITY STATEMENT
Our work can be reproduced in a straightforward way. The datasets are provided in (Miller 2022 and (Perich al. 2018), with all pre-processee prepprcessiig techniques detailed in (Chowdhury et al., 2020p. (Lawlor et al., 2018), and in the Appendix of this work. Moreover, detailed algorithms and technical details are provided for each step of the learning and inference, with comprehensive pseudo-code for the implementation in the main text and the Appendix.
ACKNOWLEDGMENTS
I.. Gong supported by the Swartz Foundation Postdoctoral Fellowship. S. Saxena supported by the National Science Foundation grant No. 2350329 and the National Institutes of Health grant No. R34DA059718. The authors especially thank Maryann Rui and Dr. Munther Dahleh for helpful discussions.
REFERENCES
Animashree Anandkumar, Rong Ge, Daniel J Hsu, Sham M Kakade, and Matus Telgarsky. Tensor decompositions for learning latent variable models. Journal of Machine Learning Research, 15(1): 2773-2832, 2014.
Vivek R Athalye, Preeya Khanna, Suraj Gowda, Amy L Orsborn, Rui M Costa, and Jose M Carmena. Invariant neural dynamics drive commands control different movements. Current Biology, 33 (14):2962-2976, 2023.
Ainesh Bakshi, Allen Liu, Ankur Moitra, and Morris Yau. Tensor decompositions meet control theory: learning general mixtures of linear dynamical systems. In International Conference on Machine Learning, pp. 1549-1563. PMLR, 2023.
Christopher M Bishop. Pattern recognition and machine learning, volume 4. Springer, 2006.
Kirsten Bulteel, Francis Tuerlinckx, Annette Brose, and Eva Ceulemans. Clustering vector autoregres- sive models: Capturing qualitative differences in within-person dynamics. Frontiers in Psychology, 7:1540, 2016.
Hengnu Chen, Lei Deng, Zheng Qu, Ling Liang, Tianyi Yan, Yuan Xie, and Guoqi Li. Tensor train decomposition for solving large-scale linear equations. Neurocomputing, 464:203-217, 2021.
Yanxi Chen and H Vincent Poor. Learning mixtures linear dynamical systems. In International conference on machine learning, pp. 3507-3557. PMLR, 2022.
Yenho Chen, Noga Mudrik, Kyle A Johnsen, Sankaraleengam Alagapan, Adam S Charles, and Christopher Rozell. Probabilistic decomposed linear dynamical systems for robust discovery of latent neural dynamics. Advances in Neural Information Processing Systems, 37:104443-104470, 2024.
Raeed H Chowdhury, Joshua Glaser, and Lee E Miller. Area of primary somatosensory cortex encodes kinematics of the whole arm. Elife, 9rc48198, 2020.
Richard De Veaux. Mixtures of linear regressions. Computational Statistics & Data Analysis, 8(3): 227-245, 1989.
Arthur P Dempster, Nan M Laird, and Donald Rubin. Maximum likelihood from incomplete data via the em algorithm. Journal of the royal statistical society: series B, 39(1):1-22, 1977.
James Durbin and Siem Jan Koopman. Time series analysis by state space methods. Oxford university press, 2012.
Emily Beth Fox. Bayesian nonparametric learning of complex dynamical phenomena. PhD thesis, Massachusetts Institute of Technology, 2009.
11

---

Published as a conference paper at ICLR 2026
Zoubin Ghahramani and Geoffrey E Hinton. Parameter estimation for linear dynamical systems. 1996. URL (RLLhhts:://// wtt.eeatttttcccca....oooooooaaaa
Zoubin Ghahramani and Geoffrey E Hinton. Variational learning for switching state-space models. Neural computation, 12(4):831-864, 2000.
Timothy D Hanks, Charles D Kopec, Bingni W Brunton, Chunyu Duan, Jeffrey C Erlich, and Carlos D Brody. Distinct relationships of parietal and prefrontal cortices to evidence accumulation. Nature, 2200754 2015.
Volodymyr Kuleshov, Arun Chaganty, and Percy Liang. Tensor factorization via matrix factorization. In Artificial Intelligence and Statistics, pp. 507-516. PMLR, 2015.
Patrick N Lawlor, Matthew G Perich, Lee Miller, and Konrad P Kording, Linear-nonlimear-imee warp-poisson models of neural activity. Journal of Computational Neuroscience, 45(3):173-191, 2018.
Yuanzhi Li and Yingyu Liang. Learning mixtures of linear regressions with nearly optimal complexity. In Conference On Learning Theory, pp. 1144. PMLR, 2018.
Scott Linderman, Matthew Johnson, Andrew Miller, Ryan Adams, David Blei, and Liam Paninski. Bayesian learning and inference in recurrent switching linear dynamical systems. In Artificial intelligence and statistics, pp. 914-922. PMLR, 2017.
Lennart Ljung. System identification. In Signal analysis and prediction, pp. 163-173. Springer, 1998.
Stefanic Luecke, Katherine M Sheu, and Alexander Hoffmaann. Stimulus-specific responses innate immunity: Multilayered regulatory eircuits. Immunity, 14(9):1115-9922, 2021.
Valerio Mante, David Sussillo, Krishna Shenoy, and William T Newsome. Context-dependent
computation by recurrent dynamics in prefrontal cortex. Nature, 503(7474):78 2013.
Lee Miller. Area2 bump: macaque somatosensory area spiking activity during reaching with per- turbations, 2022. URL httts::////::tdiiaatt ottp://ddd.arrii..e...///////e/s orrg///adddseett000007700.0000.0
ve.
0359
Noga Mudrik, Yenho Chen, Eva Yezerets, Christopher J Rozell, and Adam S Charles. Decomposed linear dynamical systems (dlds) for learning the latent components of neural dynamics. Journal of Machine Learning Research, 11559):1-44, 2024.
Samet Oymak and Necmiye Ozay. Non-asymptotic identification of Iti systems from single trajectory. In 2019 American control conference (ACC), pp. 5655-5661. IEEE, 2019.
Soumyabrata Pal, Arya Mazumdar, Rajat Sen, and Avishek Ghosh. On learning mixture of linear regressions in the non-realizable setting. In International Conference on Machine Learning, pp- 17202-17220. PMLR, 2022.
Yannis Panagakis, Jean Kossaifi, Grigorios Chrysos, James Oldfield, Mihalis A Nicolaou, Anima Anandkumar, and Stefanos Zafeiriou. Tensor methods in computer vision and deep learning. Proceedings of the IEEE, 109(5):863-88)) 890, 2021.
Liam Paninski and John P Cunningham. Neural data science: accelerating the experiment-analysis theory cycle in large-scale neuroscience. Current opinion in neurobiology, 50:232-241, 2018.
Matthew G Perich, Patrick N Lawlor, Konrad P Kording, and Lee E Miller. Extracellular neural recordings from macaque primary and dorsal premotor motor cortex during a sequential reaching task, 2018. URL
dataa sets/ wttt:///oooooooooooooooooooott// cortex/pttexteccooo pmd-
Maryann Rui and Munther Dahleh. Finite sample analysis of tensor decomposition for learning mixtures of linear systems. In 7th Annual Learning for Dynamics & Control Conference, pp. 1313-1325. PMLR, 2025.
Shreya Saxena and John P Cunningham. Towards the neural population doctrine, Current opinion in neurobiology, 55: 101-111, 2019.
12

---

Published as a conference paper at ICLR 2026
Krishna V Shenoy, Maneesh Sahani, and Mark M Churchland. Cortical control of arm movements: dynamical systems perspective. Annual review of neuroscience, 36(1): (661):337-359) 2013.
Carsen Stringer and Marius Pachitariu. Analysis methods for large-scale neuronal recordings, Science, 386(6722):7429, 2024.
Anne E Urai, Brent Doiron, Andrew M Leifer, and Anne K Churchland. Large-scale neural recordings call for new insights to link brain and behavior. Nature neuroscience, 25(1):11-19, 2022.
Saurabh Vyas, Matthew D Golub, David Sussillo, and Krishna V Shenoy. Computation through neural population dynamics. Annual review of neuroscience, 43(1):249-211 2020.
Lei Xu and Michael Jordan. On convergence properties of the algorithm for gaussian mixtures.
Neural computation, 8(1):129-151, 1996.
Xinyang Yi, Constantine Caramanis, and Sujay Sanghavi. Solving a mixture of many random linear equations by tensor decomposition and alternating minimization. arXiv preprint crXiv:1608.05749 2016.
Yongxu Zhang and Shreya Saxena. Inference of neural dynamics using switching recurrent neural networks. Advances in Neural Information Processing Systems, 17:131556-13449 2024.
Kai Zhong, Prateek Jain, and Inderjit S Dhillon. Mixed linear regression with multiple components. Advances in neural information processing systems, 29, 2016.
A LLM USAGE STATEMENT
We used LLMs to assist with language editing in parts of the manuscript. All LLM's suggestions were manually verified and edited by the authors.
B TENSOR INITIALIZATION DETAILS
In this section, provide brief introduction MoLDS's connection MLR and the censor-based moment method for MoLDS (see (Bakshi et al.) 2023 Rui & Dahleh, (2025) for more details).
B. MOLDS TO MLR REFORMULATION
The key insight enabling tensor methods for MoLDS is that any linear dynamical system can be equivalently represented through its impulse response (also called Markov parameters), which describes how the system responds unit impulses in its inputs. This representation allows to reformulate MoLDS as MLR, making algebraic moment-basee methods applicable. We now briefly introduce this process following the formulation (Rui Dahleh 2025).
Step 1: Impulse response representation. For notational brevity, we assume an LDS with m-- dimensional inputs and scalar output (p = 1), and zero feedthrough (D₃ = 0). The extension to multiple outputs (p 1) is straightforward. The system parameters are (Ax, BR, Ck). Its impulse response, also called the sequence Markov parameters, defined by
gk(j) = CLAR-1BA B₁ E R",
j == 1,2,3
so that 9k (1) Br 9u (2) = C. Ar. etc.
The sequence (Øx (j)) captures the system's memory: (j) specifies how input applied j steps in the past influences the current output. Given an input sequence {u.), J. the output can be expressed as the infinite convolution
1/+ Σ 9N (i wr + noise terms.
In practice, this sum is truncated after L terms, since (i) decays exponentially for stable systems. denote the stacked first L impulse responses by
(L) [1) ,9xx(2) .... ......99) (I. 4m
g
13

---

Published as a conference paper at ICLR 2026
Step 2: Lagged input construction. To exploit the impulse response representation, construct lagged input vectors that collect recent input history. For trajectory i time 2 L, define
a ( u 11, ... -,,--++1 TR
This vector stacks the most recent L inputs (the current input and the previous L 1), that u has dimension Im and aligns with the stacked impulse responses (L), With this construction, the truncated output becomes
Vi.t Ih (...) + noise terms,
where K; denotes the (unknown) component generating trajectory
With the above input construction, consecutive lagged vectors ūg and ūg,+1 share L-1 input entries, which induces strong statistical dependence and complicates the estimation of higher-order moments. To mitigate this, we sub-sample the time indices with stride L, i.e.,
E { 2L, 3L, ..... T}.
This construction ensures that the resulting lagged vectors non-overlapping. Under standard input assumptions (c.g., i.i.d. or persistently exciting inputs), the sub-sampled vectors are approximately independent, which crucial for consistent moment estimation.
To apply tensor methods, the covariates must have unit variance. We therefore normalize each lagged input by the input standard deviation,
Step 4: Normalization and MLR formulation.
o
2
[[]]u,113]
and flatten the dataset via the re-indexing map (i, -. This yields the mixture of linear regressions (MLR) form
Ú; (U/z Bkg) + ny Egك
where
új Vi.t is the observed output,
uy E L. is the normalized lagged input (covariate).
BR ouge, (L) E B.Lm is the scaled Markov parameter vector (regression coefficient),
ky E {I....,]
K} is the (unknown) component index,
m captures process and observation noise,
E accounts for truncation error from ignoring impulse responses beyond lag L.
Step Mixture structure. Since each trajectory originates from one of latent LDS components with probabilities {px K the regression model inherits the same mixture structure:
Plk-
k] Pk-
Thus the learning task reduces to recovering the mixture weights {pt} and regression vectors {.&,) from the dataset {(., Ø,}}. Once these regression parameters are estimated, the corresponding state-space models (Ax, Bg, Ck) can be reconstructed via the Ho-Kalman realization algorithm
(Ommkk & Ozay 2019).
This reformulation is crucial: it transforms original problem identifying a mixture of L.DSs into the algebraic problem of learning MLR model, for which polynomial-time tensor methods with identifiability and sample-complcxity guarantees are available (Rui & Dahleh 2025).
B.2 MOMENT CONSTRUCTION AND WHITENING
Let Lm denote the dimension of the lagged input vectors. We partition the sample indices into two disjoint subsets N₂ and N3 for constructing second- and third-order moments, respectively as in Ruii & Dahleh 2025). The empirical moments are defined as
Ma
07 JENZ
2 03 Ta),
M3
03 (1 E(v))),
2
6|N3
JENs
14

---

Published as a conference paper at ICLR 2026
where La is the d x d identity and
E(v) Σ (v0,, + 0,08 e, 0e,
with er the r-th standard basis vector in Rd These corrections ensure that the resulting tensors are centered and symmetric.
At the population level, these moments satisfy
K
K
M[
Σ Pk: BR
E[Ma]
Σ
Ph:
1=1
they encode the regression vectors {Ba} and mixture weights {pk)-
To obtain an orthogonally decomposable form, we perform whitening. Let M2 z USUT be the KankKK eigendecomposition, and define
1-1/2 W ((I1:K) 21:K
so that W M2W Tg. Applying W along each mode of M3 yields the whitened tensor
K
r =
M,WW M,W,WW,W) = Σ
83
Og == T BE-
A=1
This tensor symmetric and orthogonally decomposable. We then apply Simultaneous Matrix Diagonalization (SMD) (B.)) Kullehhvv al. 2015) recover {au, Pk Finally, unwhiten to obtain 3g W and recover the state-space parameters (As, Bx, Cş) via Ko-Kalman
realization OOymak Ozay 2019).
B.3 SIMULTANEOUS MATRIX DIAGONALIZATION
This section provides the Simultaneous Matrix Diagonalization (SMD) method for tensor decom- position (Kuleshov et al. 2015) Algorithm 2 SMD recovers tensor components by reducing the problem to joint matrix diagonalization, exploiting linear algebraic structure to recover all components simultaneously rather than sequentially as RTPM.
Jacobi Joint Diagonalization (JJD). Given matrices {M JJD finds an orthogonal matrix U such that U M,U as close to diagonal as possible for all simultaneously. We use the Jacobi rotation-based algorithm that iteratively applies Givens rotations to minimize the off-diagonal Frobenius norm E, ||U M,U diag(U McU) Convergence is declared when the relative change this objective falls below 1
Algorithm 2 Simultancous Matrix Diagonalization (SMD) for Tensor Decomposition
Require:
Noisy symmetric tensor / e R number of random probes Lo 2 2. Ensure: Factor estimates (å, K and weights ( ش such that 7 >X P; 6 Stage 1: random projections - simultaneous diagonalization Lo i.idd. from the unit sphere sK-1
1: Sample [w
2: Form projected matrices (0) + (T,,, w where T(I, ΣKK W. T(anr).
3: Compute an approximate joint diagonalizer via JJD: U10) JJD(M(()).
4: Set (0)) + (U1
Stage 2: inverse-guided projections refinement
5: Build M()) <- (T(I, T, )
k where (0) is the i-th column of ((0).
6: Refine with JJD: (1)) (J).
7: Let à; be the i-th column of (1)) and normalize to ||áu|l2 = 1.
8: for to K do
9: p. + (T, âg @ d; @ a.).
end for
11: return (á.).K (
15

---

Published as a conference paper at ICLR 2026
B.4 TENSOR INITIALIZATION ALGORITHM
With the MLR reformulation and the SMD tensor method, the full algorithm of tensor initialization for MoLDS can be provided now.
Algorithm 3 Tensor Initialization for MoLDS
Require: Trajectories ((u,,0TT-11, (h.0o,-)- truncation L, LDS order n. #components K
Ensure: Estimates of mixture weights PI:k and LDS params {(A,, BRs CR2 Dx)
(A) MoLDS MLR design & sub-sampling
1: for i 1:N do
2: for (T. Ti} do
3:
ūi,t f [u 1...-...... 1 L+11' E 1R
4: Append sample (v,, B5) with vj t us,e/ous Úz Vi,
5: end for
6: end for
7: Partition samples: {1....,MM M) =
where M == E:EiTTi
(B) Moment construction
8: M2 < 2| 2 je Bu, - La) where d = Lm
9: M3 < BIN Σ IjeN, Dj (u E(v)))
(C) Symmetric whiteming and tensor formation
10: Compute rank-K SVD: M2 UEUT, set W U E ET.
11: Form whitened tensor: 7 <- M3 (W, W, W) E RKxXX
(D) Tensor decomposition & recovery
Apply SMD T to recover {åx. Pu k= (Appendix Alg
13: Unwhiten: 3x +- W-T and recover Markov params 9t f- Bu/On.
(E) State-space realization
14: for k = 1:K do
15: Build Hankel matrix from (1
( 3i=1 and apply Ho-Kalman algorithm Recover state-space parameters (A,, Be, Cz, Dz.)
16:
17: end for
18: return ( and ((Ar, BR, Cr DWIlK
Under standard identifiability and excitation conditions, RTPM-based tensor decomposition provably recovers the parameters } with finite-sample guarantees (Rui & Dahleh] 2025). Simultancous matrix diagonalization (SMD) enjoys analogous guarantees the orthogonal-tensor setting and, in practice, tends to be more numerically stable and noise-robust (Kuleshov et al., (015). These advantages lead to improved parameter estimates during the tensor initialization stage (Figure 2 main paper), which in turn can provide stronger starting point for the subsequent EM refinement.
C EM FOR MOLDS: COMPLETE TECHNICAL DETAILS
C.1 OVERVIEW AND STRUCTURE
This appendix provides complete technical details for our EM formulation for MoLDS. The EM algorithm for MoLDS alternates between two phases: (i) an E-step that computes trajectory- responsibilities via Kalman filter likelihoods and extracts sufficient statistics via Kalman smoothing, and (ii) an M-step that updates all parameters using closed-for -oom maximum likelihood estimates from sipposibilltyywwighteee statistics. The algorithm is detailed in Algorithm 4
C.2 COMPLETE EM ALGORITHM
The EM procedure operates on trajectories () with current parameter esti- mates ((t) (pkk Åg, Be, Ćt, De, Qk, K k= iteration 1.
16

---

Published as a conference paper at ICLR 2026
E-step Computations. For each trajectory i and component k, we compute:
t, log 11,0::111 1 -1.0:2;11;0%)) )²)
(8)
01,k = log (1) + bi
(9)
exp(ark)
7i.k E exp(air))
Yi.k =
(10)
The likelihood & is computed via the Kalman filter, while responsibilities 7 the log- -og-sum-exp trick for numerical stability. Subsequently, the Kalman smoother computes per-trajectory sufficient statistics Si,k, which are aggregated as:
S
7t.k Si,k-
(11)
M-step Updates. Mixture weights and LDS parameters are updated via:
Pk
(1++)
1
Σ
(12)
N
((++)) MLE-LDD((E)).
(13)
where the closed-form LDS parameter updates are derived in C.3
Convergence. The algorithm monitors the observed-data log-likelihood:
N
c(++1)
log
Pk
(++1) p p(
us,0:T, 0g t
(14)
and terminates when the relative improvement falls below the threshold e.
Algorithm EM Refinement for MoLDS
Require: Trajectories {(u,00T,-11 -)) initial parameters ((0); max iterations maxx;
tolerance E
Ensure: Refined mixture weights and LDS parameters
1: for iter = 1 to EM_max do
E-step: Compute responsibilities and sufficient statistics
2: for i 1: do
3:
for = 1: K do Compute 4,, via Kalman filter on (u4,0:1 -1,00:;-11122 01,0.,.1) using @(itor-1)
4:
Set Ori,k + log ffiter-1) + fi
5:
6:
Run Kalman smoother compute sufficient statistics Si,k end for
7:
8:
Compute Ise; k exp(as.r) and responsibilities Yr, + exp(ar, Ise;) end for
9:
Aggregate responsibility statistics: S +- Σ 7i.k Si,k for each M-step: Update all parameters
Update mixture weights: pk ((itor)
1 Dim Ne.k for all k
12:
for & ==11 : K do
Update LDS parameters (Åe, Br, Oc, Dr, Quy Rie) from Sa (see C.3
14:
end for
Convergence check
15:
Compute observed log-likelline ligllikllihood: ((terr) + Σ Ise;
16: Stop if (Cite))) - Cliter C(iter-1) /) <e
end for
18: Return {pr, (Áe, Br. Cx, Dr, Rk)k
17

---

Published as a conference paper at ICLR 2026
C.3 CLOSED-FORM M-STEP UPDATES
The M-step updates follow standard LDS maximum likelihood estimation (Ghahramani & Hinton 1996). For simplicity, we present the case Di 0 (no direct feedthrough) and drop component and iteration indices (k, t) for readability:
We solve the normal equations with optional ridge regularization À each component.
to obtain system matrices for
(
S S Saaa
ST Sunu
-1
[4 BJ = 15az- Sza-]
+ x1
(15)
C Sym (Szo + X1)
(16)
And the noise covariances are
T (Snu CS yz
- CyaCT + CSZzaC)).
(17)
1 Q N (Sux,eurr AST gx- -BST (AS 1x- BS zw-)
(18)
AS A AS AS B BS A BS BS
(19)
Here, all S.. are the sufficient statistics.
C.4
COMPUTATIONAL COMPLEXITY
Each EM iteration requires ((NKn operations:
Kalman filtering: Onn²7) per trajectory-component pair, total O(NKnTT)
Kalman smoothing: TT) per trajectory-component pair, total O(NKn
M-step updates: O(Kn²) for matrix inversions
Tensor initialization significantly reduces the iteration count compared to random initialization, which substantially improves computational efficiency. The tensor initialization phase requires (d8) operations where d Lm is the MLR dimension, making it negligible in the whole Tensor-MM pipeline.
D INITIALIZING (Q, AFTER TENSOR INITIALIZATION
The tensor stage yields (Å_, B., C but not the noise covariances (Q., R.). We initialize (Qz, R2) from data for each component as follows.
Step (labels/weights: optional). Assign labels by selecting, for each trajectory i, the component with the smallest one-step prediction MSE under (A,, Bz, C²),, and w..s for that component and Wi,r for all # z.
Step (state buck-projection). For each z and trajectory i, decode provisional latents with a ridge pseudo-inverse
(CI C, + A1) CT (Vi. D_u,,)) where we recommend 1 max (diag (CT C,)) for numerical stability.
L
Step 3 (residual covariances). Define
(x) = 1,1+1 A, B B.us
(z) = Yi C, D_ug.
Nalt
With T, = 1, set
Σn Winz
2t=0 1/tt Titt
Ti -2-2 (=) (z)T
R1O0 E. Wi. DiWi, DT (=) (=]T
6 Bit
Σ; Wi,; T.r
E, W;, Tj
18

---

Published as a conference paper at ICLR 2026
Step (positive semidefinite projection). Symmetrize and project to the positive semidefinite cone: e + T (Q + ))
R < IT IT ( + )),
where o M) projects matrix M the nearest positive semidefinite matrix via eigendecomposition:: if M UAUT, then (M) U max(A, 0YUT.
This yields (Q.), ) for each component 2, providing a stable initialization for the first E-step; subsequent EM iterations refine (Q., R2) from responsibility sufficient statistics.
E
TENSOR-EM MOLDS APPLICATION ON NEURAL DATASET DETAILS
E.I
AREA2 DATASET
Dataset introduction. We evaluate our methods on a neural population dataset recorded from the motor cortex (Area2) of rhesus macaque. The task is center-out reaching paradigm where the subject applies force to manipulandum in response to cues. Neural activity was recorded with a Utah array, and spike counts extracted from thresholded multiunit activity. The dataset provides simultaneously measured behavioral covariates (applied force, hand kinematics) alongside the neural recordings.
Data preparation for MoLDS. Following the Neural Latents Benchmark preprocessing, we first apply principal component analysis (PCA) to reduce the neural activity to p latent dimensions (here we use p 6, 20). The behavioral force signals (6D) are used as exogenous inputs. For MoL.S training, we constructed input-output trials of the form (w,, (u4, where u, IR'6 corresponds to force inputs and E I are PCA-reduced neural outputs. Each reaching direction corresponds to a separate trial. We split the dataset into non-overlapping train, val, and test sets, ensuring balanced coverage across each direction.
checked the variance explained by PCA and found that the first six PCs capture over 90% of the total variance for each direction, while the first 20 PCs explain nearly all of it. This indicates that PCA preserves the essential structure the neural activity. The following figure illustrates this.
Noural sclltty Conditor!
Cumslaivve explaited
30
03
209
Hamder
Figure 7: Area2 neural activity PCA spectrum.
Tensor-EM MoLDS train/val/test setup. For model initialization, we computed input-output moment tensors from the training data and applied SMD to obtain globally consistent estimates of system parameters and mixture weights. These estimates were used to initialize a full Kalman filter-smoother EM refinement, which updates all LDS parameters in closed form. We trained MoLDS models with different numbers of mixture components
E 3, 4, 5, latent dimension n E 3, 4, 5, and lag parameters L 16, 24. Validation data is used for model selection by comparing negative log-likelihood (NLL), root mean squared error (RMSE), and Bayesian Information Criterion (BIC).
NLL: sum of Gaussian-innovatioo log- ligllikeihooos.
RMSE: root mean squared error of one-step predictions ÿr vs. yt-
BIC: BIC = -2 NLL + Po log Nabs. where Py is the total parameter count and Noss the
number of observed outputs.
19

---

Published as a conference paper at ICLR 2026
andedoeree LD * dM
Figure 8: Markov parameters representation of LDS models from different fitting methods.
The best configuration was tested on the held-out test set.
Per-trial responsibilities.
We compute responsibilities for each trial under component k
K
Tik exp
1,
where ohe is the one-step Kalman log- likeliooddo of trial under component k. Global usage is
summarized by usage; N Σ, Tik.
Analyses.
1. Validation metrics. We plot BIC/NLL/RMSE across K and initializations (Tensor vs. Random) to assess stability and choose capacity.
2. Representative predictions. We visualize Ú: vs. J/ for representative test trials (all p outputs).
3. Direction organization & component usage. In a polar wheel with one wedge per discrete direction, we color each wedge by the dominant component (highest mean Tak among trials in that wedge).
4. Markoy response comparison. For each component, we calculate g (T( = C& E IRxm O..... L (omit Dj for clarity). Plot (1)x(T)] for all components single axis to compare gain/decay.
5. Global geometry of dynamics. Vectorize all {gu (7)),, run PCA across that set, and scatter the first two PCs; use color to encode lag 7 and marker/edge to encode component. Components typically form distinct, smoothly evolving trajectories.
Supervised per-direction LDS fitting baseline. Because Area2 discrete directions, also fit a single LDS per direction using the same observation space. We compute each supervised model's impulse response, cluster the per-direction vectors, and align clusters to MoLDS components via a Hungarian assignment on centroid distances. A polar plot (inner = supervised clusters, outer dynamies.
MoLDS predictions) confirms that unsupervised MoLDS recovers
Single-trial LDS fitting baseline. In addition, we also fit an LDS for each single trial. We compute the average LDS impulse response in each direction and compare the obtained clusters as compared to the per-direction LDS. Following the same procedure as in per-dir LDS fitting, we compare the clusters of single-trial LDS fitting with MoLDS components as well. We have examined the trial- to-trial variability of the Area2 dataset by comparing their full Markov parameter vectors between single-trial, perd-dir, and single global LDS fits. As shown in Figure g by evaluating the first three PCs of the Markov vectors, single-trial LDSS within each direction form compact clusters, and the
20

---

Published as a conference paper at ICLR 2026
aa)
(b)
(c)
$1.05 trtal grouping by state Edominant st)
SLDS: mean time share by state
407
20"
997
232*
3207
186
0.4
180
3
208
278"
02
312"
40
Mtat
Figure 9: SLDS result: Inferred switch states for all test trials with mild (a) and strong stickiness (b). (a) Group trials and mean time by states with the mildly sticky SLDS fit.
pooled per-direction LDS lies near the cluster center for most directions. In contrast, these single-trial
LDS are very far from global LDS fitted on all trials. This indicates that within-direction variability mild, and the per-direction LDS provides a stable summary of the underlying dynamics rather than oversimplification. We also performed clustering directly independently fitted single-trial LDSs. These LDSs naturally cluster by movement direction, even without pooling, demonstrating that direction-specifie dynamical signatures already present the single-trial level. Moreover, by comparing with the per-dir LDS and MoLDS fit, these single-trial LDS fit- btaed clusters are very similar, with only two directions being switched between clusters (see Figure [2]. Together, these analyses show that while the per-direction LDS is not perfect "gold standard," it is a reasonable and stable reference, and our comparisons to MoLDS capture meaningful structure the data.
SLDS comparison. We fit switching LDS (SLDS) models with S E 1, 2, discrete states to the training (Gaussian emissions, 6-D force inputs, mild strong sticky transitions). Validation NLL/BIC selected S 2 SLDS, We then decoded the test set with Viterbi and summarized the inferred state sequences (Figure D. Panel (a) shows that within-trial switches are infrequent and tend to occur in similar temporal regions across trials. In panel (b, left), grouping test trials by their dominant state reveals that the same state dominates most directions, with only a mild departure around dir-90. The mean time-share view (b, right) tells the same story: state composition looks broadly similar across directions. Consistent with this, a "collapsed" SLDS that uses only each trial's dominant state achieves nearly the same predictive metrics as the fully segmented model, indicating that within-trit -riil switching adds little in this dataset. These results support our interpretation that the primary structure is between trials (by direction), not within trials.
More results.
Impulse response comparison. We compare the direction-specifie impulse response (Markov pa- rameter) curves of the assigned MoLDS component and the corresponding single LDS fit for that direction, which shows high similarity (see Figure 10 where titles report cosine similarity).
Results MoLDS fitting using 20 PCs. We also applied the MoLDS method on the Area2 dataset with 20 PCs, and found consistent results for the 6-PC fitting model shown in Figure 11
21

---

Published as a conference paper at ICLR 2026
Por-direction impulse curves vs nearest MoLDS component
componert
Imgulen nomm
5.03
g" come 77
45"
002
0.00
90"
C05.0.75
0.06
135
co500 50
0.03
0.02
0.01
0.01
Impuiee rorm
6:06
0.02
6.04
0.01
0.21
10
20
10
lag
Figure 10: Impulse response comparison between MoL.DS components and per-dir L.S fit.
(a) Validation metrics
(b) Test 1-step predictions
dim
dim
clm
The Naee Maree
4.0
14
3.5
3.0
2.5
Thai Doest ke=21
dmm
cim
dim
0.5
20
+
0.4
0.8
15.
2
0.3
Randomed 1orcar-81
Rosser
Random
Random
Recson Rarser
Method
(c1) Single-trial LDS| MoLDS
90
(d) MoLDS component usage
135
45
Unage fractinn
1.00
0.57
180
9,50
0.25
0.23
225
315
225
315
Comp
Comp
Comp
270*
330*
Figure 11: MoLDS application Area2 with 20 PCs: results consistent with those presented in the main text.
22

---

Published as a conference paper at ICLR 2026
det
dtte
dn
dms 4
Figure 12: Multi-step forward predictions of MoLDS on one trial of Area2 data: 5-step vs fully free run after 100ms warm-up.
Dominant second component diagnostics
Distribution of dominant component responsibilities
Distribution of second component responsibiliitie
10.0
10.0
7.5
0.0
0.
responsbilty
Per-direction
responsibiliity
component per
1.0
companaann (oden most)
30
0.8
25
200
0.8
0.0
03
(0-7)
Figure 13: Area2 MoLDS fit dominant and second component analysis.
Multi-step prediction In addition the one-step-ahead predictions presented in the main paper, we also inspected multi-step and free-run predictions. A representative example is provided in Figure The model follows the true neural trajectory for several steps before diverging as runs forward without correction, which is the characteristic and expected pattern for linear dynamical systems. This behavior consistent with LDS-based modeling and does not affect our use of MoLDS for identifying shared local dynamical structure across trials.
12
Responsibility distribution. We analyzed the full responsibility distribution of MoLDS components on held-out test trials to assess whether non-dominant components carry substantial probability. Across all trials, the dominant responsibility is very large 0.976, while the second-largest responsibil- ity is very small 0.023. This pattern is consistent across directions shown in Figure [3
E.2 PMD DATASET
The dorsal premotor cortex (PMd) dataset contains single-trial reaches, where time-aligned kinematics positioovvelocitylcceleaaioo and multi-unit spikes, plus reach direction angle E (-x, n] are stored. For MoL.DS, the input and observation are taken as
Dataset introduction.
U+ [Uz, Byy any ays speed] RR.
Y4 E R6
23

---

Published as a conference paper at ICLR 2026
Inpust wise cerrp
PCA vectoriend iepoine Mede
002
cemp
002
cemp
1400
PCD
camp
000
PC3 (42-00
Figure 14: Input specificity and geometry of Markov responses: Left: input-wise energies versus lags are shown per component, Right: PCA of vectorized impulse blocks (marker shape component, color lags).
where speed v0 ug and Y, denotes PCA-reduced, 2-scored firing rates. Here p 16 PCs, which explain 90% of neural activity. As in Area2 (App. E.1 angles used only for analysis/visuaalzzaton; MoLDS is trained without angle labels.
Tensor-EM MoLDS train/val/test setup. Training and selection mirror Area2: tensor or random initialization followed by EM refinement; Model selection was performed on the validation split, For PMd, we sweep (2, 3, 4), n (4, 5, 6), and Markov horizon L (12, 16, 24). The selected model in our experiments is K K44, 5= 1-12.
Analysis idea. We reuse the analysis logic from Area2. As the direction angles are not discrete as in Area2, split them into 12 bins, and each bin is colored by the dominant component (largest mean rik the bin). Component wise preferred directions are shown by responsibiliyy-weighhed arrows
UK Σ Tik [cos 6,, sin 0,],
plotted as an inset (angle Luk, length [|v,])). For each component, we examine the Markov parameter blocks
(T( CRARBE E RPX" = O...., 1-
and plot (a) the compact magnitude curves [T]]( and (b) per-input energies Ex,(T) = (19H()). ||2 for j E {Uz, Byy az ay speed]. For global view, we also vectorize (g)(T) run PCA, and scatter the first two PCs.
Figure 14 (a) shows per-input energies versus lag, revealing input selectivity and decay rates. Fig- ure (b)b) visualizes the geometry of vectorized impulse blocks, where marker shape represents component identity and color denotes the lag. These demonstrate that MoLDS components exhibit. separated clusters corresponding to component-specife dynamical subspaces.
24


Reviews:
[{"review_id": "mi7ohFg1xs", "reviewer": "Reviewer_jFVZ", "summary": "This paper develops a new approach towards fitting Mixture of Linear Dynamical Systems models, using tensor-based methods followed by EM refinement. The paper demonstrates the superiority of their model-fitting approach on simulated data. The paper then fits the model to two neural activity datasets.", "questions": "1. It took me a very long time into the paper to really understand how MoLDS was different than an SLDS. It would help readers a lot if you more explicitly compared them early on in the methodology, since SLDS is much more common in comp neuro (which will be the primary readers of this paper).\n2. I'd clarify why the RTPM method is being used for a comparison - what past paper was this used in?\n3. It seems very odd to be doing PCA on Area 2 data prior to fitting MoLDS, when generally an advantage of LDS models is that they are directly learning the underlying latent states from high-D data.\n4. Relatedly, does the PMd experiment also use PCs or all neurons?\n5. Line 72-79 - if possible, some citations for pros/cons of approaches could be useful (although I realize sometimes these insights are from personal experience)", "limitations": "", "rating": 2, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "-The authors demonstrate a model that (to my knowledge) hasn't been used in neuroscience data before\n-The authors show convincing results of the benefit of their new tensor-EM fitting approach in simulated data (and quantitatively versus standard EM approaches on one real neural dataset)", "weaknesses": "The importance of MoLDS for modeling neural data (which is the focus of this paper) never got across clearly to me. It feels like the authors take that importance as a given, but this is a method that hasn't been used (to my knowledge) in neural data, so really clearly motivating the need is essential - at the moment, it's primarily a sentence on line 60-61. When are specific examples when a neuroscience researcher would want to use MoLDS when existing approaches would fail? The real data examples don't make this clear, since it seems like you could just fit LDS models to different reach directions in those scenarios (they don't show any benefit of an unsupervised approach), and it's especially unclear what to take away from the results on the PMd dataset. \n\nAdditionally, the comparisons to alternative approaches for the real datasets were lacking. There was not a comparison to the tensor approach to fitting, so it's currently not clear whether this paper's fitting approach is advantageous relative to pure tensor methods. Additionally, while the paper shows a numerical advantage versus a pure EM approach, there's no advantage shown in terms of the model analysis (e.g. within Fig 4c). So in general, the advantages of this paper's model-fitting approach aren't currently shown strongly.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "MplFEZXzER", "reviewer": "Reviewer_Lrdy", "summary": "This paper presents a new fitting procedure for modeling heterogeneous time-series data based on Mixtures of Linear Dynamical Systems (MoLDS) by introducing a robust Tensor-EM algorithm. The authors propose combining tensor-based moment estimation with EM refinement and show that this method achieves more accurate and efficient results than alternative methods based on single LDSs (Linear Dynamical Systems) or switching LDSs. They also extend the application to neural data and exemplify the method on complex recordings of a reaching task and neural recordings from premotor cortex to demonstrate the applicability of their MoLDS framework to neuroscience.", "questions": "1) Could you clarify how you performed the SLDS comparison? Specifically, did you take the most active operator from SLDS trained on all trials together? How did you define the number of operators? What happens if you include stronger penalization on the number or frequency of switches per trial? I would assume that with strong penalization you may be able to capture the individual dominant operator as in MoLDS.\n\n2) Regarding my first concern about non-stationarity within a trial, I am curious what happens if you split each trial into multiple periods (e.g., stimulus presentation vs. reaching period) and run your method treating these as separate trials. Do you get the same operator being active in both periods, or do they switch?\n\n3) Have you examined multi-step reconstruction instead of single-step?\n\n4) It seems you focused mainly on the dominant component from MoLDS. I am curious about the distribution of probabilities for different LDSs within a single trial. Is the probability of the second component, for example, very close to that of the first? If so, what does the second component look like? Is it also shared between trials of similar direction? (I.e., for a fixed and large K, unlike Figure 4b where you appear to modulate K).\n\n5) Could you explain the rationale for choosing a test-val-train split over cross-validation?\n\n6) Could you clarify the difference between Sections 4.1 and 4.2 in terms of which version of your method was used? For instance, in Section 4.1, did you use only SMD with randomly initialized EM?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1) The authors enhanced the MoLDS model to achieve improved numerical stability and robustness to noise via Simultaneous Matrix Diagonalization (SMD) compared to prior tensor methods.\n2) The authors extend the tensor initialization paradigm by integrating it with a Kalman filter-smoother EM procedure.\n3) The method, which uses LDSs, maintains representation interpretability.\n4) The authors demonstrated their results on two synthetic experiments and two real-world neural datasets.", "weaknesses": "1) My biggest concern relates to the extension of MoLDS to neural data. Neural activity exhibits complex temporal evolution that changes over time, potentially even within a single trial. While I understand this approach enables interpretability by using a single LDS per trial, it cannot capture the temporally changing nature of brain dynamics, which is a major question in neuroscience. Hence, I wonder whether neuroscience is the right application to demonstrate the method. Alternatively, how could the method enable understanding of varying brain interactions over time? For example, could your fitting procedure be extended to more complex LDS models that change over time (e.g., SLDS or dLDS)?\n\n2) It seems you treat the \"per-direction single-LDS\" as a gold standard. I understand that in real data we lack access to the true operators, so estimation is necessary. However, I am concerned that per-direction single-LDS oversimplifies the dynamics. Given the complexity of neural data, including trial-to-trial variability due to unobserved factors even under the same hand direction, this baseline may yield biased comparisons.\n\nHave the authors examined trial-to-trial variability within a direction and its effect on your comparison? Specifically, what happens if you train both your method and the per-direction single-LDS on random subsets of trials within a direction, repeating this process with different subsets? Does performance remain stable?\n\nWith respect to that, I also think a missing comparison is to train a single LDS per trial, then perform clustering on these independently trained LDSs to verify whether they capture trial direction as you implicitly assume. After validating that LDSs cluster by direction (without forcing a single LDS on all trials from a direction), you could then compare this to your method\n\n3) I think the paper could be written more clearly, as it is currently difficult to follow. The utility and importance of certain components are unclear. For example, more background information on SMD and its contribution to stability would be helpful in the related work section. Additionally, please explain more explicitly how SMD was applied in your method. I would also suggest dedicating space to explaining and defining the moments and W in Eq. 3. Finally, in the LDS paragraph (lines 54-59), there should be discussion of other LDS models that support multiple LDSs with mixtures, like dLDS [1,2].\n\n4) The abstract is very long and difficult to follow. I suggest condensing everything from `We validate our framework on...' to the end into a single sentence, which would also save space in the main text.\n\n5) I am not sure why the method illustration appears only in Figure 3. This illustration is helpful for understanding the paper and should appear earlier\n\n5) Some smaller comments on writing and presentation issues that should be fixed: 1) Please pay attention to \\cite vs. \\citep, it seems all references appear as part of the text (no brackets), which hinders readability. 2) Hyphens are used instead of em dashes in several places. 3) In Fig. 4b, the subplots are ordered oddly (dim 4-0-1-5...). \n- Some typos and formatting issues: extra bracket in line 466, double (c)(c) in line 240, missing space before reference in line 060.\n\n[1] Mudrik, N., Chen, Y., Yezerets, E., Rozell, C. J., & Charles, A. S. (2024). Decomposed linear dynamical systems (dlds) for learning the latent components of neural dynamics. Journal of Machine Learning Research (JMLR)\n\n[2] Chen, Y., Mudrik, N., Johnsen, K. A., Alagapan, S., Charles, A. S., & Rozell, C. (2024). Probabilistic decomposed linear dynamical systems for robust discovery of latent neural dynamics. Advances in Neural Information Processing Systems (NeurIPS)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "7tyQhMzu00", "reviewer": "Reviewer_tLS4", "summary": "The paper proposes a tensor-based moment method for learning mixtures of linear dynamical systems (MoLDS). The authors use Simultaneous Matrix Diagonalization on the moment tensors constructed from the input-output data to recover globally consistent estimates of mixture weights and system parameters, which are then refined with the Kalman smoother. The paper applies the proposed method to both synthetic and real-world neural datasets.", "questions": "I am curious about how the proposed method will perform as we scale up the size of the dataset in terms of the observations and latent dimension. In general, I would be interested in stress testing the proposed method by varying different experimental axes.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "* Clarity: The paper is extremely well written, with clear motivation and definition of the problem, related work, and where the proposed method stands (i.e., Contributions). It also has a clear explanation of the method with appropriate equations and algorithm pseudocode, and the Figures are well presented with descriptive captions.\n\n* Evaluation: The method is tested extensively on both simulated and real-world neural datasets. The results are well presented, showing the advantage of the proposed method over the compared methods.\n\n* Significance: The paper is well motivated, setting up the existing problem of the difficulty of learning MoLDS despite its usefulness. The paper tackles this problem and demonstrates its effectiveness in analyzing neural data.", "weaknesses": "As the authors noted, one of the limitations of the paper is that the proposed method is only evaluated on datasets that have linear dynamics. A more detailed explanation of what could happen when there’s a model mismatch, in addition to how the method could be extended to mitigate the mismatch, would be valuable.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "J1l4gsyikV", "reviewer": "Reviewer_vVHY", "summary": "The paper proposes a two-stage pipeline for learning mixtures of linear dynamical systems (LDSs). Stage one follows the tensor-decomposition framework of Bakshi et al. to obtain consistent initial estimates under standard conditions. Stage two refines all LDS parameters via an EM procedure that uses responsibility-weighted sufficient statistics (updating dynamics, emissions, and initial states). A principled initialization for the noise covariances Q and \nR based on residual covariances is introduced. The method is applied to neural recordings from primate somatosensory cortex during reaching, illustrating the approach on a scientifically meaningful dataset.", "questions": "What further applications to the authors see for mixtures of LDS’s?  What further improvements can be made to their statistical pipeline?", "limitations": "", "rating": 6, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The paper is exceptionally well written and easy to follow.\n\nThe two-stage design is natural: tensor methods provide strong initializers with theoretical guarantees, while EM addresses the well-known limitations of pure tensor approaches by locally refining parameters.\n\nAssumptions are clearly stated and appropriate (persistent excitation; controllability/observability; component separation), with careful discussion of when they hold.\n\nA notable contribution is the explicit treatment of Q and R, which are often sidelined in prior tensor-based work; the residual-covariance initialization is principled and practical.\n\nThe overall study exemplifies well-designed statistical methodology: theory-grounded, reproducible, and attentive to the mechanics of LDS estimation (including filtering/smoothing considerations).\n\nThe application to neural data is compelling and broadens the paper’s relevance beyond synthetic studies.", "weaknesses": "I can’t identify any obvious weaknesses, other than that it’s hard for me to evaluate data collected on primates beyond evaluating the statistical methodology, which as I’ve stated before is rock solid.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)