An Algorithm for Learning Switched Linear Dynamics from Data
Authors: Guillaume O Berger
Monal Narasimhamurthy
Kandai Watanabe
Morteza Lahijanian
Sriram Sankaranarayanan
Venue: NeurIPS 2022 Accept (2022)
Keywords: Control Systems
System Identification
Machine Learning
Hybrid Systems
Dynamical Systems
Ellipsoidal Algorithm
Separation Oracles

Abstract: We present an algorithm for learning switched linear dynamical systems in discrete time from noisy observations of the system's full state or output. Switched linear systems use multiple linear dynamical modes to fit the data within some desired tolerance. They arise quite naturally in  applications to  robotics and cyber-physical systems.  Learning  switched systems  from  data is a NP-hard problem that is nearly identical  to the $k$-linear regression problem of fitting $k > 1$ linear models to the data.  A direct mixed-integer linear programming (MILP) approach  yields time  complexity that is  exponential in the number of data points. In this paper, we modify the problem formulation  to yield an algorithm that is linear in the size of the data while remaining exponential in the number of state variables and the desired number of modes. To do so, we combine classic ideas from the ellipsoidal method for solving convex optimization problems, and well-known oracle separation results in non-smooth optimization.  We demonstrate our approach on a set of microbenchmarks and a few interesting real-world problems. Our evaluation  suggests that the benefits of this algorithm can be made practical even against highly optimized off-the-shelf MILP solvers.

An Algorithm for Learning Switched Linear Dynamics from Data
Guillaume Berger* Monal Narasimhamurthy" Kandai Watanabe Morteza Lahijanian Sriram Sankaranarayanan University of Colorado Boulder, Boulder, CO, USA f firstname.. lastnameccolorado edu
Abstract
We present an algorithm for learning switched linear dynamical systems in dis- crete time from noisy observations the system's ' full state or output. Switched linear systems use multiple linear dynamical modes to fit the data within some desired tolerance. They arise quite naturally in applications to robotics and cyber- physical systems. Learning switched systems from data is NP-hard problem that is nearly identical to the k-linear regression problem of fitting k >> 1 linear models the data. A direct mixed-integer linear programming (MILP) approach yields time complexity that is exponential in the number of data points. In this paper, we modify the problem formulation yield an algorithm that is linear in the size of the data while remaining exponential in the number of state variables and the desired number of modes. To do so, combine classic ideas from the ellipsoidal method for solving convex optimization problems, and well-known oracle separa- tion results in non-smooth optimization. We demonstrate our approach on of microbenchmarks and a few interesting real-world problems. Our evaluation sug- gests that the benefits of this algorithm can be made practical even against highly optimized off off-the-shelf MIIP solvers.
1 Introduction
Switched linear systems model dynamical systems that arise in diverse areas including natural sci- ences (biological models) [4, 8, 15], robotics [2], and cyber-physical systems [1]. They are charac- terized by a finite set of modes, wherein each mode features different set of governing equations for the future states in terms the current state of the system. In this paper, we study algorithms for inferring the dynamics of switched linear system from data that consists of full-state observations (this can subsequently extended to output observations using ARX models; see Section for details). Specifically, given set of data points {(x, x)) of current and next state observations, and number of modes m, we wish to find matrices A Am such that each data point (x, x() is explained by some matrix Ay: (Vi) (3j) ||x| - A,x;I| el|x;|| + т, for given eитоr tolerances T > 0. However, the switching signal which governs the current mode associated with each data point is latent. This renders the problem NP- and the best known approaches have exponential complexity in terms of the number of data points and the desired number of modes [22]. The ex- ponential complexity in the number of data points makes this problem especially challenging since many datasets of interest have thousands of data points but with tens of modes and state variables.
The key insight in this paper is that by reformulating this problem to incorporate "gap" in the tolerance, we can significantly improve the complexity to be linear in the number of data points N and exponential in the number of modes and the dimension of the state space. Rather than providing YES/No answer for single tolerance bound e, our approach takes as inputs two levels of tolerance a ( It either finds matrices that satisfy the "upper" tolerance bound c2 0 or
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

returns No if no set of m matrices can fit the data with the "lower" tolerance 61 (< E2). However, the underlying data could have been modeled with some tolerance that lies in the gap (f1. our algorithm can either succeed or return No with no guarantees provided. Our technical approach exploits the gap in our formulation to argue that any set of solutions must have a lower bound on its volume. Furthermore, by exploring the solution space in carefully defined manner, we guarantee that each step of our approach will shrink the volume of the remaining solution space by a constant factor. This is achieved by using ideas that are similar to separation oracles in combinatorial optimization [20], or the so-called "cutting plane" argument [10]. However, our approach applies these inside tree-based branch-and-bound algorithm, wherein we use the volume contraction to prove bounds on the depth of the underlying tree. Thus, we achieve a bound on the running time that is exponential in the state dimension and the number of modes but linear in the data size.
We demonstrate our ideas using a prototype implementation that is compared against a standard mixed integer linear programming (MILP) formulation of the problem solved using a state-of-the art solver (Gurobi) [17]. Our experimental results confirm that the theoretical insights also apply in practice to yield an algorithm that is often orders of magnitude faster, especially as the amount of data increases. We also demonstrate our approach using two interesting applications that include modeling human writing on a tablet in order to predict what alphabet is being written from a small number of samples, and deriving models of mechanical systems with contact forces from data.
1.1 Related Work
The problem of switched and hybrid system identification has been widely studied using vari- ety of approaches going back to the early 1980s; see, e.g.. [32]. We discuss a few representative approaches, and refer the reader to the monograph by Lauer al [22] for further details.
The identification approach may be exacr (seeking a global solution that minimizes the error be- tween the data and the model prediction), or approximate (wherein the optimization problem is solved approximately, or assumptions about the nature the switching signal are used to simplify the problem). Our approach here is exact but the problem itself is reformulated with gap. As for the nature of the switching mechanism, many approaches (including ours) focus on identifying the dynamics of individual modes while assuming that the switching signal is exogenous; on the other hand, other approaches, such as piecewise affine system identification [5, 28], hybrid automata learning [29], or linear complementary systems learning [19], attempt to identify the dynamics of each mode along with rules for transitioning between modes.
Vidal et al [33] present an exact approach for finding matrices that fit the data in the absence of noise by posing the problem as one of finding zeros of multivariate polynomial. This subsequently refined by Ozay et al [26] to handle the case of noisy data, by using sum-of-squares relaxations to obtain a semidefinite optimization problem. Assumptions on the nature of the switching signal can drastically simplify the problem at hand. For instance, Ozay [25] shows that the problem can be solved in polynomial time using a dynamic programming algorithm if the number of mode switches in trajectory is bounded. Other approximate approaches involve greedy algorithms [5] and block- coordinate descent (similar to k-mean regression) algorithms [21, 31j. However, the performance of these algorithms can vary depending on the nature of the data and guarantees available.
Machine learning techniques have also proved useful in hybrid system identification. For instance, standard machine learning approaches (such neural networks) can be used to infer hybrid system models [16, 24]. However, the resulting number of modes can be exponential in the network size. Moreover, switched systems often involve discontinuous switching which cannot be modeled ad- equately using standard activation functions such as sigmoids or ReLU. Ly and Lipson present an approach for learning hybrid automata from data using evolutionary techniques [23]. Their approach combines symbolic regression for learning the dynamics of the individual modes with techniques for guessing the latent modes and inferring the conditions for switching from one mode another.
Also relevant to this work, Dempster et al [12] study the problem of inferring hidden Markov models with linear systems, called Jump Markov Models (JMMs), using modification of the Expectation- Maximization (EM) algorithm. However, the E-step, which infers the assignments of modes, can be quite expensive for JMMs. Subsequently, Gharamani and Hinton [14] propose a technique wherein the difficulties in the E-step are resolved using variational inference, whereas Blake et al [7] propose sampling- based scheme. A key difference between our approach and JMMs is that we focus on

---

the algorithmic efficiency for learning the dynamics of the modes. Although this requires us to infer the sequence of modes for the data, we do not learn the Markov model that generates this sequence. However, our example on handwriting recognition presented in Section 4 illustrates how off-the- shelf approaches for learning a finite-state model can combined with our approach in fruitful manner. As for our application involving mechanical systems with contact forces (Section 4), a recent work by Jin et al [19] propose a convex relaxation approach to learn linear complementary systems. These systems provide compact representations for class of piecewise linear systems, including mechanical systems with contact forces [3, 27]. Our implementation does not use this compact representation, but is nevertheless able to learn switched linear models for these systems. Finally, our approach side-steps the problem of estimating the number of modes which is often an is- sue. This problem can be addressed by adding penalty term accounting for the "model complexity" to the overall objective function to be minimized [23] or by using more sophisticated non- non-paraettric Bayesian approaches wherein the prior distribution specifies unbounded number of modes and the varying number of model parameters varying as the number of modes increases [13].
2 Problem Formulation
We will use bold, lower case symbols (x, y, w) to denote vectors, and upper case letters (A, B, C) to denote matrices, unless otherwise mentioned. Let denote the set of real numbers, and N the set of natural numbers. For
matrix with all entries equal to 1. We write A & B to denote that each entry of A is less than or equal to the corresponding entry of B (same notation will be used for vectors).
N. let [m]] = (1, (1...., m}. Let E,, denote the
linear dynamical system in discrete time defined by its state x(I) E R" and output y() E IRP, at each time / E N. The state evolves using the dynamical rule: x((+1) = Ax((), wherein A E Rax=. The output is obtained as (t) = Cx(t), wherein C' E R2xx.
Definition (Switched Linear Dynamical System). A switched linear dynamical system with state x N + R" and output y : N - IR" with m 2 1 modes is defined by m matrices each time N, there exists a mode (t) E [m] s.t. x(t + 1) = ((t)) The switching signal N + [m] chooses mode for each time such that the continuous state x(t) evolves according to this mode. The output is given by y(t = Cx(t) according to a given matrix C' E IRP"".
Am- For
Remark 1. Switched affine systems, with x(t + Ayx(t) +b, for bj R", are handled by augmenting x with an extra component that is always 1, thus treating by as a new column of Aj-
The goal of this paper is to learn a switched linear dynamical system model from observation data. The simplest type of data consists of full state observations involving pairs of states of the form {(x,, x wherein each x; = x(4+) and x = x(4; + 1) are states observed at two successive time instants and t; + 1. In particular, the switching signal is not observed. We assume that the number of modes m is known.'
Remark 2. The assumption of full state observation is often impractical. However, our tech- niques extend to output observations by using fixed-length history of output observations y(t y(, and finding a switched auto-regressive (SARX) model [35] of the form:
y(t+ E Bo k), wherein (tt)= [m], B3 E RP*P.
Therefore, in the following, will focus our description on the full-state data. Definition (Fitting With Error Bound and Tolerance c). A set of matrices A Am "fits" the data D with error bound T > 0 and tolerance iff
(V (x,, x)) D) (3j E [ml) ||x] - A,xiloc S Ee|xx]loo0 +T.
(1)
Our approach can be modified handle T.2 or Ta norm instead of Too (see Appendix A).
The switched linear system identification (SLS-ID) problem inputs the number of modes m 2 2, data D {(x,, x!) Y. error bound 7 > 0, tolerance > 0 and a magnitude bound y. It outputs set of matrices Aq...., An that fits with error bound r and tolerance and satisfies
-, < A, < En, Vj E [m]] i.e, each entry of A₁ lies in [-7,7],
(2)
'Using repeated doubling followed by binary search, we can find the smallest number of modes m" for which switched linear system fits the data by running our procedure O(log(m") )) times.
3

---

or outputs INFEASIBLE if no such matrices exist. The SLS-ID problem is known to be NP-hard [22, Theorem 5.1]. We formulate below a mixed-integer linear program (MILP) for the SLS-ID problem.
We choose mn? decision variables which form the unknown entries of the matrices AI...., Am- We also use binary indicator variables Wi,1,--.;
{0, }} wherein each data item in D. The indicator variable Wa denotes that the 7th data item in D is fitted by the matrix Aj. Each data item must be fitted by one matrix, i.e,, Σj=1 Wi,j = 1, V Vic E [M] We add
,,...,, N ranges over
constraint (2) to set limits on the magnitude of the entries of the matrices Aj. Finally, we add the constraint that each data item is fitted by the corresponding matrix:
|| A,x: S <|/>/|>>> + T + (1 M, A (x,, x E D, Vj E [m]
(3)
Here, M > 0 is a constant that is larger than >/|x:||| || + ||x| oc for all (x, (x) E D. implying that (3) trivially satisfied if w.j 0 (so-called "big-M' .M" trick in integer linear programming).
The set of constraints above along with W1, E {o. }} form the constraints of MILP. The objective can be set in many ways: for instance, minimize the sum of the absolute value of the entries of each matrix, as a regularization term.
Lemma 1. A sel of matrices AI...., Am is feasible for the MILP defined above (for some valuation of the binary variables wi,)) )) iff it fits the data D with error bound T and tolerance e.
The branch-and-a and-bound algorithm for solving MILPs will solve m' N linear programs in the worst- case, noting that for each E [N] exactly one binary variable in the set (wi,1,-... Wi.m) ) can be set to 1. This is exponential in the data size N (= DD). This bound is seemingly independent of the dimension of the underlying state-space. However, if N mn, we are estimating more unknowns than the available data. Therefore, in practice, it generally holds that N > mn.
2.1 Identification with a Tolerance Gap
We now formulate a "relaxed version" of the SLSS-D problem with a "tolerance gap", that will be called the SLS-ID-GAP problem.
Inputs: Number of modes m 2 2, data D: {(x,, )) absolute error bound T > 0, two relative error tolerances 62 > G1 2 O and a magnitude bound 7. Let (gap = <2 - 61. For technical reasons
that will be explained later, we assume that Cgap-
Output: FEASIBLE along with set of matrices A A. satisfying the bounds constraints (2) that fits D with error bound and tolerance e2, or, INFEASIBLE if no set of m matrices with magni- tude bound 7 Égap can fit the data with tolerance e and error bound
For any given data set D and tolerance 62 > O, we define émin (D, e2) E [0, oo] as the least tolerance for which D can be fitted with fixed error bound and magnitude bound Y (e2 <).² Thus, the algorithm for the SLS-ID -D--IP problem guarantees that (a) if emin (D, E2) S e then the algorithm will return FEASIBLE with matrices that fit the data with tolerance 62 and magnitude bound y (b) if Ciiu (D, (2) > (2, the algorithm will return INFEASIBLE; (c) if Ciiu (m, 72) E (C1, C2], the algorithm may return either answer without any guarantees. The main result of the paper is follows:
Theorem 1. There exists an algorithm for solving the SLS-ID-GAR GAP problem with complexity O (m Cmne [ogins//<w)) N poly(m, n)), wherein C is constant factor and poly (0,, n) is polynomial function of m and n.
The significance of this approach is that it is linear in the size of the data set N, although exponential the number of modes m and the dimension of the state space n.
3 Tree-Based Exploration Algorithm
We now present, in stages, our algorithm for finding set of matrices A Am that fits the data. First, we present the algorithm as tree-based exploration approach, wherein each node of the tree
Cmin (D, E2) well defined; for instance, it suffices to solve MILP in the previous section with as objective function to minimize, and replacing 7 with 7 - (E2 -e) in constraint (2). Note that emin (D,ez) 62) oc iff there is x; = 0 with ||x|||>> > T (indeed, this is the only situation in which no set of matrices A1,... A. can tit the data for any tolerance 0).

---

Algorithm 1: Overall algorithm for switched linear system identification
Data: m. D, 61 62, 7 (see Section 2.1).
Result:
YES with matrices A1...., Am that fit the data with tolerance c2, or No. Initialize tree T with root node (see text for details)
2 while there exist unexplored leaf nodes in I do unexplored leaf node in T 4 Mark y as explored
resull Expand using Algorithm 2 if result = (FEASIBLE, (A1......Ae ...A)) then return (YES, A1........ Am) /* Solution discovered
8 return No /* No nodes remain to be explored
Algorithm 2: Algorithm to expand a non-terminal leaf node in the tree.
Data: Leaf node with polyhedra P1. PI...., P.. data set U and assignment map μ.
Result: New leaf nodes, or matrices Q1,-.. Qm that fit all the data.
Choose feasible solutions (matrices) Q1,..... Qm s.t. for all j E [m], Qj E (Section 3.2 will
specify how to choose Q1...., Qm) .mm)
Find (x,xx) U s.t. for all j E [m], ||x Qixillo > e2 <2||x:|100 +
3 if' not found then
L return FEASIBLE, (Q Qm) /* Qu.... Qm fit all the data
else
Ü -+D U. (Xx,,() for j E [m] do
/+ Remove (x,, x() from unassigned data
/* Constrain polyhedron Pj s.t. matrix Aj "fits" (x,xx)
P3 + P, {Aj: : || A,xilloo I <2|||x]loo +r}
if P, contains Lind ball of radius Cgap (Lemma explains why) then Create new child node V/ of
10
11
Associate vj with polyhedra PI...., P...,,,P,... Py. Pm
12
ūJ + U ((x, x)) + j}
/* Associate (X,, x() with mode j.
Associate node V/ with set Ú and map úú
13
represents a set of m convex polyhedra. Subsequently, we will show how careful choice of the tree exploration strategy will guarantee bound on the maximum length of each branch of the tree. This, in turn, yields the desired algorithm and its complexity guarantee. The proofs of all the lemmas are outlined: detailed proofs are provided in Appendix B (cf. supplementary material).
The central data structure maintained by our algorithm a tree Each node of the tree has the following associated information:
1. Convex polyhedra Pi...., Pm. represented as systems of linear inequalities. Each P, C Raxx describes a set of possible solutions for matrix Aj.
subset C D containing data points that have not been "assigned" a mode yet.
3. A assignment map μ (D U' + [m] mapping each data point in D \ U to a mode in [m].
Root Node: The root the tree initialized by m convex polyhedra Pi...., .a....,Pm. wherein P. encodes the bounds each entry of the matrix Ay: -nE, A, S E,, The set U at the root is the full data set and the associated map y is the "empty map" since its domain is empty.
3.1 Expanding Node
Starting from the root node, our algorithm iteratively expands the tree until no unexplored leaf remains suitable set of matrices is found (see Algorithm 1).
5

---

Expansion: Each expansion step involves choosing an unexplored leaf v of the tree and carrying out the steps outlined in Algorithm 2. Namely, expanding v begins with choosing feasible solutions Q1..... Qm from polyhedra Pı, Pm. respectively (line 1). The choice of these solutions will be described in detail in Section 3.2. Next, we sean through all the "unassigned" data U to find a data point (x:, x.) E U that cannot be fitted by Q1..... Qm (line 2). If no such data point can be found then we have, in fact, fitted all the data using Qi,..., Qm Therefore, we can terminate (line 4). If data point (x,, x() E U is found in the previous step, then we create potential child nodes Vi,..., V' for the current node. Each potential child vj will have associated polyhedra PI...., PJ-1. P.. Ps+1..... Pm. wherein, for j' # j. Py remains the same as in node v, and P3 is defined in order force the matrix A, to fit the data point (x,x() that could not be fitted by the previous candidate matrices (line 8). Node Vy will have associated data set Ú U \ {xx,, x()) and mode assignment map fy : 1 U {(x,,x] + j}.² Finally, the potential child node v/ is actually added to the tree only if P, satisfies some lower bound on its size (line 9).
Thus, the process of expanding node of the tree results either in termination with a feasible solution QI,.... Qm that fits the entire data, or in the possible addition of at most m new leaf nodes to the tree. We will now establish important properties of this algorithm. Lemma 2. Let Qi,..., m be the chosen candidates while expanding some node v. For each child Vy of u with associated polyhedra PI....: .....:P-- PJ-1, P3: PJ+1.-.., Pmn we have that Q3 & Pj.
Qm
The proof (see Appendix B) shows that the constraints added to Pj in line 8 are not satisfied by Qj-
Recall that the induced (operator) Loe norm of a matrix: ||A]]in |ind - max st AxllAxx]ss Given A E R=*= and e 2 O, let (A, denote the ball of radius (W.r.. I1-11.ed) centered at A. We prove that the tree-based exploration algorithm is complete for the SLS-ID-GAP problem: if there is a set of matrices ....... A,m
fitting the data with tolerance 6 and satisfying the reduced magnitude bound 7 Egap. then there is an unexplored leaf whose associated polyhedra contain these matrices. Lemma 3. At each iteration of the algorithm, there exists an unexplored leaf v in the tree with associated polyhedra Pl...., , Pm such that for all j E Bind (Aj)) S P.
The proof (see Appendix B) is by induction on the number of iterations of the algorithm. At the beginning, we show that the root node satisfies the property of the lemma. Then, we prove that, at each iteration, if a leaf satisfying the property is picked to be expanded, then at least one child of the leaf is created and "inherits" the property. The requiremen that y > Egap is key for this lemma.
From the definition of the child nodes (line 9), the convex sets associated to each leaf of the tree satisfy lower bound their volume.
Lemma 4. At each iteration of the algorithm, and for each leaf v in the tree, with associated polyhedra Pi,.... Pan. it holds that for all j E [m], vol(P;) ); 2 (2 )²
(m!))
Proof. Let be a leaf node as in Lemma 3 and let j E [m]. By Lemma 3, Pj must contain (A), (A Gaap) The ball Bind (Aj,(gap) is the product of n unit Ly-norm balls (one for each row) each scaled with factor Egap- The volume of unit norm ball in n dimensions is given by 2". n!. Combining these observations, we have that (o)( > (2epp)")
m!
We have not proven a bound on the size of the tree explored by our algorithm so far. In the next subsection, describe way of selecting the candidates Qi,..., Q.. (line 1) such that the volume of P3 will be smaller than some fraction a < 1 times the volume of Pj. This, combined with the lower bound on the volume of these sets (Lemma 4), will provide bound on the length of each branch (i.e.. the depth) of the tree.
3.2 Cutting Plane Argument
We prove an effective bound on the depth of the tree using the so-called cutting-ppane argument from non-smooth optimization [10]. First, we refine line I of our algorithm wherein we choose candidates
In other words, the newly assigned data point is assigned to mode j in node v/j
"Note that | = max([|AA.] ....).,.................. ||A.. 11) (maximum among the L1 norm of each row) [18].

---

Q2 Pj for each j E [m] Specifically, we will choose Q, as the center of the maximum volume inscribed ellipsoid (MVE) of P,. The MVE center of polyhedron can be computed efficiently by using semi demi-definite programming (SDP) [6, Proposition 4.9.1]. Now consider the child node Vj of v, to which we associate the polyhedron P3 & Pj. We show that the volume reduces by at least factor = (1 --)=)= 1.
Lemma 5, vol(P;) S (1 ne )vol(P)).
The above follows from the classical cutting plane argument [30] applied on the polyhedron Ps. which is obtained by cutting Pj with hyperplane excluding the MVE center (see Appendix B).
Let Vminn (26pe)" )n denote the bound proved in Lemma (n!)"
Lemma 6. The depth of the tree is (mm" 1og(n7/<ap))
The proof (see Appendix B) observes that each time we expand node in the tree, polyhedron has its volume reduced by a factor a as stated in Lemma 5. Also, when the overall volume becomes smaller than Viin in Lemma 4, the node is pruned.
This places bound on the depth of the tree, as stated in the lemma. In fact, the n term is reduced to n² using the observation that each polyhedron Pj for each node v is the Cartesian product of polyhedra PjA.. for k E [r2], involving the variables from the 1thh row of the unknown matrix Aj. This "fine-grie mraingea analysis provided in Appendix C.
The overall size of the tree cannot exceed m O(mn log(ny/<wo)) nodes, wherein the complexity of expanding each node is linear in the size N of the data set and involves solving m SDPs each with ni variables and O(mr²) constraints.
Theorem 2.
4 Experimental Evaluation
In this section, we will describe an evaluation of our approach meant to answer two key questions: (a) Do the theoretical guarantees translate into superior empirical performance when compared to highly optimized MILP solvers? (Namely, we compare our approach against the MILP solver Gurobi [17] over a set "microbenchmarks" of varying dimensions, number modes and data sizes.); (b) Does the approach yield interesting results on real-life datasets? (Namely, we illustrate our approach on datasets from handwritten alphabets and mechanical systems with contact forces.)
Implementation: We implemented the proposed approach in Python 3.8, using Gurobi [17] to encode and solve linear programs. We assume that the number of modes m (or an upper bound on it) is given. Also, our implementation fixes e 0 and 62 Egap > 0. Fixing 6 = O implies that our algorithm either finds matrices that fit the data with tolerance e2 or concludes that no matrices exist that fit the data with zero relative етror tolerance. Therefore, throughout this section, we will report the value of e. Our implementation coincides with the description Section 3 with one important modification:
compute the Chebyshev center (center of the largest inscribed ball) instead of the center the maximum volume inscribed ellipsoid (MVE). The Chebyshev center can be computed very efficiently and reliably using Linear Programming [9], and provides a good approximation of the MVE center. Although the theoretical guarantees on the termination of the process using the Chebyshev center weaker than those with the MVE center (Lemma 5), the center in this context is a widely-used heuristic [10, $4.4). The containment check described in line 9 of Algorithm 2 is implemented by comparing the Chebyshev radius against Egap.
of Chebyshev
We compare the proposed approach against the MILP ap- proach on suite of synthetic microbenchmarks. Each such benchmark consists of randomly chosen n X Hurwitz matrices, where m. n are varied systematically. A total of N data points were generated for each experiment from = N/T trajectories, each of length T 10 time steps starting from a random initial state in [-1, 1jn. An additive noise sampled uniformly at random in the range [-0.05, 0.05] was added to each state. Figure we illustrate how the two approaches scale in terms of computation time with respect to the number of data points (N), number of states (1), and number of modes (m). All timings are averaged over 10 separate runs to account for the variability in computation times. Further details microbenchmark generation are provided in Appendix D.
Evaluation on Mierobenchmarks:
7

---

150
15
0.5
+ T5 + MILP
1.5
TS
-
MILP
MILP
0.4
100
time (s)
1.0
0.3
50
0.2
0.5
0.1
0
0.0
102
10"
2
6
8
data points, N
states,
of modes, m
(a) Performance across 10 random (b) Performance across 40 random (c) Performance across 40 random microbenchmarks with m 3, and varying values of N.
4, microbenchnarks with 100. and varying values n. 200, and varying values m.
4. microbenchmarks with
4.
of microbenchmarks. All timings are reported in seconds on a Linux server running Ubuntu 22.04 OS with 24 cores and 64 GB RAM. Each point in the plot represents the average time taken by the algorithm across 100 experiments (10 runs for each the 10 microbenchmarks). The error bars represent the minimum and maximum values of time taken across experiments. (a) The proposed approach (TS) scales better than the MILP approach as the number of data points increases and has smaller variance. Both approaches scale similarly with the dimension (b) and the number of modes (c).
Figure 1: Performance of MILP proposed approach (TS)
Table I: Performance of k-Linear Regression vs. our approach on of microbenchmarks. We N data points for training both approaches and 50 data points constitute the held-out test dataset. Each row of the table reports average/min/max over 5 runs with 7 0.05 and 1 for fixed train-test data set. All experiments were carried out on a Linux server running Ubuntu 22.04 OS with 24 cores and 64 GB RAM.
PROPOSED APPROACH
k-LINEAR REGRESSION
Error avg
Time
Error max
Time
N
max
min
avg (8)
avg
min
avg (s)
Benchmark #1
20 0.17 1.16 0.01
0.04
0.34
3.21
0.01
0.03
(n 2, m = 4)
100 0.05
0.33 0.00
8.54
0.23
2.13
0.01
0.02
500 0.05
0.37 0.00
15.61
0.19
2.15
0.00
0.02
Benchmark #2
20 0.14 1.10 0.01
2.62
0.29
2.13
0.02
0.08
(n 4, m = 2)
100 0.07 0.46 0.01
1.37
0.21
2.31
0.02
0.02
500 0.06 0.26 0.02
1.56
0.27
2.47
0.02
0.03
Benchmark #3
20 0.56
2.72 0.03
1.70
0.95
5.58
0.02
0.16
(n = 4)
100 0.23 1.65 0.03
506.29
0.48
3.74
0.01
0.41
500 0.18 1.24 0.02
1184.16
0.41
2.46
0.02
0.11
The comparison clearly shows that the MILP solver's computation time increases rapidly with the number of data points N. In contrast, the running time of our approach depends linearly on N. Both approaches scale similarly in terms of number of modes (m) and dimension (n))
We compare the accuracy of the proposed approach against k- Lineea Regression (k-LR) [21] in terms of one-step prediction errors: minf" |x| - -A, loo measured over held out test dataset. Table compares the prediction errors for a set of 3 microbenchmarks. The proposed approach has smaller error bounds compared to kLR but takes more time, as expected.
Handwriting Recognition: We now evaluate our approach on a dataset from human handwriting on a tablet: our goal is identify various modes with dynamics that describe how letters are traced out on the tablet. We generate our own handwriting dataset by having "a", "b", "c" and "d" using their fingers on the mouse pad of their laptop. We collect the (a, y) locations of the handwritten letters over time. We generate our final dataset by interpolation in order to ensure that the samples are roughly equidistant from each other. This can be done in an online fashion, if need be. It ensures that our model here is not capturing artifacts of the pressure sensor in
author trace out the letters
8

---

1: 0.92
A
start
b
b
b
b
b
b
pred:
C
C
C
C
C
0.08
1: 0.01
pred:
d
d
(A
pred:
Figure 2: Left:
Handwriting samples with identified modes shown in different colors, Middle: JMM model for letter b and Right: predicting letters from partial observations using the JMM.
Figure 3: Left: Cartpole with soft walls. Right: Acrobot with soft joint limits.
the mouse pad or how the person tracing the alphabet may speed-up or slow-down various points in time. We also scale the (a, y) values so that the letters fit within a square of side length 1. Further details are provided in Appendix We collect 10 samples for each letter (Figure 2-1eft).
We apply our approach to learn m = 3 matrices that fit the data with c = 0.05 and T = 0.1. For each sample, obtain corresponding switching sequence from this process. We apply well- known grammatical inference approach to learn a Jump Markov Model (JMM), wherein the various dynamical modes are represented on the edges of the automaton [11, 34]. The estimated mode for each data point of each letter is represented in Figure 2-1eft. Each color represents a different mode. We note that the modes naturally corresponds to different parts of the letters.
Figure 2- 2-midlle shows an example of JMM for the letter "b". Each transition in the JMM represents mode and probability of taking the edge. At state go. mode 1 is observed for 92 % of the time, and at G1, mode is observed for 93 % of the time. We then try to predict a class for an alphabet as it is being written, based on the estimated JMMs (see Figure 2-right). Each column represents 20 %, 40 %, 60 %, 80 %, 95 % portion of a letter. We use the JMMs and Bayesian inference to compute the probability of a letter given partial sequence of states. The letters "a", "b" and "d" are correctly predicted but the letter "c" is classified as an "a": a visual inspection of Figure 2-right explains why this is the case.
In conclusion, we show that our switched linear system identification algorithm works on the noisy handwriting data and can be used to identify JMM from sequences of identified modes which in turn can be useful for data classification and prediction.
Acrobot: We evaluate our approach on an acrobot benchmark with soft joint limits [3]. We sam- pled N = 30 trajectories with time step = 0.03 seconds for T = seconds. We then identified the dynamics at each mode using the proposed approach (with m 3, 0.01, and e 0.01). The average one-step prediction error" of the proposed approach on a held-out test dataset consisting of 5 trajectories is 0.005 (min: 0.0, max: 0.18). The average training time is 7.56 seconds. In compar- ison, the average one-step prediction error of k--ineer Regression( [21] on the same dataset is 1.80 (min: 0.O, max: 8.68). The average training time of &LR is 0.18 seconds. Figure 4 shows the predictions of the proposed approach on sample test trajectory. Despite the underlying system having infinite number of modes (as the time-sampled system of continuous-time hybrid linear system), our system identification technique is able to identify three main linear modes that explain most of the data both in the training and test datasets.
Onn-step prediction min.- ||x] A,x,

---

2.0
1.5
Ref
2.0
15
1.5
2.0
0.5
0.5
Ref
100
0.0
timestep
75
timestop
10-
T5
Ref
Ref
-5
50
timestep
timestep
Figure 4: Data from the Acrobot system plotted against those of the identified model.
Cartpole: We evaluate our approach on a cartpole system with soft walls [3]. We sampled 30 trajectories with time step dt = 0.05 seconds for T = 5 seconds. We then identified the dynamics at each mode using the proposed approach (with m 3, r = 0.1, and c = 0.1). The average one-step prediction error of the proposed approach on held-out test dataset consisting of 5 trajectories is 0.016 (min: 0.00, max: 0.15). The average training time 9.35s. The average prediction error of k-Linear (egression(kLLR) [21] is 0.04 (min: 0.0, max: 0.15) and the average time taken by kLR is 0.16 seconds. Additional plots are provided in Appendix E. Similarly to the acrobot example, the learned system is able to explain most of the data (both in the training and test datasets) with a few (three) linear modes, despite the fact that the actual system has an infinite number modes (as the sime-sampled system of continuous-ting hybrid linear system).
5 Conclusions: Limitations, Future Work and Broader Impacts
In summary, we have presented an approach that improves the time complexity of switched lin- system identification using novel problem formulation with "gap" and using ideas from the ellipsoidal method in combinatorial optimization. The resulting complexity is linear in the num- ber of data points. Limitations include (a) Our approach uses full state information and extends to output observations using SARX models. However, the number of modes in the SARX model can be exponentially larger when compared to the state-space model. (b) Our approach finds the matrices fitting dynamics but does not infer the "generator" for the mode switches, although we present an example using Jump Markov Models. However. the problem of learning the matrices for various modes in conjunction with the "generator" may not always lend itself to a clean separation wherein the modes are first inferred using our approach and generative model such as Markov chain is then inferred. For instance, our approach may find matrices that fit the data well but lead to poorly fitting system models. We will study the use of backtracking in our approach to generate multiple "dissimilar" solutions; (c) Finally, our implementation performs the tree exploration with- out leveraging other information about the system or using heuristics that have proved powerful in branch-and bound algorithms. Nevertheless, it provides proof of concept that it can be competitive against opttiieo MILP solvers. We will expand on these ideas in our future work.
In terms of broader impacts, model identification is very important machine learning problem that has very important and beneficial applications in areas such as controls, autonomous systems and medical applications. Better modeling of physiological processes can lead to improved closed- loop medical devices. Improved ability to predict human movements can lead to safer human-robon robott interactions. However, our approach be used to improve surveillance systems that could be themselves be used in ways that are detrimental to the freedom individuals and societies.
We thank the anonymous reviewers for their detailed comments and sugges- tions. This research was funded in part by the Belgian-American Education Foundation (BAEF) and the US National Science Foundatton (NSF) under award numbers 1836900 and 1932189.
Acknowledgments:
10

---

References
[1] Rajcev Alur. Principles of cyber-physical systems. MIT Press, 2015.
[2] Aaron D. Ames. Human-inspired control of bipedal walking robots. IEEE Transactions on Automatic Control, 5965):1115-1130, 2014.
[3] Alp Aydinoglu, Victor M Preciado, and Michael Posa. Contact-aware controller design for complementarity systems. In 2020 JEEE International Conference on Robotics and Automa- tion (ICRA), pages 1525-1531. IEEE, 2020.
[4] E. Bartocci, F. Corradini, M. R. Di Berardini, E. Entcheva, S. A. Smolka, and R. Grosu. Model- ing and simulation of cardiac tissue using hybrid I/O automata. Theoretical Computer Science, pages 31 3149-3165, aug 2009.
[5] Alberto Bemporad, Andrea Garulli, Simone Paoletti, and Antonio Vicino. A bounded-error approach to piecewise affine system identification. IEEE Transactions on Automatic Control, 50(10):1567-1580, 2005.
[6] Aharon Ben-Tal Arkadi Nemirovski. Lectures on modern convex optimization: analysis, algorithms,, and engineering applications. SIAM, 2001.
[7] Andrew Blake, Ben North, and Michael Isard. Learning multi-class dynamics. In Advances in Neural Information Processing Systems, volume 11. MIT Press, 1998.
[8] Luca Bortolussi and Alberto Policriti. Hybrid systems and biology. In Marco Bernardo, Pier- paolo Degano, and Gianluigi Zavattaro, editors, Formal Methods for Computational Systems Biology, pages 424-448, Berlin, Heidelberg, 2008. Springer Berlin Heidelberg.
[9] Stephen Boyd and Lieven Vandenberghe. Convex optimization. Cambridge University Press, 2004.
[10] Stephen Boyd and Lieven Vandenberghe. Localization and cutting-plann methods. From Stan- ford EE 364b lecture notes, 2007.
[11] Colin De la Higuera. Granmatical inference: learning automata and grammars. Cambridge University Press, 2010.
[12] A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Sociery. Series B (Methodological), 39(1):1-38, 1977.
[13] Emily Fox, Erik Sudderth, Michael Jordan, and Alan Willsky, Nonparametric Bayesian learn- ing of switching linear dynamical systems. In D. Koller, D. Schuurmans, Y. Bengio, and I.. Bottou, editors, Advances in Neural Information Processing Systems, volume 21. Curran Associates, Inc., 2008
[14] Zoubin Ghahramani and Geoffrey E. Hinton. Variational learning for switching state-space models. Neural Computation, 12(4):831-864 2000.
[15] Ronojoy Ghosh, Keith Amonlirdviman, and Claire Tomlin. hybrid system model of planar cell polarity signaling in drosophila melanogaster wing epithelium. In IEEE Conference Decision and Control (CDC), pages 1588-1594. IEEE, 2002.
[16] Ian Goodfellow, Yoshua Bengio, and Aaron Courville, Deep learning. MIT Press, 2016.
http::/www. http///wwwddeeeeeeennnnnnnnooo// org.
[17] Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2022.
[18] Roger A Horn and Charles R Johnson. Matrix analysis. Cambridge university press, 2012.
[19] Wanxin Jin, Alp Aydinoglu, Mathew Halm, and Michael Posa. Learning linear complemen- tarity systems. Learning for Dynamics and Control Conference, pages 1137-1149, PMLR, 2022.

---

[20] Bernhard Korte and Jens Vygen. Combinatorial optimization: theory and algorithms. Springer, 2010.
[21] Fabien Lauer. Estimating the probability of success of simple algorithm for switched linear regression. Nonlinear Analysis: Hybrid Systems, 8:31 -47, 2013.
[22] Fabien Lauer and Gérard Bloch. Hybrid system identification: theory and algorithms for learn- ing switching models. Springer, 2019.
[23] Daniel L. Ly and Hod Lipson. Learning symbolic representations of hybrid dynamical systems. Journal of Machine Learning Research, 13(115):3585 3618, 2012.
[24] Daniele Masti and Alberto Bemporad. Learning nonlinear stace-space models using autoen- coders. Automatica, 129:109666, 2021.
[25] Necmiye Ozay. An exact and efficient algorithm for segmentation of ARX models. In Ameri- can Control Conference (ACC), pages 38 41, 2016.
[26] Necmiye Ozay, Constantino Lagoa, and Mario Sznaier. Robust identification of switched affine systems via moments- based convex optimization. Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference, pages 4686 4691, 1, 12 2009.
[27] Arvind U Raghunathan, Devesh K Jha, and Diego Romeres. PYROBOCOP: Python-based robotic control & optimization package for manipulation. In 2022 International Conference on Robotics and Automation (ICRA), pages 981-991. IEEE, 2022.
[28] Sadra Sadraddini and Calin Belta. Formal guarantees in data-driven model identification and control synthesis. In Proceedings of the 21st International Conference on Hybrid Systems: Computation and Control (part of CPS Week), pages 147-156, 2018.
[29] Miriam García Soto, Thomas A Henzinger, and Christian Schilling. Synthesis of hybrid au- tomata with affine dynamics from time-series data. In Proceedings of the 24th International Conference Hybrid Systems: Computation and Control, pages 1-11, 2021.
[30] S. P. Tarasov, L. G. Khachiyan, and I. Erlikh. The method of inscribed ellipsoids. In Dokl. Akad. Nauk SSSR, volume 37, pages 226-230, 1988.
[31] Sacid Tizpaz-Niari. Pavol Cemn, Bor- Yuh Evan Chang, and Ashutosh Trivedi. Differential performance debugging with discriminant regression trees. Proceedings of the AAAI Con- ference on Artificial Intelligence, volume 32, 2018.
[32] Jitendra K. Tugnait. Detection and estimation for abruptly changing systems. Automatica, (K(5):6076115, 1982.
[33] René Vidal, Stefano Soatto. Yi Ma, and Sankar Sastry. An algebraic geometric approach the identification of a class of linear hybrid systems. In 42nd IEEE International Conference on Decision and Control, volume pages 167-172 Vol. 2003.
[34] Kandai Watanabe, Nicholas Renninger, Sriram Sankuranarayanan, and Morteza Lahijanian. Probabilistic specification learning for planning with safety constraints. 2021 JEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pages 6558-6565 IEEE, 2021.
[35] Siep Weiland, A Lj Juloski, and B Vet. On the equivalence of switched affine models and switched ARX models. Proceedings of the 45th IEEE Conference Decision and Control, pages 2614-2618 IEEE, 2006.
12

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you descrihe the limitations of your work? [Yes] This was discussed partly in the
related work and partly in the concluding Section 5.
(c) Did you discuss any potential negative societal impacts of your work? [Yes] see Sec- tion 5.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes] The proofs are sketched in the main paper but provided in full details as part of the appendix.
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material URL)? [Yes]
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes]
(c) Did you report error bars (e.g., with respect to the random seed after running exper- iments multiple times)? [Yes] We averaged across multiple runs and include error bars.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [N/A] Our experiments each ran on single CPU on desktop. We did use GPUs, clusters or cloud computing resources.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) your work uses existing assets, did you cite the creators? [Yes] Our approach is built on top of Gurobi [17] optimization tool.
(b) Did you mention the license of the assets? [Yes]
(c) Did you include any new assets either in the supplemental material URL? [Yes]
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] We are using publicly available data that can be down- loaded by everyone fireely
(e) Did you discuss whether the data you are using/curating contains personally identifi- able information or offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
13


Reviews:
[{"review_id": "gvM1DIQP3l5", "reviewer": "Reviewer_h6Sd", "summary": "Summary: This paper studies the identification of switched linear dynamics where the underlying mode is hidden (but the total number of modes is known). The contribution is an approach whose computation complexity scales linearly in the number of samples, as opposed to the exponential scaling in the MILP approach. However, the complexity of the proposed approach still scales exponentially with the number of modes and the state dimension. The approach is a tree search algorithm, and the core innovation is a pruning-based technique on the volume. This limits the depth of the tree. The computation complexity of the approach is provided. \n", "questions": "Line 4 of algorithm 1: 1: how are the $Q_1,\\ldots,Q_m$ chosen. Arbitrarily from $P_j$? \n\nLine 6 of algorithm 1, I assume $(x_i,x_i’)$ is from line 2. However, what if Line 2 returns multiple pairs of $(x_i,x_i’)$? \n\nLine 232: should it be the mode $j = \\mu^*(x_i,x_i’)\\in [m]$? \n\nLine 224: where does the $2^{mN}$ compute from?  Based on my read of the algorithm, at each layer, each node splits into m nodes, and the maximum depth is N. So shouldn’t the complexity be $m^N$? \n\nTable 1, benchmark 2, tau = 0.00, for the TS run time, why the run time for N=10000 is smaller than N=1000. \n\nLine 236: a “stronger” assumption is introduced here - does the main result Theorem 1 relies on this assumption? If so, this stronger assumption should be introduced in or before Theorem 1. \n\nIn both Lemma 4 and Lemma 5: there is a notion of “iteration”, which I think means each round of expansions of the current leaf nodes. Also, in Line 257, there is a redefinition of a step of the algorithm. Going back to my earlier comment on pseudo-code, I think the paper needs a formal presentation of the full algorithm so that the “iteration” is formally defined and the details of the algorithm are finalized. \n", "limitations": "The authors have addressed the limitations and potential negative societal impact of their work", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Weakness: \n\nStill, exponential dependence on the number of nodes and state dimension is needed. The author needs to provide a discussion on why the exponential dependence is difficult to avoid. For example, I don’t really see why the exponential dependence on $n$ is needed. Based on the analysis, the exponential dependence comes from Lemma 6, but is there a way around it?  Also, the MILP algorithm only scales exponentially with $m,N$, not $n$, right? If so, this means the MILP algorithm could scale to high dimensions, while the proposed algorithm cannot. This should be acknowledged and discussed. \n\nComparison: the majority of the paper focuses on comparison with the MILP approach, which is fine. However, in the related work session, a few other approaches are mentioned, e.g. the Ozay et al one and the clustering approach. How does the proposed approach compare with these related work, both theoretically and empirically? \n\nPresentation: the only pseudo-code is Algorithm 1, which is a subroutine on how to expand a node in the search tree. I suggest adding the pseudocode for the entire algorithm. \n", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "dA4tuMR-RCu", "reviewer": "Reviewer_Vd1p", "summary": "The authors study learning switched linear dynamics from transition pairs of states. They exploit a gap in tolerance (over which performance is ill-defined) to guarantee that when a set of feasible matrices becomes too small, no solution is consistent with the constraints. This is used to prune the search tree of a mixed-integer program (MIP). The size of the tree is bounded using a cutting-plane argument when candidate solutions are picked as the centers of maximum volume ellipsoids (within constraint polytopes). The method is tested on synthetic data, author-curated handwriting data, and highway driving data.", "questions": "- How are volumes of constraint polytopes computed? Are they computed, estimated, or are nodes simply pruned once the tree reaches depth $K_\\min$? This is a minor concern but it should be addressed.\n- Why exactly were handwriting data preprocessed to keep points roughly equidistant?\n- Should Lemma 9 read $\\mathsf{vol}(\\widehat{P}_j) \\leq (1 - \\frac{1}{n})^n \\cdot \\mathsf{vol}(P_j)$? It seems like the multiplicative factor $(1 - \\frac{1}{n})$ should be included for each row of the matrix $A_j$. Additionally, this would sharpen the $n^3$ in Lemma 10 to $n^2$.", "limitations": "I believe the authors have adequately captured the limitations of their work and have thought about potential negative societal impact. I think the submission could use a comment on how approaches for switched nonlinear systems may need to be developed (what kinds of approximations could be necessary, localized analysis, etc.), though this is clearly out of the current scope of the work.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n- The paper is well-written, the proofs (in the appendix) are easy to follow and the proof summaries provide good baselines from which most of the details can be inferred.\n- To my knowledge, the idea is novel and can be extended to several use cases beyond the proposed setting, such as control systems, state-based hybrid systems (with mode transitions at guards), and general mixed integer convex programming.\n\n# Weaknesses\n- It is unclear how to specifically prune leaves, as this involves computing the volumes of complex constraint polytopes.\n- Some of the design decisions and experimental conclusions were questionable. Specifically:\n  - Modern semidefinite programming (SDP) solvers can handle fairly large-scale problems, I would be interested to see how solving for maximum volume ellipsoid centers (as opposed to Chebyshev centers) would affect the performance of the overall approach.\n  - The authors claim (in the checklist) that error bars are not informative, but the timing data for TS in Benchmark #2 is couterintuitive (time is not monotonically increasing); this makes me wonder if the behavior is due to a bug, third-party implementation details, randomness, etc., and standard deviations over more trials could help make this behavior more obvious. \n  - The decision to space data roughly equally (in Euclidean distance) for the handwriting data set isn't completely justified for me; sometimes real-world systems are affected by dynamics with very different transition behavior throughout the state space. I would be interested in experiments without this preprocessing step. \n  - While the predicted velocities for highway data seem good, the \"behaviors\" extracted from the trajectory modes (in Fig. 2 Right) are not clearly distinguished or obvious to me.\n\n# Typos/Other Issues\n- There is a period in the title\n- Line 88 and 324, Markov should be capitalized\n- The sentence on lines 101-102 is confusing\n- Line 137 should read \"and *an* absolute limit $\\gamma$\"\n- For Lemma 1 to be necessary and sufficient as worded, the matrices should explain the data away with the mentioned error bound and tolerance, but also the absolute bound $\\gamma$\n- Line 160, data should be lowercase\n- Page 6 was the first place I noticed nonstandard use of colons (to indicate setting a variable?)\n- I believe the cited proposition on lines 260-261 should be Proposition 3.7.1\n- Line 273, Cartesian should be capitalized\n- Line 336-337 should read \"which in *turn* can be useful\"\n- Line 348 should read \"turning for changing *lanes*\"\n- Line 495 should read \"*computing* convex sets\"\n- I believe the citation in Line 508 should be Section 8.4.2\n- Should Line 553 (also 592) say that the path from the root to the leaf is of length *at least* $mK$? This seems to agree with the use of $T$ later (which is not introduced)\n- Line 569 should read \"as in *the* proof of Lemma 7\"\n- The proof of Lemma 10 is repetitive after the proof of Lemma 7, consider modifying this to highlight only the differences\n- Line 617, \"linear line\" is redundant\n- Line 619, I believe $i$ should take values from 0 to $T_{\\mathrm{raw}} - 1$\n- Page 17, the figure is hard to read. Moreover, the weights on each edge are all listed as 1.0", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "YrxuQjuSXoh", "reviewer": "Reviewer_FbXU", "summary": "This work deals with the problem of learning a switched linear system model from state-trajectory data, where the underlying mode switches are not observed. The goal is to learn state transition matrices and an assignment of each transition in a given data set to the corresponding mode, i.e., linear transition function, such that a certain error criterion is fulfilled.\nSince the naive problem setting based on Mixed-Integer-Linear-Programming is computationally infeasible, the authors propose a relaxation of the problem. For the relaxed setting, a concrete algorithm is presented, its correctness is proven and its time complexity is investigated. Finally, the algorithm is applied to synthetic task, a handwriting prediction task and a traffic data set.", "questions": "## Questions\n* How exactly could the SARX model be used to deal with output measurements? Does the algorithm and its analysis need adaption to this case?\n* Page 5, line 203: \"then we create $m$ child nodes\" This seems incorrect, since in Algorithm 1 a child node $\\nu_j$ is only created if $\\hat{P}_j$ is feasible\n* Proof of Lemma 4, lines 531-532: Does such a leave node (with $\\mu=\\mu^\\ast$) always exists?\n* Proof of Lemma 5: According to the proof, the lower bound should be $\\epsilon_\\text{gap}^n2^{n^2}/(n!)^n$ How is the lower bound stated in the lemma then derived?\n* Experiments:\n  - Why the choice of $\\epsilon_1=0$? Does it make the problem easier?\n  - the MVE is replaced by Chebychev centers. What exactly is lost when using Chebychev centers, i.e., which guarantees cannot be given anymore?\n  - What is the MILP approach? No details are given\n* Page 7, line 303: Why $N=kT$? If the trajectories have $T$ time steps, then there are only $T-1$ transitions (and each data point correponds to a transition)\n\n## Additional detailed remarks\n* p1, l26: Why \"explain away\"? This is somewhat of a technical term in e.g. Bayesian networks\n* p1, l28: Reference for NP-hard would be good here\n* p1, l30: Reference for exponential complexity would be good here\n* p2, l37: Unclear, what is $\\epsilon_2$?\n* p2, l39: The rest of the paragraph is very unclear. E.g., it is unclear what is meant by \"tree\" here.\n* p2, 48: Which approach is used to generate the optimization problem solved by GUROBI? The latter is just a general-purpose MILP solver\n* p2, l52: What does it mean to \"predict the alphabet\"?\n* p2, l56: Is reference 28 the earliest one in this area?\n* p2, l61: Wording unclear, \"over a formulation\"\n* p2, l74: Unclear how minimizing the number of mode switches is a relaxation / simplification of the problem\n* p2, l79-84: This paragraph needs to be improved. In particular, no references to existing works using the approach here are given. Are Deep Learning approaches used at all in this area? Furthermore, while at the start of the paragraph ReLUs are mentioned, at the end of the paragraph it is argued that Deep Learning approaches are problematic due to the use of smooth activation functions like sigmoids.\n* p2, l87: The technical terms \"traces\" and \"guards\" haven't been defined\n* p3, l91-92: Unclear, is variational inference used in the EM algorithm or is the EM algorithm regarded as a variational inference technique?\n* p3, Definition 1: The formulation is rather unusual and a bit imprecise. For example, $A_1,\\ldots,A_m$ are not $m$ different dynamical systems, but rather $m$ matrices inducing $m$ different linear dynamics. Furthermore, the last sentence \"We can generalize this to have output maps that also depend on the mode\" should be outside the definition.\n* p3, l121: The term \"observable\" has a technical meaning in systems theory, better to use \"not observed\" here\n* p3, l121-122: This sentence is unclear, more details on this approach are needed.\n* p3, l135: Typo, \"with\" instead of \"will\" needed. Furthermore, the definition of $E_n$ should be moved to the notation section.\n* p4, l136: Abbreviation SLS-ID not introduced\n* p4, l136-137: Better not to used colons in the sentence\n* p4, formula (2): Here $\\leq$ seems to be mean entrywise order. Either define this or rewrite the formula in an elementary manner (e.g. $\\forall j=1,\\ldots,m$ we have $(A_j)_{k,l}\\in[-\\gamma,\\gamma]$ for all $k,l=1,\\ldots,n$)\n* p4, l148: Is \"big M method\" a technical term in the field of switched linear systems?\n* p4, l149: Abbreviation MILP not introduced\n* p4, l150-151: There should be at least one reference for \"a common approach\"\n* p5, l158: This is not really a definition, but rather an algorithm description and should be formatted and labelled accordingly\n* p4, l159: Use $:=$ instead of $:$ \n* p4, l162-164: Are the outputs really the only two cases? This needs to be clarified\n* p4, l166: \"It is easy to show\". At least some details should be provided, maybe in the supplementary material. Furthermore, this statement should be made more precise. For example, what about the case $x_i=0$ and $\\|x_i'\\|_\\infty >> 0$?\n* p4, l167: What does \"ignoring the technicality related to the limit\" mean?\n* p5, Algorithm 1, l9: What does $\\hat{P}_j$ feasible mean? That this polyhedron is non-empty?\n    Furthermore, what happens if this condition is False for all $j=1,\\ldots,m$? \n* p5, l186: In which space are the polyhedra?\n* p7, l300: What is meant by traces? Trajectories?\n* p7, l303: Why $N=kT$? If the trajectories have $T$ time steps, then there are only $T-1$ transitions (and each data point correponds to a transition)\n* p7, l304: What are \"predicted and actual states\"? The setting introduced is identification, not prediction\n* p8, l314-321: Why was this setting chosen? Is there related work on such tasks? What would be the state of the art?\n* p8, l322-327: Details missing. For example, how are the JMM used?\n* p9, Figure 2: Better labelling required. In particular, axis labels are missing in both subplots\n* p9, l1360-361: It hasn't been explained how the approach can be extended using SARX models\n\n## Supplementary material\n* p13, l522-528: This paragraph could be made a bit easier to follow, e.g., to be more explicit.\n* p15, Section D: Details on the MILP method and the implementation of the proposed method should be provided here", "limitations": "The authors addressed adequately any potential societal or ethical implications.", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## 1. Originality\nThe problem setting seems to be established, but the proposed relaxation approach is new. The algorithm and its analysis are also novel, but based on standard arguments. In particular, in the relaxed problem setting no technical complications appear and no novel techniques seem to be necessary. The algorithm uses a pruning approach based on a treshold of the volume of certain polyhedrons which seems to be new, but due to lack of familiarity with this area this evaluation is not certain.\n\n\n## 2. Quality\nThe presentation needs to be improved. The article has language and grammar issues, as well as formatting problems (e.g., using a colon instead of an equality sign). Furthermore, some technical arguments could be explained better, e.g., in the Lemmas in Section 3 or Figure 2. Additionally, there should be more references. Finally, some relevant details on the experiments are missing. There seems to also some potential technical problems, cf. Questions section below. \n\n## 3. Clarity\nAs already mentioned, the presentation needs to be improved. For example, there is no pseudocode description of the whole algorithm (only some parts of it are presented in Algorithm 1). Some arguments and explanations are also difficult to follow.\n\n## 4. Significance\nThe identification of switched linear systems is an established and relevant area. The proposed algorithm seems to be working and computationally efficient (in the number of data points, not the state dimension and number of modes). However, the comparison with prior work and state of the art (in particular, in the experimental section) needs to be more thorough, so its difficult to confidently assess the overall significance.", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "WBnvWU9TW3", "reviewer": "Reviewer_3C4e", "summary": "This paper presents an algorithm for learning switched linear dynamical systems. An algorithm is proposed which is linear in the size of the problem and the proposed algorithm is demonstrated on some benchmark systems.", "questions": "Please try and answer these:\n1.Switched systems can also be modeled using complementarity constraints. It would be good to compare or contrast against those techniques. Do you think that your method could be more efficient than the methods using complementarity constraints? For example, see Contactnets by Post et al. or https://arxiv.org/pdf/2203.10013.pdf.\n2. The proposed method has been shown for autonomous systems only. Would it work for controlled systems, e.g. where x_k+1=Ax_k+Bu_k?\n3. One the practical problems of estimating or learning in hybrid systems is dealing with imbalance in data between different modes? To make any of these algorithms useful for, say robotic systems, that is a problem that needs to be solved. How do you suggest your algorithm will work if there is an imbalance in data between different modes?\n4. Have you tested you method when there are lots of possible modes in the system? A lot of practical systems like robotic systems usually have a lot of modes that can arise? How does your method scale with the number of modes? Please explain.\n", "limitations": "Some limitations have been described. I would like authors to think about the limitations I pointed out.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 1, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. the problem is well motivated, and still remains open in so many aspects.\n2. The analysis of the proposed method is interesting.\nWeaknesses:\n1. the comparison of the proposed work is weak and makes evaluation of the proposed work difficult.\n2. Switched systems can also be modeled using complementarity constraints. It would be good to compare or contrast against those techniques. \n3. The proposed method has been shown for autonomous systems only. Would it work for controlled systems, e.g. where x_k+1=Ax_k+Bu_k?\n", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}]

Author rebuttal:

Decision: Accept