Learning linear state-space models with sparse system matrices
Authors: Yasen Wang
Kaiqi Fang
Guijun Ma
Junlin Li
Mengyu Sun
Zhilan Huang
Gang Lu
Venue: ICLR 2026 Poster (2026)
Keywords: linear state-space models
expectation-maximization algorithm
system identification
state estimation
Primary area: learning on time series and dynamical systems

Abstract: Due to tractable analysis and control,  linear state-space  models (LSSMs) provide a fundamental mathematical tool for time-series data modeling in various disciplines. In particular, many LSSMs have sparse system matrices because interactions among variables are limited or only a few significant relationships exist.  However, current learning algorithms for LSSMs lack the ability to learn system matrices with the sparsity constraint due to the similarity transformation. To address this issue, we impose sparsity-promoting priors on system matrices to balance modeling error and  model complexity. By taking hidden states of LSSMs as latent variables, we then explore the expectation-maximization (EM) algorithm to derive a maximum a posteriori (MAP) estimate of both hidden states and system matrices from noisy observations. Based on the Global Convergence Theorem, we further demonstrate that the proposed learning algorithm yields a sequence  converging to a local maximum or saddle point of the joint posterior distribution. Finally, experimental results on simulation and real-world problems illustrate that the proposed algorithm can preserve the inherent topological structure among variables and significantly improve prediction accuracy over classical learning algorithms.

Published as a conference paper at ICLR 2026
LEARNING LINEAR STATE-SPACE MODELS SPARSE SYSTEM MATRICES
WITH
Yasen 2 Kaiqi Fang", Guijun Ma 3,4*
Junlin Li", Mengyu Sun Zhilan Huang', Gang Lu
CCina Telecom Research Institute, Guangzhou, China
Schooo of Future Technology, South China University of Technology, Guangzhou, China
3 Schooo of Artificial Intelligence and Automation,
Huazhong University of Science and Technology, Wuhan, China
"State Key Laboratory of Intelligent Manufacturing Equipment and Technology,
Huazhong University of Science and Technology, Wuhan, China
"Scoool of Mathematics and Statistics, Fuyang Normal University, Fuyang, China
"China Telecom Research Institute, Beijing, China
ABSTRACT
Due to tractable analysis and control, linear state-space models (LSSMs) provide fundamental mathematical tool for time-series data modeling in various disci- plines. In particular, many LSSMs have sparse system matrices because inter- actions among variables are limited or only a few significant relationships exist. However, current learning algorithms for LSSMs lack the ability to learn sys- tem matrices with the sparsity constraint due to the similarity transformation. To address this issue, impose ance modeling error and model complexity. By taking hidden states of LSSMs as latent variables, we then explore the pepecttaioonmaximmzzttttt derive maximum a posteriori (MAP) estimate of both hidden states and system matrices from noisy observations. Based on the Global Convergence Theorem, we further demonstrate that the proposed learning algorithm yields a sequence converging to local maximum or saddle point of the joint posterior distribution. Finally, experimental results on simulation and real-world problems illustrate that the proposed algorithm can preserve the inherent topological structure among vari- ables and significantly improve prediction accuracy over classical learning algo- rithms.
priors on system matrices to bal-
(EM) algorithm
INTRODUCTION
Linear state-space models (LSSMs) are fundamental mathematical tools for analyzing time-series data with applications in robotics (Mamakoukas et al.] 2019, 2020), systems biology (Jin et al.]
2020b Pillonetto & Ljung 2023), and natural language processing (Smith al., 1999 |Belanger & Kakade 2015). Generally, LSSMs deseribe time-series data {(4., y:) through the following stochastic difference equation:
Ax + Bu, t
(1)
y: Cr, Du: + wi.
(2)
where Un IR is the input signal, yı Rd is the noisy observation, 2; E R" is the hidden state, A E Ruxx B E RTxp, C E mp and E 1p=xp are the unknown system matrices, and E, (O, R) and u. N0, Q) the diagonal process and measurement noise, respectively. In addition, LSSMs are also widely used to approximate complex non-linear systems in industrial processes given their relative simplicity (Yuan al.) 2017 Lusch
2018). Due to a complete rigorous theory available on LSSMs, learning them from noisy observations can enable us to make
tractable analysis and control of systems (Chen & Poor 2022 Bakshi et al. 2023)
"Correspondence to: Yasen (angg Gujjun Ma (mgj @ hust.edu.cn).. Source codes are available on GitHub
Saa.s//paaannnnaaaa-a--ee- Sparse- 1.SMM

---

Published as a conference paper at ICLR 2026
In this paper, we focus learning LSSMs with sparse system matrices for two important reasons. First, the learned LSSMs should include the minimally required parameters to explain time-series data following the Occam's 's razor principle, which favors explanations constructed with the smallest possible set of elements. Additionally, many real-world systems indeed have a sparse topology, as each state measurement variable only depends on few other state variables and inputs (Efroni iet al. 2022). For example, gene only regulates the expression of a limited number of other genes in gene regulatory networks (He et al. (2024b). In industry, communication systems usually have sparse topology to reduce energy consumption (Jin et al. 2020a b). However, available learning algorithms lack the ability to learn LSSMs with the sparsity constraint on system matrices due the similarity transformation.
To learn LSSMs with sparse system matrices, we impose sparsity-promotte priors on them to bal- ance model complexity and modeling error. Following the Bayes' rule, we can combine the marginal likelihood and prior functions to derive the joint posterior distribution of all unknown variables. However, directly maximizing such posterior distribution to estimate system matrices is intractable because the hidden states of LSSMs are unknown. To address this issue, we explore the expectation- maximization (EM) algorithm to give an alternate maximum a posteriori (MAP) estimate of hidden states and system matrices by taking hidden states as latent variables. In the expectation step, we use the Rauch-Tung-Stre Strrebbee (RTS) smoother give a closed-form update rule for the hidden states. In the maximization step, we leverage the block coordinate descent method to analytically update the system matrices in turn. By alternately performing the expectation and maximization steps un- til convergence, the proposed algorithm can determine the sparse system matrices of LSSMs from noisy observations. summary, the contributions of this paper are threefold:
Leveraging pparsit-prroooting techniques, we propose an algorithm to learn LSSMs with sparse system matrices from noisy observations. Following the Global Convergence Theo- rem (Luenberger et al.) 1984). we also demonstrate that the proposed algorithm is guaran- teed to converge to a local maximum or saddle point of the posterior distribution composed of the marginal likelihood and prior functions.
Because available learning algorithms only learn LSSMs up to a similarity transformation, the learned system matrices usually differ from the true ones in both numerical values and topological structure. However, the proposed algorithm learns system matrices by bal- ancing model complexity and modeling error. As a result, the learned system matrices can preserve the inherent topological structure among variables, which is a significant improve- ment over classical learning algorithms.
Experimental results on simulation and real-world datasets demonstrate that the proposed algorithm outperforms classical ones on learning L.SMMs with sparse system matrices. In addition, the learned system matrices of the proposed algorithm are more valuable for ex- ploring the interaction laws of systems.
RELATED WORK
Least-squares minimization. Basically, least-squares minimization (LSM) learns unknown models by minimizing the sum of the squares of residuals (Faradonbeh et af.] 2018 2020) Modi et al. 2024]. By taking the one-step prediction error as the objective function in LSM, prediction error minimization (PEM) is proposed to learn LSSMs via gradient-based optimization methods (Ljung) 2002 Katayama et al. 2005). Given a symmetric transition matrix, Hazan et al. (2017) design an efficient method for the online prediction of LSSMs by formulating system identification as an online PEM problem. Recently, combining the Kl-Klamn (HK) algorithm with LSM, Oymak & Ozay 2019) propose a method to learn system matrices of I LSSMs with sample complexity analysis. However, it is well known that LSM sensitive to noise and cannot characterize the sparsity of system matrices (Tibshirani 1996 Martens 2010).
Subspace state-space system identification. Subspace state-space system identification (4SID) al- gorithms project data Hankel matrices onto certain subspaces to estimate the extended observability matrix and hidden states using linear algebra tools (Larimore) 1990; Verahegen & Dewilde 1992. Van Overschee & De Moor 1994; He et al.) 22224)). Subsequently, system matrices can be recovered from either the extended observability matrix or hidden states (Favoreel et al., 2000). Based on prin- cipal component analysis, [Wang& & Qin (2002) present a new 4SID algorithm to learn LSSMs under
2

---

Published as a conference paper at ICLR 2026
the errors-in- varriable situation. By choosing different weighting matrices perform the singular value decomposition. Van Overschee & De Moor (2012) provide a geometric framework to unify almost all classical 4SID methods. Further, [Huang et al. (2016) present the Weight-Least-Squar method to learn stable LSSMs by multiplying the unstable component with a weight matrix. How- ever, it is widely recognized that such algorithms generally cannot obtain accurate system matrices as required (Martens) 2010; Qin 2006).
Maximum likelihood estimation. Because the joint likelihood function of LSSMs involves hidden states, the EM algorithm is employed to give the maximum likelihood estimation (MLE) of system matrices (Shumway & Stoffer 1982, Ghahramani & Hinton| 1996). Leveraging the EM algorithm, the distribution of hidden states can be explicitly derived using the Kalman smoother based on the current estimate of system matrices. then updates system matrices by maximizing the expected log-likelihood with respect to the hidden states. To present robust MLE for LSSMs, Gibson & Ninness (2005) implement the expectation and maximization steps via the LR and Cholesky factori- sation, respectively. To increase the efficiency of EM for learning LSSMs, Martens (2rtens)(010) proposes an approximate second-order statistics (ASOS) scheme to approximate the expectation step. Com- bining EM and Lagrangian relaxation, |Umenberger et al. (2018) use semidefinite programming to optimize the tight bounds on the likelihood to learn LSSMs with model stability constraints. How- ever, such learning algorithms lack the ability to deal with sparse system matrices.
Oparsity-promotie methods. By adding a penalty term on model parameters, methods can balance model complexity and modeling error learn systems from data (Brunton et al. 2016). Leveraging the { regularization term, Tibshirani (1996) proposes method named Lasso estimate parameters in linear models. Further, reweighted ε minimization is proposed to enhance sparsity (Wipf& & Nagarajan] 2007 Candes et al., [0008] Garrigues & Olshausen 2010)) However, solving an & minimization problem is challenging due to dioodifffrentiaaiiiiy at the origin, and these methods also require careful fine-tuning of hyperparameters. To address such issues, sparse Bayesian learning (SBL) imposes sparsity-promoting priors on model parameters to enforce sparsity (Samanta et al., 2022 Chakraborty et al. (023). Subsequently, it maximizes the posterior distribution consisting of the likelihood function and priors to estimate model parameters and hyperparameters (Tipping 2001 Wipf & Rao 2004). Recently, SBL has been applied to learn various systems from data, with system states being measurable yet potentially corrupted by process noise (Pan et al., 2013 Yuan al. 2019: Wang al. 2024). However, leveraging such sparsity- promoting methods to learn LSSMs with sparse system matrices remains an elusive and challenging problem because system states are unavailable and observed data are corrupted by both process and
measurement noise (Course & Nair 2023).
3 METHODOLOGY
Due to the similarity transformation, LSSMs admit many equivalent representations with different levels of sparsity, where the corresponding transformed system matrices are given by $ BB, C -1 and D, with & E T being a nonsingular matrix. However, we focus on learning the LSSMs with sparse system matrices that include minimally required parameters in accordance with the Occam's 's razor principle. Hence, we define the identifiability of LSSMs with sparse system matrices ensure that the resulting ambiguities can only be permutations and scaling, formalized below.
Definition 3.1. (Identifiabilitt) For LSSMs with nonzero system matrices A, B,C.C and D, if any nonsingular matrix & E R satisfying
| ll0 = A B
|Bllo and || ||C'llo.
(3)
must be a generalized permutation matrix, then such systems are said to be essentially identifiulle, up to permutation and scaling.
3.1 STUDENT'S L-DISTRIBUTION PRIOR
Here, we impose the Student's ' t-distribution prior on the system matrices A, B, C, and D to pro- mote model sparsity, because it can be sharply peaked at zero compared to other priors (Tipping, 2001). Generally, the Student's t-distribution prior is implemented in a hierarchical way (Wang al. 2024 Zhou et al. 202T). It imposes zero-mean Gaussian prior the system matrices and
3

---

Published as a conference paper at ICLR 2026
then adopts Inverse-Gamma distribution the unknown variance. For example, we can impose the Student's t-distribution prior on A to promote its sparsity as follows:
A 2T a.j
T
(Asg
Ta.y)
exp
(4)
2
=1
=1j=1
p(T.)
40 r P ((ao)
exp
bo
(5)
T a,
where r(.) is the gamma function, and Air and r a. are the ijth components of A and Fa. respec- tively. To generate con-informttive hyperprior on r Fa.ig ag and by are typically set very small values (e.g.. 10-9). In addition, Tb, Te, Td. I b.ii, I ee,ij and Ta,ú defined in similar man- ner (see Appendix A). For the process noise R and measurement noise Q, we impose a uniform distribution prior on them derive a flat prior.
3.2 LOSS FUNCTION
Following the Bayes' rule, we can combine the marginal likelihood and prior functions to derive the joint posterior distribution of all the unknown variables as follows:
p Y) x p(Y
X
((O))
(6)
manginal likelihood
prior
is the set of unknown variables and Y [y1, y2, 3/r] Because the system state 2, unobserved, it is hard to explicitly compute p(Y )). Hence, directly maximizing equation t estimate O is generally intractable. To tackle this problem, we explore the EM algorithm iteratively improve equation by regarding at as the latent variable. Instead of directly maximizing equation 6 the EM algorithm focuses on improving the expected value of the log posterior function of θ with respect to the state vector X [11 x2 .... æTl as follows:
where ④ = ,,
H Oh) = Ex [log (p(Y, X (e)(,
where X p X Y, O+), and X and O denote the estimates of X and O the kth iteration, respectively. It is known that iteratively maximizing equation 7 is equivalent to iteratively
maximizing equation 6 (Little & Rubin 2019).
3.3 EXPECTATION STEP: RAUCH-TUNG-STRIEES SMOOTHER
Because equation 7 involves (X] Y, OY, we first need to derive the conditional distribution of given Y and current *# {A Bk C D R Q* r IÍ r rt), which can be formulated as a classical smoothing problem. For LSSMs, the RTS smoother provides a closed-form smoothing solution for p (xx Y, O*).
Lemma 3.1. (RTS smoother
Särkkä & Svensson 2023) For LSSMs, the RTS smoother states that
p px; Y, @*) N P!))
(8)
where O.... T. Here, m and P are derived via the reverse-time recursions as follows:
m μ + G (mi+1 pl+1),
(9)
ph = > G (pL+ Eit1) (Gh)",
(10)
where G
Σ -Ei (At) (x*+1 The quantities pti phe Ej, and Σ coupled in equation
and equation T0
pre-computed using the Kalman filter follows:
A 11 - + Bius,
Σ = A >/ (At) + +R*
(11)
K = Σ (C#)' (C*** (O*) + Q+)
(12)
μ = k (Y: Ch D'u,), D/ (In kk:C Σ
(13)
4

---

Published as a conference paper at ICLR 2026
where I, is an identity matrix of dimension n. Note that the reverse-time recursions of equation and equation 10 start from the initial conditions m = μt and P = ET, and the recursions of 13 start from the mean pg and covariance 2 of the initial state æ0-
equation 11 equation
Besides p (xt Y, O*), we also need to derive the covariance matrix between the adjacent states 3t and 2t-1 given Y and @ to compute equation 7 gives necessary recursions.
To address this issue, the following lemma
Lemma 3.2. (The lag-one covariance smoother Särkkä & Svensson (2023) ) For LSSMs, the covari- ance matrix P between the adjacent states It and 2t-1 given Y and @ can be recursively derived as follows:
Pk fi,-1
E + G. Pt G. A*)) (G!
(14)
with Pk
TT--1 (I, C* A -T--
Based on Lemmas 3.1 and 3.2 we are able to calculate the loss function in equation 7 as follows:
H ( + O*) H1 (A, B, R D.F. D.C.C.D.....AAAA................
H.
(15)
where
H1 R) = xx (log pl A,
(16)
H2 D, Ex [log pl Y X, C, D.Q)).
(17)
b I Fa) = x [log (p( A, T. B, I (,,CC,T))(DE..)))
(18)
Due to the limited space, the detailed derivation of equation 15 and explicit mathematical expres- sions of (A, B, R), H2(C, Q). and H. (A B, B,C, D, Co F Fa) are given in Appendix
B.
3.4 MAXIMIZATION STEP: BLOCK COORDINATE DESCENT
Obviously, H (O @*) is a non-convex function and unknown variables are highly coupled. To provide analytical update formulas, we leverage the block coordinate descent method to sequentially optimize the model parameters.
Update procedures of A, B, C, and D. Due to the introduction of sparsity promoting priors, intractable derive closed-form solutions for A, B, C, and by setting the derivatives of H (@ @*) with respect to these variables to zero directly. To address this issue, we adopt row- wise update rule for A, B, C. and D. For example, the derivative of H the rth row of A, denoted as A,, is as follows:
@*) with respect to
a ( B R
Ha A, B Ch D* DK. D.
T
0A,
R
((m A,m B: (m 1) + (PI, A,P
A,
a,r
(19)
where R is the rrth component of R Pk ptt-1, 1.r. m and B are the rth rows of f,t- 1. m
and B&, respectively. In particular, Tkd diag a.r. with I²t a.r being the rth row of F Set- ting equation 19 to zero leads to the following update rule for A, at the kth iteration:
A k+1
(m B:u..) (m + Pkk Pi
P + m (m) Rh Rkd
(20)
5

---

Published as a conference paper at ICLR 2026
The detailed derivation of equation can be found in Appendix B.2 Similarly, update the rth row B, C, and at the kth iteration as follows:
B=
T
(m - + m) ui
Σ
T
ur + R2 I hid b,r
(21)
-1
C
Σ Uttr Ds (mit)
P + m (ml)) + Q² Frid Σ
(22)
op
D Ut.r c ui
Σ ueuy + Q Ftud
(23)
where yt. is the rth component of yes Qr is the rrth component of Q", and T b, are defined as that of T a,
and F kd
Update procedures of R and Q. The derivative of И ( @t) with respect to the rrth component
of R. denoted as R,,, is as follows (see Appendix B.3
0H, (Ah+1, R+1 1 R)
T
E
+ + II,
(24)
RRr
2
2R2,
where
π = m - A k+ m Bu+1
(25)
II = Pl A k+1 P P... P...-1 (A 1 + A+1 pk_ (A (A
k+1
(26)
with o. and II tor being the rth and rrth components of πł and II, respectively. Setting equa- tion 24 to zero yields the following update rule for Ryy the kth iteration:
2 + EE= II
(Tr)
(27)
T
Similarly, we can update the rrth component of Q, denoted as Q at the kth iteration as follows:
Q++1 Er (vtt )² + +I
1==1 w
(28)
T
where
/b = 1 m.
D un -Z= = --Ctlpp P. (Ch+l)
(29)
with vtt,, and V f being the rth and rrth components of wh and V respectively.
Update procedures of Ta, To, Pe, and F,, Because each component of Ta. Pb. To. and Ta is independent of the others, we can update them individually. For example, we can calculate the derivative of H ( @*) with respect to Ta,i at the kth iteration as follows:
OH, (A++1 B++1 C+++1 D T D P.
20p +
(A
(Ak+1) k++
1) 2bo
(30)
a..i
2r a..i
2T
a.1j
Setting equation 30 and solving for T Faji leads to:
(A + 20g
p a.ij
(31)
20p +3
Similarly, we can update each component of Tú, Te. and Fa at the kth iteration as follows:
1)²
2bo
(O 2 2bp
(D)+1)² 2bg
p++ b,gj
r++1 c,éf
r++1 d,i
(32)
2ap +3
20g + 3
2ap +3
Finally, we summarize the procedure for learning LSSMs with sparse system matrices Appendix

---

Published as a conference paper at ICLR 2026
3.5 GLOBAL CONVERGENCE ANALYSIS
Given @ the proposed learning algorithm presents the analytical mathematical expressions to de- define the proposed algorithm point-to-point mapping A(.). Lever- aging the Global Convergence Theorem (Luenberger et al. 1984), we can demonstrate that the proposed algorithm globally convergent.
rive @k+1 Hence,
Theorem 3.3. From any valid initialization point @º, the limit point of the sequence {0*) Ik=1 generated via @ = A (@*)
( is a local maximum (or saddle point) of equation 6
Proof. See Appendix D
SIMILARITY TRANSFORMATION OF LSSMS
The similarity transformation provides an equivalent realization original LSSMs by transforming states into different coordinate systems. For LSSMs, it is an important mathematical operation to analyze system properties like controllability, observability, and stability. Specifically, we can transform the state vector Dt into a new state vector 2; through the relation a; tr,, where & E Raxx is nonsingular matrix. As such, we can derive an equivalent realization of the original
LSSMs as follows (see Appendix E
It A + Bu, Et.
(33)
y: Ca. + Du: + Wt:
(34)
where
A = 1=
B= C=
(35)
and e, ~N,, (0, FR However, the similarity transformation makes it particularly difficult to accurately learn system matrices. Given the input signals {u,)-=1- the transformed L.SSMs can pro- duce the same output data {y:)7= as the original LSSMs. Hence, classical learning algorithms for LSSMs only learn the system matrices up similar transformation (Viberg 1994). For LSSMs with sparse system matrices, such a transformation changes not only the numerical values but, more importantly, the topological structure of the system matrices, resulting in misinterpretation of intrin- sic working mechanisms.
4.1 BENEFIT OF SPARSITY- PROMOTING PRIORS
Unlike classical learning algorithms, the proposed algorithm learns LSSMs with sparse system ma- trices by adopting sparsity-promoting prior to balance model complexity and modeling error Given the sparsity constraint of system matrices, the similarity transformation cannot be applied using any arbitrary nonsingular matrix &. For the identifiable LSSMs, the nonsingular matrix is restricted a generalized permutation matrix; otherwise, the transformed LSSMSs will include redundant parameters to describe the systems. For example, consider the LSSMs with sparse system matrices as follows:
0
0
0.87
2
0
0
2
0
0
0
0
=
0.8
0
B =
0
: C =
0
0
2 0
D
=
(36)
0
0.8
2
0
easy to verify that such system is identifiable, as the rank of system matrices is equal to the number of nonzero components. Hence, we can derive that the nonsingular matrix Φ must be one of the following generalized permutation matrices:
Fa
07 00000
a
07
a
000
07
0
00
a
0
a
0
0
a
000
6
0
b
0
b
(37)
000
0
0
0
where a, b, and c are arbitrary constants. As such, the transformed system matrices A, B, and C do not introduce additional parameters to increase model complexity.
Note that applying the similarity transformation with a generalized permutation matrix to the original state variables will scale their magnitudes and reorder them. However, it will scale the nonzero
7

---

Published as a conference paper at ICLR 2026
Dura 7.%E-
LBM. PEM
LSM HK
48ID
MLE
Ground Trush
MRE--10 44%
MRE--21 1.00
MRE-22.5R
MRE 88%
Time=240 14
Time=0.36
Time-18.45
Time-3.02 02
Time-152.344
0.0
0.5
0.6
0.0
0.
04
02
0.2
x110
0.5
0.5
1.5
100
50
0.5
15
1.5
0.5
05
0.5
0.5
Figure 1: Experimental results of all the algorithms the 10-dimensional synthetic system. Com- pared to the ground truth, only the proposed algorithm preserves the topological structure among the variables with + 27 In addition, the proposed algorithm obtains the lowest MRE among all the algorithms. Experimental results on the more complex non-diagonal and non-invertible cases can be found in Appendix I.I and Appendix I.2 respectively.
components and permute the rows columns of system matrices accordingly. Hence, an additional advantage of the pporsity prior is its ability to preserve the inherent topological structure among the variables. While the learned system matrices differ from the true ones in numerical values, such difference is only caused by the scaled definition of state variables, rather than failure to capture the underlying system dynamics.
EXPERIMENT
In this section, we validate the proposed algorithm on simulation and real-world datasets. addi- tion, we compare the proposed algorithm with classical ones mentioned previously to demonstrate its superior performance, including LSM PEM, LSM_HK, 4SID, and MLE. In experiments, the dataset is split into training and testing sets with a 2:1 ratio, where 66.7% of the data is used for training and 33.3% for testing. Due to the lack of ground truth for real-world datasets, we use the mean relative enror (MRE) to evaluate the performance of all the algorithms defined follows:
IIv- MRE = Σ 1=1 T
where { is the sequence of data points generated by the learned conducted using MATLAB
LSSMs in response to the same input signals. The experiments 2022b on the PC with an Apple MI Pro chip with 10-core CPU and 32GB RAM. Experimental re- sults on simulation and real-world problems illustrate that the proposed algorithm can preserve the inherent topological structure among variables and significantly improve prediction accuracy over the classical ones. However, because the closed-form updates entail the inversion of large number of matrices, the proposed algorithm involves high computational complexity (see Appendix F).
5.1 SYNTHETIC SYSTEMS
First, we test all the algorithms on a 10-dimensional synthetic system, indicating that each system matrix has 100 unknown parameters. Specifically. consider A to be anti- mitiodiaoonal matrix with nonzero components equal to 0.8, and set B C = D = 2I and R = Q = 0.81Inn- Because these system matrices are extremely sparse, accurately learning their topological structures is partic-
8

---

Published as a conference paper at ICLR 2026
Table 1: Comparison results the real-world industrial process systems
Dataset
Ours
LSM PEM
LSM HK
4SID
MLE
14.93%
17.90%
Inf
43.77%
20.14%
Evaporation system
(249.14 s)
(9.35 s)
(18.45 s)
(3.92
(152.34 s)
23.63% % (45.62
62.62%
31.21%
24.32%
30.21%
Glass furnace
(0.52 s)
(6.20
(0.29 s)
(33.36 s)
20.70%
22.70%
39.80%
29.26%
22.23%
Steam generator
(441.88 s)
(4.94 s)
(84.52 s)
(0.58 s)
(299.10 s)
ularly challenging. To generate data points, the value of zo is drawn from the Gaussian distribution with the mean diag [1 and identity covariance matrix, and the input signal is drawn from the uniform distribution on [0. 2]. As for algorithm implementation, we collect 2,100 data points and set the initial value of A, B, C, D, R. and Q both to be 10- Because the iterative numerical op- timization method rarely yields exact zeros, the learned parameters below the predefined threshold truncated zero accelerate convergence. The threshold selection procedure can be found in Appendix G To ensure fair comparison, the learned parameters of all the other algorithms that fall below this threshold are likewise set to zero.
Because the system matrices of the synthetic system are known, we can compare the learned system matrices with real ones directly. Figure T shows the learned system matrices of all the algorithms, together with the associated MRE and running time. Obviously, the learned system matrices of classical algorithms are completely different with the original ones in both numerical values and topological structures due to the similarity transformation, making it difficult for to understand the system. However, pparsity-promotinn transformation to be a generalized permutation matrix for this system. By comparing the learned B and C' of the proposed algorithm with real ones, we can derive O = 20 which indeed consistent with the analysis in Section 4.1 Hence, the learned system matrices of the proposed algorithm preserve the inherent topological structure among the variables, differing only in numerical values due to the scaled definition of state variables.
priors will restrict the nonsingular matrix of the similarity
We test the proposed algorithm on different initial values to demonstrate its robustness in Ap- pendix H.] We also conduct 20 independent trials report the success rates of all the algorithms in learning the topological structure of the system matrices and the average MRE in Appendix H.2 experimental results illustrate that the proposed algorithm stable across multiple runs.
In addition, we further demonstrate the effectiveness of the proposed algorithm on the non-diagonal and non-invertible cases in Appendix I.1 and Appendix I.2 respectively.
5.2 INDUSTRIAL PROCESS SYSTEMS
Next, we validate the proposed algorithm on the real-world datasets obtained from the Database for the Identification of Systems, which are standard datasets used for learning LSSMs (Zhu al. 1994 Martens 2010). While underlying physical systems may be non-linear, the learned LSSMs can provide a linear approximation of systems, enabling tractable analysis and control. To empirically ensure the performance of all the algorithms, we prefer to select datasets that contain at least 1 1,000 sample points and have multi-dimeessonal inputs and outputs.
Evaporation System. In industry, multiple-stage evaporators are widely used to reduce the water content product such as milk. The dataset is composed of 3-dimensional time-series with length of 6,305. The inputs consist of the feed flow, vapor flow to the first evaporator stage, and cooling water flow to the condenser, while the outputs include the dry matter content, flow rate, and temperature of the product.
Glass Furnace. The second dataset comes from the Philips glass furnace, which is used to melt materials into glass. The glass furnace has two burners and one ventilator. Hence, the dataset includes two heating inputs and one cooling input with a length of 1,247. In addition, we collect three outputs from temperature sensors in cross section of the furnace.

---

Published as a conference paper at ICLR 2026
Steam Generator. The dataset comes from a boiler at Abbott Power Plant in Champaign IL, which a dual-fuel (oil and gas) combustion unit used for both heating and electricity generation. consists of four imputs (i.e., air flow rate, fuel flow rate. feedwater flow rate, and disturbance) and four outputs (i.e., drum pressure, excess oxygen, drum water level, and steam flow) with a length of 9,600.
Table displays the MRE between the predicted outputs of all the learned LSSMs and the real one, and records the running time of each algorithm. Due the lack of ground truth, the learned system matrices of all the algorithms are not depicted for comparison. Because the proposed algorithm obtains minimum MRE across all the datasets, experimental results demonstrate its superiority over the classical ones in real-world applications. Note that the hidden states can be defined in different coordinate systems as discussed in Section + As a result, the learned system matrices of classical algorithms may differ from the ground truth in both numerical values and topological structures. However, the proposed algorithm learns the system matrices by balancing model complexity and modeling error. Unlike classical algorithms, it will thus learn LSSMs with only the minimally required parameters to explain the time-series data. In particular, if the real-world systems are identifiable, the proposed algorithm preserves the topologiaa structure among variables, which is more valuable for exploring interaction laws of systems.
6 DISCUSSION
To learn the LSSMs with sparse system matrices, we impose matrices to balance model complexity and modeling error. Following the MAP principle, we then learn system matrices by exploring the EM algorithm to maximize the joint posterior distribution composed of the priors and marginal likelihood function. Based on the Global Convergence The- orem, we demonstrate that the sequence generated by the proposed algorithm converges local maximum or saddle point of the posterior distribution. In addition, we explain why the sparsity- promoting prior capable of retaining the inherent topological structure of LSSMs, as the non- singular matrix of the similarity transformation is limited to be a generalized permutation matrix. Hence, the proposed algorithm is more useful for us explore the interaction laws of LSSMs com- pared to the classical ones.
priors system
There still remain some potential limitations associated with the proposed algorithm. First, the similarity transformation may shrink many parameters in system matrices to very small values, potentially leading to numerical errors. However, note that is a common issue faced by all the learning algorithms for LSSMs. The other limitation is that the proposed algorithm is hard to deal with large-scale problems currently due to its high computational requirements. Hence, our future work will focus reducing computation time make the proposed algorithm applicable to large- scale settings. For instance, we can adopt the stochastic EM algorithm to reduce the computational cost by using a mini-batch of data instead of the full batch during the expectation step (Chen al. 2018). In addition, we hope to explore how to ensure that the learned system matrices are exactly the same as the true ones by imposing additional constraints on system matrices beyond the sparsity constraint. Overall, we believe that the proposed algorithm sheds light on the learning of LSSMs with sparse system matrices.
ACKNOWLEDGMENTS
This work is supported by the National Natural Science Foundation of China (62503185, 62303118) and Doctoral Foundation of Fuyang Normal University (2021KYQD0034).
ETHICS STATEMENT
This paper presents theoretical research whose goal is to advance the field of Machine Learning. Because it does not involve human participants, animals, sensitive personal data, or other foresecable ethical concerns outlined in the ICLR Code of Ethics, no specific ethical issues arise from this paper.
10

---

Published as a conference paper at ICLR 2026
REPRODUCIBILITY STATEMENT
To ensure the reproducibility of this research, the main text and Appendix provide detailed theo- retical derivations and proofs of the proposed algorithm. Source codes are available on GitHub at https: / com/ Art.inttttht/// LerannssLeannnnn-Sppaaiaee datasets used in our experiments are publicly available https::////// wuthttss//hhoeee hww.e//h/eeeeeeee esat kuleuven. b da /sy
LSSM Moreover, the real-world
REFERENCES
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.
David Belanger and Sham Kakade. linear dynamical system model for text. In International Conference Machine Learning, pp. 833 842. PMLR, 2015.
Steven L Brunton, Joshua L Proctor, and Nathan Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences, 113(15):3932-39)) 2016.
Emmanuel J Candes, Michael B Wakin, and Stephen P Boyd. Enhancing sparsity by reweighted & minimization. Journal of Fourier Analysis and Applications, 14:877-905, 2008.
Nilanjana Chakraborty, Kshitij Khare, and George Michailidis. A Bayesian framework for sparse estimation in high-dimensional mixed frequency vector autoregressive models. Statistica Sinica, 33:1629-1652, 2023.
Jianfei Chen, Jun Zhu, Yee Whye Teh, and Tong Zhang. Stochastic expectation maximization with variance reduction. In Advances in Neural Information Processing Systems, volume 31, 2018.
Yanxi Chen and Vincent Poor. Learning mixtures of linear dynamical systems. In International Conference on Machine Learning, pp. 3507-3557. PMLR, 2022.
Kevin Course and Prasanth B Nair. State estimation of a physical system with unknown governing equations. Nature, 622(7982):211-26,, 2023.
Yonathan Efroni, Sham Kakade, Akshay Krishnamurthy, and Cyril Zhang. Sparsity in partially controllable linear systems. In International Conference on Machine Learning, pp. 5851-5860. PMLR, 2022.
Mohamad Kazem Shirani Faradonbeh, Ambuj Tewari, and George Michailidis. Finite time identifi- cation in unstable linear systems. Automatica, 06:342-353, 2018.
Mohamad Kazem Shirani Faradonbeh, Ambuj Tewari, and George Michailidis. Optimism-baa based adaptive regulation of linear-quadratic systems. IEEE Transactions Automatic Control, 66(4): 1802-1808, 2020.
Wouter Favoreel, Bart De Moor, and Peter Van Overschee. Subspace state space system identifica- tion for industrial processes. Journal of Process Control, 10(2-3):149-155, 2000.
Pierre Garrigues and Bruno Olshausen. Group sparse coding with a Laplacian scale mixture prior. In Advances neural information processing systems, volume 23, 2010.
Zoubin Ghahramani and Geoffrey E Hinton. Parameter estimation for linear dynamical systems. 1996.
Stuart Gibson and Brett Ninness. Robust maximum-likelihood estimation of multivariable dynamic systems. Automatica, 41(10)::667-1682, 2005.
Elad Hazan, Singh Karan, and Zhang Cyril. Learning linear dynamical systems via spectral filtering. In Advances in Neural Information Processing Systems, volume 30, 2017.

---

Published as a conference paper at ICLR 2026
Jiabao He, Ingvar Ziemann, Cristian R Rojas, and Häkan Hjalmarsson. Finite sample analysis for a class of subspace identification methods. In IEEE Conference on Decision and Control, pp- 2970-2976. IEEE, 2024a.
Xin He, Yasen Wang, and Junyang Jin. Bayesian inference and optimisation of stochastic dynamical networks. International Journal of Systems Science, pp. 1-15, 2024b.
Wenbing Huang, Lele Cao, Fuchun Sun, Deli Zhao, Huaping Liu, and Shanshan Yu. Learning stable linear dynamical systems with the weighted least square method. In International Joint Conferences on Artificial Intelligence, pp. 1599-1605, 2016.
Junyang Jin, Ye Yuan, and Jorge Gonçalves. A full Bayesian approach to sparse network inference using heterogeneous datasets. IEEE Trunsactions on Automatic Control, 66(7):3282-2288, 2020a.
Junyang Jin, Ye Yuan, and Jorge Gonçalves, High precision variational Bayesian inference of sparse linear networks. Automatica, 118:109017, 2020b.
Tohru Katayama et al. Subspace methods for system identification, volume 1. Springer, 2005.
Wallace E Larimore. Canonical variate analysis in identification, filtering, and adaptive control. In IEEE Conference on Decision and Control, pp. 596-559644 IEEE, 1990.
Roderick JA Little and Donald B Rubin. Statistical analysis with missing data, volume 793. John Wiley & Sons, 2019.
Lennart Ljung. Prediction error estimation methods. Circuits, Systems and Signal Processing, 21: 11-21, 2002.
David G Luenberger, Yinyu Ye, et al. Linear and nonlinear programming, volume 2. Springer, 1984.
Bethany Lusch, Nathan Kutz, and Steven L Brunton. Deep learning for universal linear embeddings of nonlinear dynamies. Nature Communications, 9()::995,, 2018
Giorgos Mamakoukas, Maria Castano, Xiaobo Tan, and Todd Murphey. Local Koopman operators for data-driven control of robotic systems. In Robotics: Science and Systems, 2019.
Giorgos Mamakoukas, Orest Xherija, and Todd Murphey, Memory-efficiee learning of stable linear dynamical systems for prediction and control. In Advances in Neural Information Processing Systems, pp. 13527-13538, 2020.
James Martens. Learning the linear dynamical system with ASOS. In International Conference Machine Learning, pp. 743-750. Citeseer, 2010.
Aditya Modi, Mohamad Kazem Shirani Faradonbeh, Ambuj Tewari, and George Michailidis. Joint learning of linear time-invariant dynamical systems. Automatica, 164:111635, 2024.
Samet Oymak and Necmiye Ozay. Non- asymptooi identification of Iti systems from single trajec- tory. In American Control Conference, pp. 3655-5661. IERE, 2019.
Wei Pan, Ye Yuan, Jorge Gonçalves, and Guy-Bart Stan. A sparse Bayesian approach to the identifi- cation of nonlinear state-space systems. IEEE Transactions on Automatic Control, 61(1):182-187, 2015.
G Pillonetto and Lennart Ljung. Full Bayesian identification of linear dynamic systems using stable kemels. Proceedings of the National Academy of Sciences. 120( 120X18):2218819711 2023.
S Joe Qin. An overview of subspace identification. Computers Chemical Engineering, 30:1502- 1513, 2006.
Srijata Samanta, Kshitij Khare, and George Michailidis. A generalized likelihood-based Bayesian approach for scalable joint regression and covariance selection in high dimensions. Statistics and Computing, 32(3):47, 2022.
12

---

Published as a conference paper at ICLR 2026
Simo Särkkä and Lennart Svensson. Bayesian filtering and smoothing, volume 17. Cambridge University Press, 2023.
Robert H Shumway and David S Stoffer. An approach to time series smoothing and forecasting using the EM algorithm. Journal of Time Series Analysis, 3(4):253-264, 1982.
Gavin Smith, João de Freitas, Tony Robinson, and Mahesan Niranjan. Speech modelling using subspace and EM techniques. In Advances in Neural Information Processing Systems, volume 12, 1999.
Robert Tibshirani. Regression shrinkage and selection via the Lasso. Journal of the Royal Statistical Society Series B: Statistical Methodology, 58(1):267-22) 1996.
Michael E Tipping. Sparse Bayesian learning and the relevance vector machine. Journal of Machine Learning Research, 1:21 1-244, 2001.
Jack Umenberger, Johan Wägberg, Ian R Manchester, and Thomas B Schön. Maximum likelihood identification of stable linear dynamical systems. Automatica, 96:280-292, 2018.
Peter Van Overschee and Bart De Moor. N4SID: Subspace algorithms for the identification of combined eeterministic-stochaeee systems. Automatica, 30(1):75-9))) 93, 1994.
Peter Van Overschee and BL De Moor. Subspace identification for linear systems: The- Springer Science & Business Media, 2012.
org-lmppeeennationn-Apiiiicc
M Verahegen and Patrick Dewilde. Subspace model identification. Part i: The output-error state- space model identification class of algorithm. Internation Journal of Control, 56:1 16:1187-1210, 1992.
Mats Viberg. Subspace methods in system identification. IFAC Proceedings Volumes, 27(8):1-12, 1994.
Jin Wang and S Joe Qin. A new subspace identification approach based on principal component analysis. Journal of Process Control, 1288):841-855 2002.
Yasen Wang, Junlin Li, Zuogong Yue, and Ye Yuan. An iterative Min-Min optimization method for sparse Bayesian learning. In International Conference on Machine Learning, volume 235, pp- 50859-50873, 2024.
David Wipf and Srikantan Nagarajan. A new view of automatic relevance determination. In Ad- vances Neural Information Processing Systems, volume 20, 2007.
David P Wipf and Bhaskar D Rao. Sparse Bayesian learning for basis selection. IEEE Transactions on Signal Processing, 52(8):2153-)144 2004.
CF Jeff Wu. On the convergence properties of the EM algorithm. The Annals of Statistics, pp- 95-103, 1983.
Xiaofeng Yuan, Yalin Wang, Chunhua Yang, Zhiqiang Ge, Zhihuan Song, and Weihua Gui. Weighted linear dynamic system for feature representation and soft sensor application in nonlin- ear dynamic industrial processes. IEEE Transactions on Industrial Electronics, 6582)11508-1517, 2017.
Ye Yuan, Xiuchuan Tang, Wei Zhou, Wei Pan, Xiuting Li, Hai-Tao Zhang, Han Ding, and Jorge Goncalves. Data driven discovery of cyber physical systems. Nature Communications, 10(1): 4894, 2019.
Wei Zhou, Hai-Tao Zhang, and Jun Wang. An efficient sparse Bayesian learning algorithm based on Gaussian-scale mixtures. IEEE Transactions Neural Networks and Learning Systems, 33(7): 3065-3078, 2021.
Yucai Zhu, Peter Van Overschee, Ban de Moor, and Lennan Ljung. Comparison three classes of identification methods. IFAC Proceedings Volumes, 27(8):169-111 1144, 1994.
13

---

Published as a conference paper at ICLR 2026
A SPARSITY-PROMOTING PRIOR
Besides A, we also impose the sparsity-promont priors on B, C, and D as follows:
P
n To,ay)
B
P rb) =
(By
P
cxp
(38)
2TT
2T b,új
1=
<= 1=1
C2
p(C r.)
P
CCy Tea))
exp
(39)
2TTe. c,
2T c,
i=1j=1
=1j=1
P
P
D a.i
P PJ)
P (D. (Ta.)
exp
(40)
2 d.j
==1j=1
where To r eag and Fa.y are the ijth components of To, Fe: and Fa, respectively. To complete the hierarchy, the Inverse-Gammi distribution prior imposed on each component of Tb. and Ta as follows:
pTr))
P
dg
he
du-1
exp
by
(41)
((ao)
To,ij
bp
(TT))
((ao)
c
exp
Te,ij
(42)
bg 4
bo
a(Ta)
r
00
11
exp
(43)
((ag)
T
14

---

Published as a conference paper at ICLR 2026
B DETAILED MATHEMATICAL DERIVATION
B.1 DERIVATION OF EQUATION 15
Given the conditional independence between the variables, derive
H O )*)
x [log (p( Y, X O ())((())) [))]
E x llog (P(Y)
X, (X Q ()))
x Dlg (p(Y)
X, C, X A, B,
x [log (p(Y
X,C,
A,B,
E [log p X
A, R)
E po X.C.D.Q)
H)(A,B.R)
(.D.))
Ex [log (p(A, P P B, D.)
(44)
H.L.BBCC..D............
Explicit mathematical expression of Hi B, R). Based on equation and the chain rule in probability, we can derive
( A. R.
p(zo)
I p(zz xt-1, ABBRR)
T x [[ 1=1
cxp (- (_(x) A
Bu. RR
A Bu))
(45)
2
Hence,
H, R)
= -xx [lgg((
p( X
A, B, R))
Tlog
R +E (x: Ax;-1 - Bu:) RR (xx - Ax Bu:)
=Ex²
2
Tlog
R
T Dt= Exi [xx Ax,-1
Bu;) R. 1 Ax.-1 Bui)
2
Tlog
R
ΣT
(m
Am
-Bu)) R-1 {m, Am -Bu)
2
Σ [Tr
(R pt) Tr (R-1
APk
Tr (r''
( + (AAR--AAP...))
AAPL1)
2
(46)
Explicit mathematical expression of H2,[C, D, Q). Based equation 2 we can derive
p( X.C.D.Q)
T
II p(u: Et. C.D)
1=1
Q - exp (_ (y) Ca.
Du.YQ-' - Cr: Dus))
[[ 1=1
(47)
2
15

---

Published as a conference paper at ICLR 2026
Hence,
H,(C,D,Q)
1Exxllogpp
X, X, C, D,Q)
Ex* Tog
Q + ++= (y/
Cx. Du.)Q Q Cx, Du,)
2
Tlog Q + Σi= Ex (y: Cx.
Du, 1 Ca. Due)
2
Tlog
Q + ET
((y) Cml -u, (y: Cm Du.) TT(C-
2
(48)
Explicit mathematical expression of H3 T, Based on the priors im- posed the system matrices and corresponding hyperparameters, derive:
P(A.T ( B, Dd)
P(A T (T ( (T)p(T (P)(D))(DD Ta)p(Ta)
70
11
A + 2bo
n
P
B + 2bu
a.
exp
2T
exp
2T
Cl 2bg
20p
D 260
x
C
exp
2D c,
exp
2r 2T
(49)
=1
=1
Hence, we have
D.A.AB.C.D.Da.I..........
Ex [log ((A, p B, Co)p(CC
= Exx [log (p(A I. TT)p(T (P To)p(TT)p(C P.)p(T.)((D [a)e(T.))
(2ap + 3) log
Ta.i5
A?₁ + 260
2
2T a.55
(2ag + 3) log
To,ij
B3 + 260
2T b.5
(2a0 + 3) log
Tos
C2 2ho
2r c.ij
(2ao + 3) log
Ta.9
D3, + 2bo
(50)
B.2 DERIVATION OF EQUATION 19
Because H, (A, B, R) and H2 (A, B. C. a. P. Pa, involve A, we can derive the term
related A, the kth iteration as follows:
H1 A, k Rd + H3 (A Ck DE D r
Σx Σ (R..)) (Tr 2A, (P%)-1,)) Er=1 A.F.L
2
ΣF Σn=
R -1 (m
A, -Buuu) (m,r - A - Blus)
+c, (51)
16

---

Published as a conference paper at ICLR 2026
where è is the term independent of A,. Hence, we can calculate the derivative of H @*) with respect to A, as follows:
@ )*)
0A,
0H, (A, B², R²)
0Ha (A, B* Oh D l N h.r
at
Ö.,,
oA,
T -1 (P + (m, A,m B) (m,-)) A, A,
R
A, a,r (52)
Setting equation 52 to zero leads to
T
Ar k:
(m-_1)= + Pİ-1) + A, AR a.
Σ
pk Pi,-rr + (mirr - B: (m-_1))
1=1
t=1
(53)
Hence, we can update A, at the kth iteration as follows:
T
A h+1
Σ
(m - B) (m/-1) + P P..-1.))
X
P + m (m + Re rkd
air
(54)
B.3. DERIVATION OF EQUATION 24
Given that R is a diagonal matrix, we can re-express H1 (A*+1, BB++ R) as follows:
H1 (A k+1 Bk+1 R
T T'" 2r=1 log Rer
r=l
21=1
R
(.r)
Σ Σ t=
IT
Gr
(55)
2
where
k = m k A k+
m B++ B
u4:
(56)
II: = P k+1
pk 61-1
P...- (A++))' + A ph_ (Aht1)'
(57)
with ,Êi and II err being the rth and rrth components of कड and II. respectively. Hence, can calculate the derivative of H(@) +*) with respect to R,, as follows:
0H(@ 10 OH, (A++1 k+1 A R) RR,
T Σt ET, (TT))
ΣT II
(58)
OR
2R,,
2R2
2R
Setting equation 58 to zero leads to the update rule for R,y at the kth iteration as follows:
Σt (T )² + -T T IIA
R
(59)
T
17

---

Published as a conference paper at ICLR 2026
C PSEUDOCODE FOR LEARNING LSSMS WITH SPARSE SYSTEM MATRICES
Algorithm The proposed algorithm for learning LSSMs with sparse system matrices
Input: Time-series data {(u,, 8t)) initial guess of O and maximum number of iterations
kmaz
for k to kmaz do
for t = 1 to T do
Update the mean m and covariance P. of ot via equation and equation 10
respectively
Update the covariance P between 2e and # via equation
end for
Update system matrices A, B, C. and D via equation 20 equation 23 respectively
Update noise covariance matrices and Q via equation (7 and equation 28 respectively
Update hyperparameter matrices Fa. Pi, Te, and Pa via equation 31 and equation tively
32 respec-
if stopping criterion is satisfied then Break end if end for
Output: System matrices A, B, C, and D, noise covariance matrices R and
D PROOF OF THEOREM 3.3
Proof. To illustrate the global convergence of the proposed algorithm, need to demonstrate that it satisfies the three necessary conditions required by the Global Convergence Theorem (Luenberger
al. 1984). For A(-) and Ω {O Vep ( Y) 0}, & (() -p(O Y) descent
function. For the point @ in Ω. it is straightforward to conclude that L (@ k+1)
< C (@*)
following the basic property of the EM algorithm. For the point @ outside 52, H +*) is continuous in both arguments,
have C (@*+1) L (+)) (see Theorem 2 in Wu 1983 Hence,
()) is descent function for 52 and A(-). The sequence {0] is contained in a compact set. If any component of O is unbounded, (@) tends infinity. Given L (@*+1)
(e*), there must exists a compact set containing the sequence { ]= The mapping A(.) is closed at points
outside Ω. Because the proposed algorithm updates O++1 from @ via analytical mathematical expressions, and all elementary functions involved are continuous, A(-) can thus be regarded
continuous function. addition, A(-) is point-to-point mapping. Hence, A(.) is closed at points outside S (see Example 2 in Section 7.6 in Luenberger al. (1984) The proof is completed. O
18

---

Published as a conference paper at ICLR 2026
E EQUIVALENT REALIZATION OF LSSMS
Based on the transformed coordinates, we can derive
F 3e. A + Bu, + te: ($ .1 + ( u, + DEr:
(60)
1/: Ca, + Du, + w, = (c& (C Zt + Du, + w.
(61)
Hence, we can derive an equivalent realization the original LSSM as follows:
ã = A + Bu, + Et.
(62)
C3. + Du; + Wt:
(63)
where A -A $B. C C$ -1 and ε DE Because Er N (O, R), we can derive the mean of Eg as follows:
E E = +E, DE[s,] = O.
(64)
In addition, the covariance of e, can be derived as follows:
[(E, E( (E, EEE))) E [ DE[e,, DRD
(65)
Hence, we have E N (0, R
COMPUTATIONAL COMPLEXITY ANALYSIS
The high computational cost of the proposed algorithm primarily stems from the matrix inversion operations involved in closed-form updates. In each iteration, the expectation step entails the inver- sion of matrices of size n and m m total of T times following Lemma 3.1 Hence, the computational complexity of the expectation step is mainly determined by OT 3 + Tmm). In the maximization step, we derive row-wise update rules for A, B, C, and D enable analytical update formulas. To update each row of A, B. C, and D, we need to calculate the inverses of the and p x matrices as shown in equation 20 equation 23 As a result, the computational complexity of the maximization step is dominated by ( (n p²)). Because the computational cost of the proposed algorithm scales at least cubically with respect to one of m, n, and p. it is hard to deal with the large-scale problems currently. Hence, our future work will explore how mitigate the computational bottleneck to make it more applicable.
19

---

Published as a conference paper at ICLR 2026
G THRESHOLD SELECTION
While the proposed algorithm is globally convergent, the generated sequence converges to a local maximum or saddle point of the posterior distribution only in the limit of infinite sequence length. Hence, many learned parameters in system matrices approach zero but never reach it exactly during algorithm implementation. To achieve accurate topology recovery, it is necessary set the learned parameters below predefined threshold to zero to avoid numerical errors. In particular, the learned parameters corresponding to the true zero components in the system matrices are typically several orders of magnitude smaller than those learned for the nonzero components. Hence, this pronounced difference magnitude provides wide and stable range for selecting the threshold.
Specifically, the threshold can be selected by visualizing how the number of nonzero components in learned system matrices, denoted as N, varies with the threshold. When the threshold is too small, all components of learned system matrices remain nonzero, N is thus relatively large. As the thresh- old increases, the learned parameters corresponding to the true zero components are gradually set to zero, and N decreases accordingly. However, due to the pronounced difference in magnitude be- tween the learned parameters corresponding to the true zero components and those corresponding to the nonzero components, N remains stable old exceeds this range, N begins to decrease again because some learned parameters corresponding to nonzero components are also set to zero. Therefore, the threshold can be safely selected from this stable interval to ensure reliable topology recovery. For example, Figure 2 illlustrates the number of nonzero components in the learned system matrices under different threshold settings for the synthetic system in Section [.1 The orange markers indicate the threshold values at which the proposed algorithm successfully learns the inherent topological structure among the variables, and the MRE remains unchanged. As such, we set the threshold to 0.005 in the experiment.
certain range of threshold values. Once the thresh-
400
350
300
Number of nonzaraa
250
50
10
10
10
10
10
Threshold
Figure 2: The number of nonzero components in the learned system matrices across various thresh- old settings. When the threshold lies within the interval marked by the orange dots, the proposed algorithm can successfully recover the inherent topological structure among the variables, and the MRE remains unchanged.
20

---

Published as a conference paper at ICLR 2026
H ADDITIONAL EXPERIMENTAL RESULTS ON THE SYNTHETIC SYSTEM
H. EXPERIMENTAL RESULTS UNDER DIFFERENT INITIAL VALUES
To assess the sensitivity of the proposed algorithm to initial values, we further evaluate its perfor- mance on the synthetic system described in Section 5.1 by initializing the A, B, C, D, R. and Q as 1o where a varies from 0.6 1.4 increments of 0.2.
Table 2 records whether the proposed algorithm successfully learns the topological structure of the system matrices and MRE defined in the paper. As observed from the table, the experimental results are consistent with those in Section [.T] indicating that the proposed algorithm is robust to initial values. Remarkably, even when the imitial state transition matrix is unstable (i.c. a > 1), the proposed algorithm is still able to accurately learn the topological structure of the true system.
Table 2: Experimental results of the proposed algorithm the 10-dimensional synthetic system across different initial values
a
0.6
0.8
1.2
1.4
Success?
v v
V
V
MRE
5.70%
5.70%
5.70%
5.70%
5.69%
H.2 EXPERIMENTAL RESULTS OF INDEPENDENT TRIALS
Here, conduct 20 independent trials demonstrate that the proposed algorithm is stable across multiple runs. Specifically, we set the random seed to increase evenly from 1 to 20. Table reeport the success rates of all the algorithms in learning the topological structure among the variables and the average MRE. Compared to the classical algorithms, only the proposed algorithm successfully learns the inherent topological structure among the variables in almost all cases. In addition, the proposed algorithm can achieve 100% success rate by slightly increasing the threshold below which learned parameters are set zero.
Table 3: Experimental results of all the algorithms on 20 independent trails
Method
Ours
LSM PEM LSM HK
4SID
MLE
Success rate
95 %
0 19.63%
0
0
0
Average MRE
5.67 %
21.56%
20.37%
5.97%
21

---

Published as a conference paper at ICLR 2026
I EXPERIMENTAL RESULTS ON THE NON-DIAGONAL AND NON-INVERTIBLE SYNTHETIC SYSTEMS
1.1 NON-DIAGONAL SYSTEM
Here, we further test the proposed algorithm on a non-diagonal system to illustrate its effectiveness. To generate non-diagonal system matrices, we randomly set one component per row to 0.8 in A, and to in B, C, and D, with all other elements set to zero. Particularly, the nonzero elements are deliberately placed to ensure that A, B, C, and maintain full rank. All other experimental settings remain the same in Section 5.1
Ours
LSW PEM
LBM HK
4510
MLE MRE-5- .2%
MRE--d 58% Ground Truth
MRE-11p 78%
MRE-21 .2%
MRRE-19.81 Timent 30
lmme22st.. 19
Timeer 7.00
1mme17.22 17.22
Times 151 17
05
66
0.2
405
1.5
0.5
405
1.5
Figure 3: Experimental results of all the algorithms on the non-diagonal systems.
Figure records the experimental results of all the algorithms the non-diagonal system. Due the similarity transformation as discussed in Section 4 the learned A, B, and C' of all the algorithms differ form from the ground truth. Unlike classical algorithms. however, the proposed algorithm restricts the nonsingular matrix of the similarity transformation to be a generalized permutation matrix. By comparing the learned B and C' with ground truth, we can derive the nonsingular matrix as follows:
0
0
0
2
0
0
07
0
0
0
2
0
0
0
0
0
2
0
0
0
0
0
0
2
S
(66)
0
0
0
0
0
0
0
0
0
0
0
0000
2
0
0
2
0
0
0
0
0
2
0
0
0
0
This implies that the discrepancy between the learned system matrices of the proposed algorithm and ground truth stems from different ordering and scaling of system states. Hence, it is easy to check that the proposed algorithm still accurately captures the underlying system dynamics and preserves the inherent topological structure among the variables.
22

---

Published as a conference paper at ICLR 2026
I.2 NON-INVERTIBLE SYSTEM
To evaluate the proposed algorithm on the non-invertible system, we replace the system matrix A in the synthetic system of Section 5.1 with singular matrix. Specifically, first construct an anti- diagonal matrix with all nonzero components set to 0.6, and then replace its seventh and eighth rows 0.6 0.6 O 000000 0J. As such, the seventh and eighth rows of A are
with the following vector: identical, making the matrix non-im -on-invertible. Besides, all other experimental settings remain the same as in S.c reports the experimental results of all the algorithms on the non-invertible system. Similarly, we can derive $ = 2T by comparing the learned system matrices B and C with the ground truth. Hence, the proposed algorithm still successfully preserves the topological structure of system matrices, and obtains the lowest MRE compared to the other algorithms.
LSM
LBM NRE--8899' Time-22:11
481D
Ground
18%
MREE-15 28%
MREE-20 92%
T5me-222
Time=11
0-5
0.5
044
02
.0.
0.5
1.5
-100
1.8
Figure 4: Experimental results of all the algorithms on the non-invertible system.
23

---

Published as a conference paper at ICLR 2026
J THE USE OF LARGE LANGUAGE MODELS (LLMS)
In accordance with the ICLR policy on the use of large language models (LLMs), we declare that LLMs are only employed writing assistant polish the language of this paper. They do not contribute to the research ideation, experimental design, data analysis, or interpretation of the results.
24


Reviews:
[{"review_id": "FqwfA09xeq", "reviewer": "Reviewer_ttg7", "summary": "This paper proposes a new EM algorithm for learning linear state-space models (LSSMs) with explicitly sparse system matrices, a critical aspect for interpretability and efficiency in various scientific domains. The authors impose sparsity-promoting Student’s t-distribution priors on system matrices and develop a principled MAP estimation framework. The resulting EM procedure leverages the Rauch-Tung-Striebel smoother in the E-step and block coordinate descent in the M-step, with global convergence to local MAP optima proven via the Global Convergence Theorem. The paper claims its approach avoids the notorious similarity transformation pitfalls of classical LSSM learning algorithms, thus better recovering the true sparse topology. Experiments on synthetic and several real-world benchmark datasets show improvements in prediction accuracy and topological interpretability compared to established baselines.", "questions": "I would wish author(s) can clearify my concerns around thresholding for sparsity, hyperparameter selection, and comparison to Lasso/regularized MLE-style estimators.\n\nThe constraints on similarity transformations in the presence of sparsity-promoting priors are presented for special cases. Can these insights be generalized, or are there pathological systems where the identifiability argument fails?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper demonstrates methodological rigor through a principled MAP-EM formulation for LSSMs under sparsity priors, with complete derivations, update steps, and connections to established statistical techniques. Its theoretical soundness is reinforced by a formally proven convergence guarantee via the Global Convergence Theorem. The work explicitly addresses the similarity transformation issue, ensuring identifiability and interpretability by restricting transformations to generalized permutation matrices. The experiments on synthetic and real-world industrial datasets validate the method’s effectiveness, with figures and tables clearly illustrating preservation of topological structures and improved quantitative performance. The paper’s mathematical transparency and detailed derivations enhance reproducibility and clarity.", "weaknesses": "- The algorithm's reliance on full EM with smoothing and coordinate updates (Algorithm 1) and complexity discussed in Section 6, makes it computationally intensive. The paper only briefly mentions the limitation with respect to large-scale settings, without offering profiling, timing, or strategies for acceleration (e.g., stochastic EM, parallelization).\n\n- As stated in Section 5.1 and experimental appendices, sparsity is only enforced \"up to a predefined threshold\" with truncation to accelerate convergence. While this is common in Bayesian learning, an explicit analysis of how the choice of this threshold affects success rates, stability, and meaningful recovery of structures is not present, nor is there a discussion of methods that enforce exact zeros by construction.\n\n- The argument that sparsity-promoting priors restrict similarity transformations to generalized permutation matrices is supported only for very specific classes of topologies; most of the argument is illustrated rather than proven generally. Real-world systems with denser interconnections may still pose identifiability challenges; this should be discussed further.\n\n- While the paper is mathematically rigorous, much of the machinery (Student’s $t$ prior for sparsity, hierarchical Bayesian treatment, EM updates, block coordinate descent, even the use of convergence theorems) is largely a synthesis/extension of established sparse Bayesian learning and system identification literature. The main differentiating insight (identifiability under sparsity priors limiting transformations to permutation/scaling) though interesting, is not deeply explored beyond special-case analysis.\n\n- Although synthetic and 3 real-world datasets are experimented on, comparison to more recent or advanced sparse LSSM learning methods, such as those using Lasso or regularized MLE approaches (e.g., [1]), is missing. Moreover, evaluation is limited mostly to MRE; the case for learning the \"correct\" topology is largely qualitative except for synthetic scenarios, which could be expanded quantitatively (e.g., with precision/recall on support recovery).\n\n**Refs:**\n\n[1] Fattahi, Salar, Nikolai Matni, and Somayeh Sojoudi. \"Learning sparse dynamical systems from a single sample trajectory.\" 2019 IEEE 58th Conference on Decision and Control (CDC). IEEE, 2019.", "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": "O97Uz24Oek", "reviewer": "Reviewer_ggGN", "summary": "This paper introduces a method to learn sparse linear state-space models. This is achieved by a hierarchical Student's t-distribution prior on system matrices which results in sparsity bias. The system matrices, hyperprior parameters and noises are then iteratively optimized via an expectation-maximization procedure. For each minimization step, the authors derive closed-form update rules based on block-coordinate descent on each parameter. They validate their results on some synthetic and real-world systems.", "questions": "1. Can the authors discuss technical novelties (if there is any) compared to the related work?\n2. Can you discuss the computational requirements a bit more in depth compared to other algorithms? Isn't it normal that each iteration takes $T$ time as that's the size of the input. Why is it more computationally expensive? Is it the inverse operations in the closed-form updates?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper seems to be a natural extension of prior work on sparse bayesian learning (e.g., of the ideas in Wang et al., 2024). I am not sure how original the technical ideas are but the paper is well-executed and derives the necessary results in a more complicated setting.\n\nThe technical results are complete in the sense that there seems to be nothing left on the table given the introduced setting. Additional improvements can only be made by imposing other constraints or relaxed versions of the algorithm for computational efficiencies which I find out of the scope. Therefore, I consider the theoretical results complete and good quality. The ideas and the results are presented very clearly.", "weaknesses": "1. The results are promising empirically but more can be done if the authors think the algorithm presented is practical. For example, they can include runtime comparisons.\n2. I think the paper puts the right results together from previous literature and has rather minimal technical contributions.", "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": "wL1Lhmu8yq", "reviewer": "Reviewer_jpHJ", "summary": "This paper addresses the problem of learning linear state-space models (LSSMs) with sparse system matrices, a setting relevant to many real-world dynamical systems where interactions are limited. The authors propose a maximum a posteriori (MAP)–EM algorithm that introduces Student-t sparsity-promoting priors on the system matrices and derives closed-form row-wise updates for all parameters. The algorithm is theoretically supported by a convergence guarantee in the EM sense (Theorem 3.3) and empirically evaluated on both synthetic and industrial datasets.", "questions": "1. **Thresholding fairness:** If the same post-hoc thresholding used in your method is applied to the MLE estimates, how do the resulting structure plots and MRE values change? Please report results for “MLE + threshold” using the identical threshold schedule.\n\n2. **Threshold sensitivity:** How sensitive are the topology recovery and MRE results to the chosen truncation threshold? How is this threshold selected in practice?\n\n3. **Non-invertible systems:** The examples presented involve invertible $A$ matrices. Can your proposed algorithm be applied to systems where $A$ is singular (non-invertible)? If yes, please provide an experiment to illustrate this case.\n\n4. **Hyperparameter robustness:** How sensitive are your results to the Student-t hyperparameters ($a_0, b_0$) and to the initialization of $\\Gamma$? Do you have recommended defaults or practical guidance for choosing these values?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. **Well-motivated problem.**\nLearning LSSMs with sparse system matrices is both practically and theoretically relevant, as sparsity naturally arises in many real systems with limited or localized interactions.\n\n2. **Clear algorithm design and convergence guarantee.**\nThe proposed MAP–EM algorithm is clearly described, with detailed derivations and an explicit convergence guarantee (Theorem 3.3), which adds credibility to the approach.", "weaknesses": "1. **Lack of a formal sparsity definition.**\n “Sparse system matrices” are introduced implicitly through element-wise priors rather than through a clear structural definition or quantitative metric for topology recovery. As a result, the structural claims rely mainly on visual inspection rather than measurable topology-level metrics.\n\n2. **Limited theoretical support for the “generalized permutation” claim.** The conclusion that the admissible similarity transformation $\\Phi$ reduces to a generalized permutation matrix is demonstrated only through an example, without a general proof or stated conditions under which this holds.\n\n3. **Potential unfairness in sparsity comparisons.** The authors apply post-hoc thresholding (“parameters below a predefined threshold are truncated to zero”) to their own results to emphasize sparsity, but it is unclear whether the same was done for the MLE baseline. Given that in Figures 1 and 2 the MLE method already yields a structure and MRE very close to the proposed approach, applying the same thresholding to MLE might lead to similarly sparse and accurate results, thus weakening the claimed structural advantage.", "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)