## Linear Dynamics: Clustering without identification

Chloe Ching-Yun Hsu † University of California, Berkeley

Moritz Hardt † University of California, Berkeley

Michaela Hardt †

Amazon

## Abstract

Linear dynamical systems are a fundamental and powerful parametric model class. However, identifying the parameters of a linear dynamical system is a venerable task, permitting provably efficient solutions only in special cases. This work shows that the eigenspectrum of unknown linear dynamics can be identified without full system identification. We analyze a computationally efficient and provably convergent algorithm to estimate the eigenvalues of the state-transition matrix in a linear dynamical system.

When applied to time series clustering, our algorithm can efficiently cluster multidimensional time series with temporal offsets and varying lengths, under the assumption that the time series are generated from linear dynamical systems. Evaluating our algorithm on both synthetic data and real electrocardiogram (ECG) signals, we see improvements in clustering quality over existing baselines.

## 1 Introduction

Linear dynamical system (LDS) is a simple yet general model for time series. Many machine learning models are special cases of linear dynamical systems [Roweis and Ghahramani, 1999], including principal component analysis (PCA), mixtures of Gaussians, Kalman filter models, and hidden Markov models.

When the states are hidden, LDS parameter identification has provably efficient solutions only in special cases, see for example [Hazan et al., 2018, Hardt et al., 2018, Simchowitz et al., 2018]. In practice, the expectation-maximization (EM) algorithm [Ghahramani and Hinton, 1996] is often used

† This work was done at Google.

Proceedings of the 23 rd International Conference on Artificial Intelligence and Statistics (AISTATS) 2020, Palermo, Italy. PMLR: Volume 108. Copyright 2020 by the author(s).

for LDS parameter estimation, but it is inherently non-convex and can often get stuck in local minima [Hazan et al., 2018]. Even when full system identification is hard, is there still hope to learn meaningful information about linear dynamics without learning all system parameters ? We provide a positive answer to this question.

We show that the eigenspectrum of the state-transition matrix of unknown linear dynamics can be identified without full system identification. The eigenvalues of the state-transition matrix play a significant role in determining the properties of a linear system. For example, in two dimensions, the eigenvalues determine the stability of a linear dynamical system. Based on the trace and the determinant of the state-transition matrix, we can classify a linear system as a stable node, a stable spiral, a saddle, an unstable node, or an unstable spiral.

To estimate the eigenvalues, we utilize a fundamental correspondence between linear systems and Autoregressive-Moving-Average (ARMA) models. We establish bi-directional perturbation bounds to prove that two LDSs have similar eigenvalues if and only if their output time series have similar auto-regressive parameters. Based on a consistent estimator for the autoregressive model parameters of ARMA models [Tsay and Tiao, 1984], we propose a regularized iterated least-squares regression method to estimate the LDS eigenvalues. Our method runs in time linear in the sequence length T and converges to true eigenvalues at the rate O p ( T -1 / 2 ) .

As one application, our eigenspectrum estimation algorithm gives rise to a simple approach for time series clustering: First use regularized iterated least-squares regression to fit the autoregressive parameters; then cluster the fitted autoregressive parameters.

This simple and efficient clustering approach captures similarity in eigenspectrums, assuming each time series comes from an underlying linear dynamical system. It is a suitable similarity measure where the main goal for clustering is to characterize state-transition dynamics regardless of change of basis, particularly relevant when

there are multiple data sources with different measurement procedures. Our approach bypasses the challenge of LDS full parameter estimation, while enjoying the natural flexibility to handle multi-dimensional time series with time offsets and partial sequences.

To verify that our method efficiently learns system eigenvalues on synthetic and real ECG data, we compare our approach to existing baselines, including model-based approaches based on LDS, AR and ARMA parameter estimation, and PCA, as well as model-agnostic clustering approaches such as dynamic time warping [Cuturi and Blondel, 2017] and k-Shape [Paparrizos and Gravano, 2015].

Organization. We review LDS and ARMA models in Sec. 3. In Sec. 4 we discuss the main technical results around the correspondence between LDS and ARMA. Sec. 5 presents the regularized iterated regression algorithm, a consistent estimator of autoregressive parameters in ARMA models with applications to clustering. We carry out eigenvalue estimation and clustering experiments on synthetic data and real ECG data in Sec. 6. In the appendix, we describe generalizations to observable inputs and multidimensional outputs, and include additional simulation results.

## 2 Related Work

Linear dynamical system identification. The LDS identification problem has been studied since the 60s [Kalman, 1960], yet the theoretical bounds are still not fully understood. Recent provably efficient algorithms [Simchowitz et al., 2018, Hazan et al., 2018, Hardt et al., 2018, Dean et al., 2017] require setups that are not best-suited for time series clustering, such as assuming observable states and focusing on predition error instead of parameter recovery.

On recovering system parameters without observed states, Tsiamis et al. recently study a subspace identification algorithm with non-asymptotic O T ( -1 / 2 ) rate [Tsiamis and Pappas, 2019]. While our analysis does not provide finite sample complexity bounds, our simple algorithm achieves the same rate asymptotically.

Model-based time series clustering. Common model choices for clustering include Gaussian mixture models [Biernacki et al., 2000], autoregressive integrated moving average (ARIMA) models [Kalpakis et al., 2001], and hidden Markov models [Smyth, 1997]. Gaussian mixture models, ARIMA models, and hidden Markov models are all special cases of the more general linear dynamical system model [Roweis and Ghahramani, 1999].

Linear dynamical systems have used to cluster video trajectories [Chan and Vasconcelos, 2005,

Afsari et al., 2012, Vishwanathan et al., 2007], where the observed time series are higher dimensional than the hidden state dimension. Our work is motivated by the more challenging sitatuion with a single or a few output dimensions, common in climatology, energy consumption, finance, medicine, etc.

Compared to ARMA-parameter based clustering, our method only uses the AR half of the parameters which we show to enjoy more reliable convergence. We also differ from AR-model based clustering because fitting AR to an ARMA process results in biased estimates.

Autoregressive parameter estimation. Existing spectral analysis methods for estimating AR parameters in ARMA models include high-order Yule-Walker (HOYW), MUSIC, and ESPRIT [Stoica et al., 2005, Stoica et al., 1988]. Our method is based on iterated regression [Tsay and Tiao, 1984], a more flexible method for handling observed exogenous inputs (see Appendix C) in the ARMAX generalization.

## 3 Preliminaries

## 3.1 Linear dynamical systems

A discrete-time linear dynamical system (LDS) with parameters Θ = ( A,B,C,D ) receives inputs x , 1 · · · , x T ∈ R k , has hidden states h , 0 · · · , h T ∈ R n , and generates outputs y , 1 · · · , y T ∈ R m according to the following time-invariant recursive equations:

$$h _ { t } = A h _ { t - 1 } + B x _ { t } + \zeta _ { t } \\ y _ { t } = C h _ { t } + D x _ { t } + \xi _ { t }.$$

Assumptions. We assume that the stochastic noise ζ t and ξ t are diagonal Gaussians. We also assume the system is observable, i.e. C, CA, CA , 2 · · · , CA n -1 are linearly independent. When the LDS is not observable, the ARMA model for the output series can be reduced to lower AR order, and there is not enough information in the output series to recover all the full eigenspectrum.

The model equivalence theorem (Theorem 4.1) and the approximation theorem (Theorem 4.2) do not require any additional assumptions for any real matrix A . When additionally assuming A only has simple eigenvalues in C , i.e. each eigenvalue has multiplicity 1, we give a better convergence bound.

Distance between linear dynamical systems. With the main goal to characterize state-transition dynamics, we view systems as equivalent up to change of basis, and use the glyph[lscript] 2 distance of the spectrum of the transition matrix A , i.e. d (Θ 1 , Θ ) = 2 ‖ λ A ( 1 ) -λ A ( 2 ) ‖ 2 , where λ A ( 1 ) and λ A ( 2 ) are the spectrum of A 1 and A 2 in sorted order. This distance definition satisfies nonnegativity, identity, symmetry, and triangle inequality.

Two very different time series could still have small distance in eigenspectrum. This is by design to allow the flexiblity for different measurement procedures, which mathematically correspond to different C matricies. When there are multiple data sources with different measurement procedures, our approach can compare the underlying dynamics of time series across sources.

glyph[negationslash]

Jordan canonical basis. Every square real matrix is similar to a complex block diagonal matrix known as its Jordan canonical form (JCF). In the special case for diagonalizable matrices, JCF is the same as the diagonal form. Based on JCF, there exists a canonical basis { e i } consisting only of eigenvectors and generalized eigenvectors of A . A vector v is a generalized eigenvector of rank µ with corresponding eigenvalue λ if ( λI -A ) µ v = 0 and ( λI -A ) µ -1 v = 0 .

## 3.2 Autoregressive-moving-average models

The autoregressive-moving-average (ARMA) model combines the autoregressive (AR) model and the moving-average (MA) model. The AR part involves regressing the variable with respect its lagged past values, while the MA part involves regressing the variable against past error terms.

Autoregressive model. The AR model describes how the current value in the time series depends on the lagged past values. For example, if the GDP realization is high this quarter, the GDP in the next few quarters are likely high as well. An autoregressive model of order p , noted as AR ( p ) , depends on the past p steps,

$$y _ { t } = c + \Sigma _ { i = 1 } ^ { p } \varphi _ { i } y _ { t - i } + \epsilon _ { t },$$

where ϕ , 1 · · · , ϕ p are autoregressive parameters, c is a constant, and glyph[epsilon1] t is white noise.

When the errors are normally distributed, the ordinary least squares (OLS) regression is a conditional maximum likelihood estimator for AR models yielding optimal estimates [Durbin, 1960].

Moving-average model. The MA model, on the other hand, captures the delayed effects of unobserved random shocks in the past. For example, changes in winter weather could have a delayed effect on food harvest in the next fall. A moving-average model of order q , noted as MA ( q ) , depends on unobserved lagged errors in the past q steps,

$$y _ { t } = c + \epsilon _ { t } + \Sigma _ { i = 1 } ^ { q } \theta _ { i } \epsilon _ { t - i },$$

where θ , 1 · · · , θ q are moving-average parameters, c is a constant, and the errors glyph[epsilon1] t are white noise.

ARMA model. The autoregressive-moving-average (ARMA) model, denoted as ARMA ( p, q ) , merges AR ( p )

and MA ( q ) models to consider dependencies both on past time series values and past unpredictable shocks,

$$y _ { t } = c + \epsilon _ { t } + \Sigma _ { i = 1 } ^ { p } \varphi _ { i } y _ { t - i } + \Sigma _ { i = 1 } ^ { q } \theta _ { i } \epsilon _ { t - i }.$$

ARMAX model. ARMA can be generalized to autoregressive-moving-average model with exogenous inputs (ARMAX).

$$\underset { \omega \,. } { \text{case} } \quad y _ { t } = c + \epsilon _ { t } + \Sigma _ { i = 1 } ^ { p } \varphi _ { i } y _ { t - i } + \Sigma _ { i = 1 } ^ { q } \theta _ { i } \epsilon _ { t - i } + \Sigma _ { i = 0 } ^ { r } \gamma _ { i } x _ { t - i },$$

where { x t } is a known external time series, possibly multidimensional. In case x t is a vector, the parameters γ i are also vectors.

Estimating ARMA and ARMAX models is significantly harder than AR, since the model depends on unobserved variables and the maximum likelihood equations are intractable [Durbin, 1959, Choi, 2012]. Maximum likelihood estimation (MLE) methods are commonly used for fitting ARMA and ARMAX [Guo, 1996, Bercu, 1995, Hannan et al., 1980], but have converge issues. Although regression methods are also used in practice, OLS is a biased estimator for ARMA models [Tiao and Tsay, 1983].

Lag operator. We also introduce the lag operator, a concise way to describe ARMA models [Granger and Morris, 1976], defined as Ly t = y t -1 . The lag operator could be raise to powers, or form polynomials. For example, L y 3 t = y t -3 , and ( a L 2 2 + a L 1 + a 0 ) y t = a y 2 t -2 + a y 1 t -1 + a y 0 t . The lag polynomials can be multiplied or inverted. An AR( ) model can be p characterized by

$$\Phi ( L ) y _ { t } = c + \epsilon _ { t },$$

where Φ( L ) = 1 -ϕ L 1 -··· -ϕ L p p is a polynomial of the lag operator L of degree p . For example, any AR(2) model can be described as (1 -ϕ L 1 -ϕ L 2 2 ) y t = c + glyph[epsilon1] t .

Similarly, an MA( q ) can be characterized by a polynomial Ψ( L ) = θ L q q + · · · + θ L 1 +1 of degree q ,

$$y _ { t } = c + \Psi ( L ) \epsilon _ { t }.$$

For example, for an MA(2) model the equation would be y t = c +( θ L 2 2 + θ L 1 +1) glyph[epsilon1] t .

Merging the two and adding dependency to exogenous input, we can write an ARMAX( p, q, r ) model as

$$\Phi ( L ) y _ { t } = c + \Psi ( L ) \epsilon _ { t } + \Gamma ( L ) x _ { t } \text{ \quad \ \ } ( 2 )$$

where Φ Ψ , , and Γ are polynomials of degree p, q and r . When the exogenous time series x t is multidimensional, Γ( L ) is a vector of degreer polynomials.

## 4 Learning eigenvalues without system identification

This section provides theoretical foundations for learning LDS eigenvalues from autoregressive parameters without full system identification. While general model equivalence between LDS and ARMA(X) is known [Åström and Wittenmark, 2013, Kailath, 1980], we provide detailed analysis of the exact correspondence between the LDS characteristic polynomial and the ARMA(X) autoregressive parameters along with perturbation bounds.

## 4.1 Model equivalence

We show that the output series from any LDS can be seen as generated by an ARMAX model, whose AR parameters contain full information about the LDS eigenvalues.

Theorem 4.1. Let y t ∈ R m be the outputs from a linear dynamical system with parameters Θ = ( A,B,C,D ) , hidden dimension n , and inputs x t ∈ R k . Each dimension of y t can be generated by an ARMAX( n, n, n -1 ) model, whose autoregressive parameters ϕ , 1 · · · , ϕ n can recover the characteristic polynomial of A by χ A ( λ ) = λ n -ϕ λ 1 n -1 -··· -ϕ n .

In the special case where the LDS has no external inputs, the ARMAX model is an ARMA ( n, n ) model.

## See Appendix A for the full proof.

As a high-level proof sketch: We first analyze the hidden state projected to (generalized) eigenvector directions in Lemma A.2. We show that for a (generalized) eigenvector e of the adjoint A ∗ of the transition matrix with eigenvalue λ and rank µ , the time series obtained from applying the lag operator polynomial (1 -λL ) µ to 〈 h , e t 〉 can be expressed as a linear combination of the past k inputs x , t · · · , x t -k +1 . Since A is real-valued, A and its adjoint A ∗ share the same characteristic polynomial χ A .

We then consider the lag operator polynomial χ † A ( L ) = L χ n A ( L -1 ) , and show that the time series obtained from applying χ † A ( L ) applied to any (generalized) eigenvector direction is a linear combination of the past k inputs. We use this on the Jordan canonical basis for A ∗ that consists of (generalized) eigenvectors. From there, we conclude χ † A is the autoregressive lag polynomial that contains the autoregressive coefficients.

The converse of Theorem 4.1 also holds. An ARMA( p, q ) model can be seen as a ( p + ) q -dimensional LDS where the state encodes the relevant past values and error terms.

Corollary 4.1. The output series of two linear dynam- ical systems have the same autoregressive parameters if and only if they have the same non-zero eigenvalues with the same multiplicities.

Proof. By Theorem 4.1, the autoregressive parameters are determined by the characteristic polynomial. Two LDSs of the same dimension have the same autoregressive parameters if and only if they have the same characteristic polynomials, and hence the same eigenvalues with the same multiplicities. Two LDSs of different dimensions n 1 &lt; n 2 can have the same autoregressive parameters if and only if χ A 1 ( λ ) = χ A 2 ( λ λ ) n 2 -n 1 and ϕ n 1 +1 = · · · = ϕ n 2 = 0 , in which case they have the same non-zero eigenvalues with same multiplicities.

It is possible for two LDSs with different dimensions to have the same AR coefficients, if the higher-dimensional system has additional zero eigenvalues. Whether overparameterized models indeed learn additional zero eigenvalues requires further empirical investigation.

## 4.2 Approximation theorems for LDS eigenvalues

We show that small error in the AR parameter estimation guarantees a small error in the eigenvalue estimation. This implies that an effective estimation algorithm for the AR parameters in ARMAX models leads to effective estimation of LDS eigenvalues.

## General (1 /n ) -exponent bound

Theorem 4.2. Let y t be the outputs from an n -dimensional linear dynamical system with parameters Θ = ( A,B,C,D ) , eigenvalues λ , 1 · · · , λ n , and hidden inputs. Let ˆ Φ = (ˆ ϕ , 1 · · · , ϕ ˆ n ) be the estimated autoregressive parameters for { y t } with error ‖ ˆ Φ -Φ ‖ = glyph[epsilon1] , and let r , 1 · · · , r n be the roots of the polynomial 1 -ˆ ϕ z 1 -··· -ˆ ϕ z n n .

Assuming the LDS is observable, the roots converge to the true eigenvalues with convergence rate O ( glyph[epsilon1] 1 /n ) . If all eigenvalues of A are simple (no multiplicity), then the convergence rate is O ( glyph[epsilon1] ) .

Without additional assumptions, the 1 n -exponent in the above general bound is tight. As an example, z 2 -glyph[epsilon1] has roots z ± √ glyph[epsilon1] . The general phenomenon that a root with multiplicity m can split into m roots at rate O glyph[epsilon1] ( m ) is related to the regular splitting property [Hryniv and Lancaster, 1999, Lancaster et al., 2003].

## Linear bound for simple eigenvalues

Under the additional assumption that all the eigenvalues are simple (no multiplicity), we derive a better O glyph[epsilon1] ( ) bound instead of O glyph[epsilon1] ( 1 /n ) . We show small perturbation

in AR parameters results in small perturbation in companion matrix, and small perturbation in companion matrix results in small perturbation in eigenvalues.

We defer the full proofs to Appendix B, but describe the proof ideas here.

For a monic polynomial Φ( u ) = z n + ϕ z 1 n -1 + · · · + ϕ n -1 z + ϕ n , the companion matrix , also known as the controllable canonical form in control theory, is the square matrix

$$C ( \Phi ) = \begin{bmatrix} 0 & 0 & \dots & 0 & - \varphi _ { n } \\ 1 & 0 & \dots & 0 & - \varphi _ { n - 1 } \\ 0 & 1 & \dots & 0 & - \varphi _ { n - 2 } \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & - \varphi _ { 1 } \end{bmatrix}.$$

The matrix C (Φ) is the companion in the sense that it has Φ as its characteristic polynomial.

In relation to an autoregressive AR( p ) model, the companion matrix corresponds to the transition matrix in the linear dynamical system when we encode the values form the past p lags as a p -dimensional state

$$h _ { t } = \left [ y _ { t - p + 1 } \ \cdots \ y _ { t - 1 } \ y _ { t } \right ] ^ { T }.$$

If y t = ϕ y 1 t -1 + · · · + ϕ y p t -p , then h t =

$$& \wedge \ s _ { n } \ r _ { 1 } y _ { n - 1 } \cdot \cdot \cdot \cdot \ r p y _ { n - p }, \dots \cdots \cdot \cdot \cdot \cdot \\ & \begin{bmatrix} y _ { t - p + 1 } \\ y _ { t - p + 2 } \\ \dots \\ y _ { t - 1 } \\ y _ { t } \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 & \dots & 0 \\ 0 & 0 & 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & 1 \\ \varphi _ { p } & \varphi _ { p - 1 } & \varphi _ { p - 2 } & \dots & \varphi _ { 1 } \end{bmatrix} \begin{bmatrix} y _ { t - p } \\ \dots \\ y _ { t - 2 } \\ y _ { t - 1 } \end{bmatrix} & \text{the $a$} \\ & = C ( - \Phi ) ^ { T } h _ { t - 1 }. \\ & & & & ( 3 ) & \text{bias} \end{bmatrix}$$

We then use matrix eigenvalue perturbation theory results on the companion matrix for the desired bound.

Lemma 4.1 (Theorem 6 in [Lancaster et al., 2003]) . Let L λ, glyph[epsilon1] ( ) be an analytic matrix function with semisimple eigenvalue λ 0 at glyph[epsilon1] = 0 of multiplicity M . Then there are exactly M eigenvalues λ i ( glyph[epsilon1] ) of L λ, glyph[epsilon1] ( ) for which λ i ( glyph[epsilon1] ) → λ 0 as glyph[epsilon1] → 0 , and for these eigenvalues

$$\lambda _ { i } ( \epsilon ) = \lambda _ { 0 } + \lambda _ { i } ^ { \prime } \epsilon + o ( \epsilon ). \quad \quad ( 4 ) \quad \begin{smallmatrix} & W \epsilon \\ & \ s i o \end{smallmatrix}$$

## Explicit bound on condition number

When the LDS has all simple eigenvalues, we provide a more explicit bound on the condition number.

Theorem 4.3. In the same setting as above in Theorem 4.2, when all eigenvalues of A are simple, | r j -λ j | ≤ κglyph[epsilon1] + ( o glyph[epsilon1] 2 ) , then the condition number κ is bounded by

glyph[negationslash]

$$\frac { 1 } { \prod _ { k \neq j } | \lambda _ { j } - \lambda _ { k } | } \leq \kappa \leq \\ \frac { \sqrt { n } } { \prod _ { k \neq i } | \lambda _ { j } - \lambda _ { k } | } \left ( \max ( 1, | \lambda _ { j } | ) \right ) ^ { n - 1 } ( 1 + \rho ( A ) ^ { 2 } ) ^ { \frac { n - 1 } { 2 } },$$

glyph[negationslash]

$$\mathfrak { e } _ { \shortparallel } \frac { \sqrt { n } } { \prod _ { k \neq j } | \lambda _ { j } - \lambda _ { k } | } \left ( \max ( 1, | \lambda _ { j } | ) \right ) ^ { n - 1 } ( 1 + \rho ( A ) ^ { 2 } ) ^ { \frac { n - 1 } { 2 } },$$

where ρ A ( ) is the spectral radius, i.e. largest absolute value of its eigenvalues.

glyph[negationslash]

In particular, when ρ C ( ) ≤ 1 , i.e. when the matrix is Lyapunov stable, then the absolute difference between the root from the auto-regressive method and the eigenvalue is bounded by | r j -λ j | ≤ √ n ( √ 2) n -1 ∏ k = j | λ j -λ k | glyph[epsilon1] + ( o glyph[epsilon1] 2 ) .

In Appendix B, we derive the explicit formula in Theorem 4.3 by conjugating the companion matrix by a Vandermonde matrix to diagonalize it and invoking the explicit inverse formula of Vandermonde matrices.

## 5 Estimation of ARMA autoregressive parameters

In general, learning ARMA models is hard, since the output series depends on unobserved error terms. Fortunately, for our purpose we are only interested in the autoregressive parameters, that are easier to learn since the past values of the time series are observed.

The autoregressive parameters in an ARMA( p, q ) model are not equivalent to the pure AR( p ) parameters for the same time series. For AR( ) models, p ordinary least squares (OLS) regression is a consistent estimator of the autoregressive parameters [Lai and Wei, 1983]. However, for ARMA( p, q ) models, due to the serial correlation in the error term glyph[epsilon1] t + ∑ q i =1 θ glyph[epsilon1] i t -i , the OLS estimates for autoregressive parameters can be biased [Tiao and Tsay, 1983].

Regularized iterated regression. Iterated regression [Tsay and Tiao, 1984] is a consistent estimator for the AR parameters in ARMA models. While iterated regression is theoretically well-grounded, it tends to over-fit and results in excessively large parameters. To avoid over-fitting, we propose a slight modification with regularization, which keeps the same theoretical guarantees and yields better practical performance.

We also generalize the method to handle multidimensional outputs from the LDS and observed inputs by using ARMAX instead of ARMA models, as described in details in Appendix C as Algorithm 2.

The i -th iteration of the regression only uses error terms from the past i lags. The initial iteration is an ARMA( n, 0 ) regression, the first iteration is an ARMA( n, 1 ) regression, and so forth until ARMA( n, n ) in the last iteration.

Algorithm 1: Regularized iterated regression for au- toregressive parameter estimation

Input: Time series { y t } T t =1 , target hidden state dimension n , and regularization coefficient α . Initialize error term estimates ˆ glyph[epsilon1] t = 0 for t = 1 , . . . , T ; for i = 0 , · · · , n do Perform glyph[lscript] 2 -regularized least squares regression to estimate ˆ ϕ j , ˆ θ j , and ˆ c in y t = ∑ n j =1 ˆ ϕ y j t - j + ∑ i j =1 ˆ ˆ θ glyph[epsilon1] j t - j +ˆ c with regularization strength α only on the ˆ θ j terms; Update ˆ glyph[epsilon1] t to be the residuals from the most recent regression;

end

Return ˆ ϕ , 1 · · · , ϕ ˆ n .

Time complexity. The iterated regression involves n +1 steps of least squares regression each on at most 2 n +1 variables. Therefore, the total time complexity of Algorithm 1 is O n T ( 3 + n 4 ) , where T is the sequence length and n is the hidden state dimension.

Convergence rate. The consistency and the convergence rate of the estimator is analyzed in [Tsay and Tiao, 1984]. Adding regularization does not change the asymptotic property of the estimator.

Theorem 5.1 ([Tsay and Tiao, 1984]) . Suppose that y t is an ARMA( p, q ) process, stationary or not. The estimated autoregressive parameters ˆ = ( ˆ Φ ϕ , 1 · · · , ϕ ˆ n ) from iterated regression converges in probability to the true parameters with rate

$$\hat { \Phi } = \Phi + O _ { p } ( T ^ { - 1 / 2 } ),$$

or more explicitly, convergence in probability means that for all glyph[epsilon1] , lim T → inf Pr( T 1 / 2 | ˆ Φ -Φ | &gt; glyph[epsilon1] ) = 0 .

## 5.1 Applications to clustering

The task of clustering depends on an appropriate distance measure for the clustering purpose. When there are multiple data sources for time series with comparable dynamics but measured by different measurement procedures, one might hope to cluster time series based only on the state-transition dynamics.

We can observe from the LDS definition (1) that two LDSs with parameters ( A,B,C,D ) and ( A , B , C , D ′ ′ ′ ′ ) are equivalent if A ′ = P -1 AP,B ′ = P -1 B,C ′ = CP,D ′ = DP , and h ′ t = P -1 h t under change of basis by some non-singular matrix P . Therefore, to capture the state-transition dynamics while allowing for flexibility in the measurement matrix C , the distance measure should be invariant under change of basis. We choose to use the eigenspectrum distance between statetransition matrices, a natural distance choice that is invariant under change of basis.

In previous sections, our theoretical analysis shows that AR parameters in ARMA time series models can effectively estimate the eigenspectrum of underlying LDSs. We therefore propose a simple time series clustering algorithm: 1) first use iterated regression to estimate the autoregressive parameters in ARMA models for each times series, and 2) then apply any standard clustering algorithm such as K-means on the distance between autoregressive parameters.

Our method is very flexible. It handles multidimensional data, as Theorem 4.1 suggests that any output series from the same LDS should share the same autoregressive parameters. It can also handle exogenous inputs as illustrated in Algorithm 2 in Appendix C. It is scale, shift, and offset invariant, as the autoregressive parameters in ARMA models are. It accommodates missing values in partial sequences as we can still perform OLS after dropping the rows with missing values. It also allows sequences to have different lengths, and could be adapted to handle sequences with different sampling frequencies, as the compound of multiple steps of LDS evolution is still linear.

## 6 Experiments

We experimentally evaluate the quality and efficiency of the clustering from our method and compare it to existing baselines. The source code is available online at https://github.com/chloechsu/ldseig .

## 6.1 Methods

- · ARMA: K-means on AR parameters in ARMA( n, n ) model estimated by regularized iterated regression as we proposed in Algorithm 1.
- · ARMA\_MLE: K-means on AR parameters in ARMA( n, n ) model estimated by the MLE method using statsmodels [Seabold and Perktold, 2010].
- · AR: K-means on AR parameters in AR( n ) model estimated by OLS using statsmodels.
- · LDS: K-means on estimated LDS eigenvalues. We estimate the LDS eigenvalues with the pylds package [Johnson and Linderman, 2018], with 100 EM steps initialized by 10 Gibbs iterations.
- · k-Shape: A shape-based time series clustering method [Paparrizos and Gravano, 2015], using the tslearn package [Tavenard, 2017].
- · DTW: K-medoids on dynamic time warping distance, using the dtaidistance [Meert, 2018] and pyclustering [Novikov, 2019] packages.
- · PCA: K-means on the first n PCA components, using sklearn [Pedregosa et al., 2011].

| # Clusters   | Method                               | Adj.                                                                                     | Mutual Info.                                                                                                                                                                                                                                                                                                                                         | Rand                                                                                                                                                                                       | V-measure                                                                                                                                                                                                                |
|--------------|--------------------------------------|------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| 2            | AR ARMA ARMA_MLE DTW k-Shape LDS PCA | 0.06 0.13 0.02 0.02 0.03 0.09 -0.00 0.11 0.18 0.04 0.04                                  | Adj. Score (0.04-0.08) 0.07 (0.05-0.09) (0.11-0.16) 0.16 (0.13-0.19) (0.01-0.03) 0.02 (0.01-0.03) (0.01-0.03) 0.02 (0.01-0.03) (0.02-0.04) 0.03 (0.02-0.05) (0.06-0.12) 0.09 (0.06-0.12) (-0.00-0.00) -0.00 (-0.00-0.00) 0.09 (0.07-0.10) 0.16 (0.14-0.18) 0.04 (0.03-0.05) 0.03 (0.02-0.03) 0.04 (0.04-0.05) (0.18-0.23) 0.17 0.00 0.11 (0.10-0.12) | 0.07 (0.05-0.09) 0.14 (0.11-0.16) 0.03 (0.02-0.04) 0.03 (0.02-0.04) 0.04 (0.03-0.05) 0.10 (0.07-0.12) 0.02 (0.02-0.02) 0.12 (0.11-0.14) 0.19 (0.17-0.21) 0.06 (0.05-0.07) 0.07 (0.06-0.08) | Runtime (secs) 1.09 (1.01-1.17) 0.44 (0.40-0.47) 70.64 (68.01-73.28) 6.60 (6.34-6.86) 28.58 (25.15-32.01) 341.08 (328.24-353.93) 0.45 (0.43-0.47) 1.02 (0.93-1.10) 0.42 (0.38-0.46) 72.58 (69.63-75.52) 6.63 (6.37-6.90) |
| 3            | AR ARMA ARMA_MLE DTW k-Shape LDS PCA | (0.09-0.12) (0.16-0.20) (0.03-0.05) (0.03-0.05) 0.06 (0.05-0.07) 0.20 0.00 (-0.00-0.00)  | (0.15-0.20) (-0.00-0.00)                                                                                                                                                                                                                                                                                                                             | 0.08 (0.07-0.09) 0.22 (0.19-0.24) 0.04 (0.04-0.04)                                                                                                                                         | 40.10 (34.65-45.55) 338.67 (325.66-351.67) 0.47 (0.45-0.49) 0.91 (0.83-1.00) 0.40 (0.35-0.45)                                                                                                                            |
|              | AR ARMA ARMA_MLE DTW                 | (0.16-0.18) (0.21-0.23) (0.07-0.09) (0.04-0.06) (0.07-0.09) (0.23-0.26)                  | 0.15 (0.14-0.16) 0.05 (0.04-0.05) 0.03 (0.02-0.03) 0.05 (0.04-0.05) 0.17 (0.16-0.18) 0.00 (0.00-0.01)                                                                                                                                                                                                                                                | 0.22 (0.21-0.23) 0.26 (0.25-0.28) 0.14 (0.13-0.15) 0.11 (0.10-0.12)                                                                                                                        | 74.00 (71.70-76.30) 6.20 (5.86-6.55) 97.83 (84.97-110.68)                                                                                                                                                                |
| 5            | k-Shape                              | 0.17 0.22 0.08 0.05 0.08 0.25                                                            |                                                                                                                                                                                                                                                                                                                                                      | 0.14 (0.13-0.15) 0.29 (0.28-0.30) 0.08 (0.07-0.08)                                                                                                                                         | 321.40                                                                                                                                                                                                                   |
| 10           | LDS PCA AR ARMA                      | 0.01 (0.01-0.01) 0.22 0.24 0.11 0.06 0.08 0.23                                           | 0.11 (0.11-0.12) 0.14 0.06 (0.05-0.06) 0.03 (0.02-0.03) 0.04 (0.03-0.04) 0.13 (0.12-0.13)                                                                                                                                                                                                                                                            | 0.38 (0.37-0.38) 0.39 (0.39-0.40) 0.29 (0.28-0.30) 0.25 (0.24-0.25) 0.26 (0.26-0.27) 0.39 (0.38-0.40) 0.14 (0.13-0.15)                                                                     | (304.17-338.64) 0.52 (0.49-0.55) 0.87 (0.77-0.98) 0.42 (0.37-0.48) 63.39 (60.49-66.30) 5.51 (5.03-5.98) 108.79 (93.41-124.18) 277.45 (253.38-301.51) 0.47 (0.44-0.51)                                                    |
|              | ARMA_MLE DTW k-Shape LDS PCA         | (0.21-0.22) (0.23-0.25) (0.10-0.12) (0.06-0.07) (0.07-0.08) (0.22-0.24) 0.02 (0.01-0.02) | (0.13-0.15) 0.01 (0.00-0.01)                                                                                                                                                                                                                                                                                                                         |                                                                                                                                                                                            |                                                                                                                                                                                                                          |

Table 1: Performance of clustering 100 random 2-dimensional LDSs based on their output series of length 1000, with 95% confidence intervals from 100 trials. AMI, Adj. Rand, and V-measure are the adjusted mutual information score, adjusted Rand score, and V-measure between ground truth cluster labels and learned clusters. The runtime is on an instance with 12 CPUs and 40 GB memory running Ubuntu 18.

Figure 1: Absolute glyph[lscript] 2 -error in eigenvalue estimation for 2-dimensional and 3-dimensional LDSs, with 95% confidence interval from 500 trials.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01039v3-figure-71a8dd00-d0f0-46a3-ab0a-9b64a37cfbdc.jpg)

Table 2: Clustering performance on electrocardiogram (ECG) data separating segments of normal sinus rhythm from supraventricular tachycardia. 95% Confidence intervals are from 100 bootstrapped samples of 50 series.

| Method   | AMI              |   Adj. | Rand Score   | V-measure        |
|----------|------------------|--------|--------------|------------------|
| ARMA     | 0.12 (0.10-0.13) |   0.12 | (0.11-0.14)  | 0.14 (0.12-0.16) |
| AR       | 0.10 (0.09-0.12) |   0.1  | (0.09-0.12)  | 0.13 (0.11-0.14) |
| PCA      | 0.04 (0.03-0.05) |   0.02 | (0.02-0.03)  | 0.07 (0.06-0.08) |
| LDS      | 0.09 (0.07-0.10) |   0.11 | (0.10-0.13)  | 0.10 (0.09-0.11) |
| k-Shape  | 0.08 (0.07-0.10) |   0.09 | (0.07-0.11)  | 0.10 (0.09-0.12) |
| DTW      | 0.03 (0.02-0.04) |   0.02 | (0.01-0.03)  | 0.06 (0.05-0.07) |

## 6.2 Metrics of Cluster Quality

We measure cluster quality using three metrics in sklearn: V-measure [Rosenberg and Hirschberg, 2007], adjusted mutual information [Vinh et al., 2010], and adjusted Rand score [Hubert and Arabie, 1985].

## 6.3 Simulation

Dataset. We generate LDSs representing cluster centers with random matrices of i.i.d. Gaussian entries. From the cluster centers, we derive LDSs that are close to the cluster centers. From each LDS, we generate time series of length 1000 by drawing inputs from standard Gaussians and adding noise to the output sampled from N (0 0 01 , . 2 ) . More details in Appendix D.1.

Clustering performance. The iterated ARMA regression method and the LDS method yield the best clustering quality, while the iterated ARMA regression method is significantly faster. These results hold up for choices of different cluster quality metrics and number of clusters.

Eigenvalue Estimation. Good clustering results rely on good approximations of the LDS eigenvalue distance. Our analyses in Theorem 5.1 and Theorem 4.2 proved that the iterative ARMA regression algorithm can learn the LDS eigenvalues with converge rate O p ( T -1 / 2 ) . In Figure 1, we see that the observed convergence rate in simulations roughly matches the theoretical bound.

Each EM step in LDS runs in O n T ( 3 ) . When running a constant number of EM steps, LDS has the same total complexity as iterated ARMA . We chose 100 steps based on empirical evaluation of convergence for sequence length 1000. However, longer sequences may need more EM steps to converge, which would explain the increase in LDS eigenvalue estimation error for sequence length 50000 in Figure 1. Depending on different implementations and initialization schemes, it is possible that the LDS performance can be further optimized.

ARMA and LDS have comparable eigenvalue estimation error for most configurations. While the pure AR approach also gives comparable estimation error on relatively short sequences, its estimation is biased, and the error does not go down as sequence length increases.

## 6.4 Real-world ECG data

While our simulation results show the efficacy of our method, the data generation process satisfy assumptions that may not hold on real data. As a proof-ofconcept, we also test our method on real electrocardiogram (ECG) data.

Dataset. The MIT-BIH [Moody and Mark, 2001] dataset in PhysioNet [Goldberger et al., 2000] is the most common dataset for evaluating algorithms for ECG data [De Chazal et al., 2004, Yeh et al., 2012, Özbay et al., 2006, Ceylan et al., 2009].

It contains 48 half-hour recordings collected at the Beth Israel Hospital between 1975 and 1979. Each twochannel recording is digitized at a rate of 360 samples per second per channel. 15 distinct rhythms are annotated in recordings including abnormalities of cardiac rhythm (arrhythmias) by two cardiologists.

Detecting cardiac arrhythmias has stimulated research and product applications such as Apple's FDA-approved detection of atrial fibrillation [Turakhia et al., 2018]. ECG data have been modeled with AR and ARIMA models [Kalpakis et al., 2001, Corduas and Piccolo, 2008, Ge et al., 2002], and more recently convolutional neural networks [Hannun et al., 2019].

We bootstrap 100 samples of 50 time series; each bootstrapped sample consists of 2 labeled clusters: 25 series with supraventricular tachycardia and 25 series with normal sinus rhythm. Each series has length 500 which adequately captures a complete cardiac cycle. We set the ARMA glyph[lscript] 2 -regularization coefficient to be 0.01, chosen based on our simulation results.

Results. Comparing methods outlined in Section 6.1, Table 2 shows that our method achieves the best quality closely followed by the AR and LDS methods, according to adjusted mutual information, adjusted Rand score and V-measure, while being computationally efficient.

## 7 Conclusion

We give a fast, simple, and provably effective method to estimate linear dynamical system (LDS) eigenvalues based on system outputs. The algorithm combines statistical techniques from the 80's with our insights on the correspondence between LDSs and ARMA models.

As a proof-of-concept, we apply the eigenvalue estimation algorithm to time series clustering. The resulting clustering approach is flexible to handle varying lengths, temporal offsets, as well as multidimensional inputs and outputs. Our efficient algorithm yields high quality clusters in simulations and on real ECG data.

While LDSs are general models encompassing mixtures of Gaussian and hidden Markov models, they may not fit all applications. It would be interesting to extend the analysis to non-linear models, and to consider model overparameterization and misspecification.

Acknowledgments We thank the anonymous reviewers for thoughtful comments, and Scott Linderman, Andrew Dai, and Eamonn Keogh for helpful guidance.

## References

[Afsari et al., 2012] Afsari, B., Chaudhry, R., Ravichandran, A., and Vidal, R. (2012). Group action induced distances for averaging and clustering linear dynamical systems with applications to the analysis of dynamic scenes. In 2012 IEEE Conference on Computer Vision and Pattern Recognition , pages 2208-2215. IEEE.

[Åström and Wittenmark, 2013] Åström, K. J. and Wittenmark, B. (2013). Computer-controlled systems: theory and design . Courier Corporation.

[Beauzamy, 1999] Beauzamy, B. (1999). How the roots of a polynomial vary with its coefficients: A local quantitative result. Canadian Mathematical Bulletin , 42(1):3-12.

[Bercu, 1995] Bercu, B. (1995). Weighted estimation and tracking for armax models. SIAM Journal on Control and Optimization , 33(1):89-106.

[Biernacki et al., 2000] Biernacki, C., Celeux, G., and Govaert, G. (2000). Assessing a mixture model for clustering with the integrated completed likelihood. IEEE transactions on pattern analysis and machine intelligence , 22(7):719-725.

[Ceylan et al., 2009] Ceylan, R., Özbay, Y., and Karlik, B. (2009). A novel approach for classification of ecg arrhythmias: Type-2 fuzzy clustering neural network. Expert Systems with Applications , 36(3):67216726.

[Chan and Vasconcelos, 2005] Chan, A. B. and Vasconcelos, N. (2005). Probabilistic kernels for the classification of auto-regressive visual processes. In 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR'05) , volume 1, pages 846-851. IEEE.

[Choi, 2012] Choi, B. (2012). ARMA model identification . Springer Science &amp; Business Media.

[Corduas and Piccolo, 2008] Corduas, M. and Piccolo, D. (2008). Time series clustering and classification by the autoregressive metric. Computational Statistics &amp; Data Analysis , 52:1860-1872.

[Cuturi and Blondel, 2017] Cuturi, M. and Blondel, M. (2017). Soft-dtw: a differentiable loss function for time-series. In Proceedings of the 34th International Conference on Machine Learning-Volume 70 , pages 894-903. JMLR. org.

[De Chazal et al., 2004] De Chazal, P., O'Dwyer, M., and Reilly, R. B. (2004). Automatic classification of heartbeats using ecg morphology and heartbeat interval features. IEEE transactions on biomedical engineering , 51(7):1196-1206.

[Dean et al., 2017] Dean, S., Mania, H., Matni, N., Recht, B., and Tu, S. (2017). On the sample complexity of the linear quadratic regulator. arXiv preprint arXiv:1710.01688 .

[Durbin, 1959] Durbin, J. (1959). Efficient estimation of parameters in moving-average models. Biometrika , 46(3/4):306-316.

[Durbin, 1960] Durbin, J. (1960). Estimation of parameters in time-series regression models. Journal of the Royal Statistical Society: Series B (Methodological) , 22(1):139-153.

[El-Mikkawy, 2003] El-Mikkawy, M. E. (2003). Explicit inverse of a generalized vandermonde matrix. Applied mathematics and computation , 146(2-3):643-651.

[Ge et al., 2002] Ge, D., Srinivasan, N., and M Krishnan, S. (2002). Cardiac arrhythmia classification using autoregressive modeling. Biomedical engineering online , 1:5.

[Ghahramani and Hinton, 1996] Ghahramani, Z. and Hinton, G. E. (1996). Parameter estimation for linear dynamical systems. Technical report, Technical Report CRG-TR-96-2, University of Totronto, Dept. of Computer Science.

[Goldberger et al., 2000] Goldberger, A. L., Amaral, L. A. N., Glass, L., Hausdorff, J. M., Ivanov, P. C., Mark, R. G., Mietus, J. E., Moody, G. B., Peng, C.-K., and Stanley, H. E. (2000). PhysioBank, PhysioToolkit, and PhysioNet: Components of a new research resource for complex physiologic signals. Circulation , 101(23):e215-e220. Circulation Electronic Pages: http://circ.ahajournals.org/content/101/23/e215.full PMID:1085218; doi: 10.1161/01.CIR.101.23.e215.

[Granger and Morris, 1976] Granger, C. W. and Morris, M. J. (1976). Time series modelling and interpretation. Journal of the Royal Statistical Society: Series A (General) , 139(2):246-257.

[Greenbaum et al., 2019] Greenbaum, A., Li, R.-c., and Overton, M. L. (2019). First-order perturbation theory for eigenvalues and eigenvectors. arXiv preprint arXiv:1903.00785 .

[Guo, 1996] Guo, L. (1996). Self-convergence of weighted least-squares with applications to stochastic adaptive control. IEEE transactions on automatic control , 41(1):79-89.

[Hannan et al., 1980] Hannan, E. J., Dunsmuir, W. T., and Deistler, M. (1980). Estimation of vector armax models. Journal of Multivariate Analysis , 10(3):275295.

[Hannun et al., 2019] Hannun, A. Y., Rajpurkar, P., Haghpanahi, M., Tison, G. H., Turakhia, M. P., Bourn, C., and Ng, A. Y. (2019). Cardiologist-level arrhythmia detection and classification in ambulatory electrocardiograms using a deep neural network. Nature Medicine , 15.

[Hardt et al., 2018] Hardt, M., Ma, T., and Recht, B. (2018). Gradient descent learns linear dynamical systems. The Journal of Machine Learning Research , 19(1):1025-1068.

[Hazan et al., 2018] Hazan, E., Lee, H., Singh, K., Zhang, C., and Zhang, Y. (2018). Spectral filtering for general linear dynamical systems. In Advances in Neural Information Processing Systems , pages 4634-4643.

[Hryniv and Lancaster, 1999] Hryniv, R. and Lancaster, P. (1999). On the perturbation of analytic matrix functions. Integral Equations and Operator Theory , 34(3):325-338.

[Hubert and Arabie, 1985] Hubert, L. and Arabie, P. (1985). Comparing partitions. Journal of classification , 2(1):193-218.

[Ipsen and Rehman, 2008] Ipsen, I. C. and Rehman, R. (2008). Perturbation bounds for determinants and characteristic polynomials. SIAM Journal on Matrix Analysis and Applications , 30(2):762-776.

[Johnson and Linderman, 2018] Johnson, M. and Linderman, S. (2018). Pylds: Bayesian inference for linear dynamical systems. https://github.com/ mattjj/pylds .

[Kailath, 1980] Kailath, T. (1980). Linear systems , volume 156. Prentice-Hall Englewood Cliffs, NJ.

[Kalman, 1960] Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Journal of basic Engineering , 82(1):35-45.

[Kalpakis et al., 2001] Kalpakis, K., Gada, D., and Puttagunta, V. (2001). Distance measures for effective clustering of arima time-series. In Proceedings 2001 IEEE international conference on data mining , pages 273-280. IEEE.

[Kalpakis et al., 2001] Kalpakis, K., Gada, D., and Puttagunta, V. (2001). Distance measures for effective clustering of arima time-series. In Proceedings 2001 IEEE International Conference on Data Mining , pages 273-280.

[Lai and Wei, 1983] Lai, T. and Wei, C. (1983). Asymptotic properties of general autoregressive models and strong consistency of least-squares estimates of their parameters. Journal of multivariate analysis , 13(1):1-23.

[Lancaster et al., 2003] Lancaster, P., Markus, A., and Zhou, F. (2003). Perturbation theory for analytic matrix functions: the semisimple case. SIAM Journal on Matrix Analysis and Applications , 25(3):606-626.

[Meert, 2018] Meert, Wannes, e. a. (2018). Dtaidistance. https://dtaidistance.readthedocs.io/ en/latest/ .

[Moody and Mark, 2001] Moody, G. and Mark, R. (2001). The impact of the mit-bih arrhythmia database. IEEE engineering in medicine and biology magazine : the quarterly magazine of the Engineering in Medicine &amp; Biology Society , 20:45-50.

[Novikov, 2019] Novikov, A. (2019). Pyclustering: data mining library. Journal of Open Source Software , 4(36):1230.

[Özbay et al., 2006] Özbay, Y., Ceylan, R., and Karlik, B. (2006). A fuzzy clustering neural network architecture for classification of ecg arrhythmias. Computers in Biology and Medicine , 36(4):376-388.

[Paparrizos and Gravano, 2015] Paparrizos, J. and Gravano, L. (2015). k-shape: Efficient and accurate clustering of time series. In Proceedings of the 2015 ACM SIGMOD International Conference on Management of Data , pages 1855-1870. ACM.

[Pedregosa et al., 2011] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al. (2011). Scikit-learn: Machine learning in python. Journal of machine learning research , 12(Oct):2825-2830.

[Rosenberg and Hirschberg, 2007] Rosenberg, A. and Hirschberg, J. (2007). V-measure: A conditional entropy-based external cluster evaluation measure. In Proceedings of the 2007 joint conference on empirical methods in natural language processing and computational natural language learning (EMNLPCoNLL) .

[Roweis and Ghahramani, 1999] Roweis, S. and Ghahramani, Z. (1999). A unifying review of linear gaussian models. Neural computation , 11(2):305-345.

[Seabold and Perktold, 2010] Seabold, S. and Perktold, J. (2010). Statsmodels: Econometric and statistical modeling with python. In Proceedings of the 9th Python in Science Conference , volume 57, page 61. Scipy.

[Simchowitz et al., 2018] Simchowitz, M., Mania, H., Tu, S., Jordan, M. I., and Recht, B. (2018). Learning without mixing: Towards a sharp analysis of linear system identification. arXiv preprint arXiv:1802.08334 .

[Smyth, 1997] Smyth, P. (1997). Clustering sequences with hidden markov models. In Advances in neural information processing systems , pages 648-654.

[Stoica et al., 1988] Stoica, P., Friedlander, B., and Söderström, T. (1988). A high-order yule-walker method for estimation of the ar parameters of an arma model. Systems &amp; control letters , 11(2):99-105.

[Stoica et al., 2005] Stoica, P., Moses, R. L., et al. (2005). Spectral analysis of signals . Pearson Prentice Hall Upper Saddle River, NJ.

[Tavenard, 2017] Tavenard, R. (2017). tslearn: A machine learning toolkit dedicated to time-series data. URL https://github. com/rtavenar/tslearn .

[Tiao and Tsay, 1983] Tiao, G. C. and Tsay, R. S. (1983). Consistency properties of least squares estimates of autoregressive parameters in arma models. The Annals of Statistics , pages 856-871.

[Tsay and Tiao, 1984] Tsay, R. S. and Tiao, G. C. (1984). Consistent estimates of autoregressive pa- rameters and extended sample autocorrelation function for stationary and nonstationary arma models. Journal of the American Statistical Association , 79(385):84-96.

[Tsiamis and Pappas, 2019] Tsiamis, A. and Pappas, G. J. (2019). Finite sample analysis of stochastic system identification. arXiv preprint arXiv:1903.09122 .

[Turakhia et al., 2018] Turakhia, M., Desai, M., Hedlin, H., Rajmane, A., Talati, N., Ferris, T., Desai, S., Nag, D., Patel, M., Kowey, P., Rumsfeld, J., M. Russo, A., True Hills, M., B. Granger, C., W. Mahaffey, K., and Perez, M. (2018). Rationale and design of a large-scale, app-based study to identify cardiac arrhythmias using a smartwatch: The apple heart study. American Heart Journal , 207.

[Vinh et al., 2010] Vinh, N. X., Epps, J., and Bailey, J. (2010). Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance. Journal of Machine Learning Research , 11(Oct):2837-2854.

[Vishwanathan et al., 2007] Vishwanathan, S., Smola, A. J., and Vidal, R. (2007). Binet-cauchy kernels on dynamical systems and its application to the analysis of dynamic scenes. International Journal of Computer Vision , 73(1):95-119.

[Yeh et al., 2012] Yeh, Y.-C., Chiou, C. W., and Lin, H.-J. (2012). Analyzing ecg for cardiac arrhythmia using cluster analysis. Expert Systems with Applications , 39(1):1000-1010.

## A Proofs for model equivalence

In this section, we prove a generalization of Theorem 4.1 for both LDSs with observed inputs and LDSs with hidden inputs.

## A.1 Preliminaries

Sum of ARMA processes It is known that the sum of ARMA processes is still an ARMA process.

Lemma A.1 (Main Theorem in [Granger and Morris, 1976]) . The sum of two independent stationary series generated by ARMA( p , m ) and ARMA( , q n ) is generated by ARMA( x , y ), where x ≤ p + q and y ≤ max( p + n, q + m ) .

In shorthand notation, ARMA ( p, m ) + ARMA ( q, n ) = ARMA ( p + q, max( p + n, q + m )) .

When two ARMAX processes share the same exogenous input series, the dependency on exogenous input is additive, and the above can be extended to ARMAX ( p, m, r ) + ARMAX ( q, n, s ) = ARMAX ( p + q, max( p + n, q + m , ) max( r, s )) .

glyph[negationslash]

Jordan canonical form and canonical basis Every square real matrix is similar to a complex block diagonal matrix known as its Jordan canonical form (JCF). In the special case for diagonalizable matrices, JCF is the same as the diagonal form. Based on JCF, there exists a canonical basis { e i } consisting only of eigenvectors and generalized eigenvectors of A . A vector v is a generalized eigenvector of rank µ with corresponding eigenvalue λ if ( λI -A ) µ v = 0 and ( λI -A ) µ -1 v = 0 .

Relating the canonical basis to the characteristic polynomial, the characteristic polynomial can be completely factored into linear factors χ A ( λ ) = ( λ -λ 1 ) µ 1 ( λ -λ 2 ) µ 2 · · · ( λ -λ r ) µ r over C . The complex roots λ , 1 · · · , λ r are eigenvalues of A . For each eigenvalue λ i , there exist µ i linearly independent generalized eigenvectors v such that ( λ I i -A ) µ i v = 0 .

## A.2 General model equivalence theorem

Now we state Theorem A.1, a more detailed version of Theorem 4.1.

Theorem A.1. For any linear dynamical system with parameters Θ = ( A,B,C,D ) , hidden dimension n , inputs x t ∈ R k , and outputs y t ∈ R m , the outputs y t satisfy

$$\nu ^ { \cdot \nu } \chi ^ { \dagger } _ { A } ( L ) y _ { t } = \chi ^ { \dagger } _ { A } ( L ) \xi _ { t } + \Gamma ( L ) x _ { t }, \quad \ \ ( 5 ) \quad \text{app}.$$

where L is the lag operator, χ † A ( L ) = L χ n A ( L -1 ) is the reciprocal polynomial of the characteristic polynomial of A , and Γ( L ) is an m -byk matrix of polynomials of degree n -1 .

This implies that each dimension of y t can be generated by an ARMAX( n, n, n -1 ) model, where the autoregressive parameters are the characteristic polynomial coefficients in reverse order and in negative values.

To prove the theorem, we introduce a lemma to analyze the autoregressive behavior of the hidden state projected to a generalized eigenvector direction.

Lemma A.2. Consider a linear dynamical system with parameters Θ = ( A,B,C,D ) , hidden states h t ∈ R n , inputs x t ∈ R k , and outputs y t ∈ R m as defined in (1). For any generalized eigenvector e i of A ∗ with eigenvector λ and rank µ , the lag operator polynomial (1 -λL ) µ applied to time series h ( ) i t := 〈 h , e t i 〉 results in

$$\cdot, n \right ) = \quad ( 1 { - \lambda } L ) ^ { \mu } h _ { t } ^ { ( i ) } = l i n e a r \ t r a n s f o r m a t i o n \ o f \ x _ { t }, \cdots, x _ { t - \mu + 1 }.$$

Proof. To expand the LHS, first observe that

$$\begin{array} { l l l } \cdot \\, & ( 1 - \lambda L ) h _ { t } ^ { ( i ) } & = ( 1 - \lambda L ) \langle h _ { t }, e _ { i } \rangle \\ & = \langle h _ { t } ^ { ( i ) }, e _ { i } \rangle - \lambda L \langle h _ { t } ^ { ( i ) }, e _ { i } \rangle \\ & = \langle A h _ { t - 1 } + B x _ { t }, e _ { i } \rangle - \langle h _ { t - 1 }, \lambda e _ { i } \rangle \\ \cdot \\ \vdots & \end{array} \cdot \\ \quad \begin{array} { l l l l } ( 1 - \lambda L ) h _ { t } ^ { ( i ) }, e _ { i } \rangle \\ ( A ^ { * } - \lambda I ) e _ { i } \rangle + \langle B x _ { t }, e _ { i } \rangle. \\ \end{array} \\ \quad \begin{array} { l l l } \cdot \\ \end{array} \\ \quad \begin{array} { l l l } \cdot \\ \end{array} \\ \cdot \end{array}$$

We can apply (1 -λL ) again similarly to obtain

$$( 1 - \lambda L ) ^ { 2 } h _ { t } ^ { ( i ) } = \langle h _ { t - 2 }, ( A ^ { * } - \lambda I ) ^ { 2 } e _ { i } \rangle \\ + \langle B x _ { t - 1 }, ( A ^ { * } - \lambda I ) e _ { i } \rangle + ( 1 - \lambda L ) \langle B x _ { t }, e _ { i } \rangle,$$

and in general we can show inductively that

$$& ( 1 - \lambda L ) ^ { k } h _ { t } ^ { ( i ) } - \langle h _ { t - k }, ( A ^ { * } - \lambda I ) ^ { k } e _ { i } \rangle = \\ & \sum _ { j = 0 } ^ { k - 1 } ( 1 - \lambda L ) ^ { k - 1 - j } L ^ { j } \langle B x _ { t }, ( A ^ { * } - \lambda I ) ^ { j } e _ { i } \rangle,$$

where the RHS is a linear transformation of x , t · · · , x t -k +1 .

Since ( λI -A ∗ ) µ e i = 0 by definition of generalized eigenvectors, 〈 h t -µ , ( A ∗ -λI ) µ e i 〉 = 0 , and hence (1 -λL ) µ h ( ) i t itself is a linear transformation of x , t · · · , x t -µ +1 .

Proof for Theorem A.1 Using Lemma A.2 and the canonical basis, we can prove Theorem A.1.

Proof. Let λ , 1 · · · , λ r be the eigenvalues of A with multiplicity µ , 1 · · · , µ r . Since A is a real-valued matrix, its adjoint A ∗ has the same characteristic polynomial and eigenvalues as A . There exists a canonical basis { e i } n i =1 for A ∗ , where e , 1 · · · , e µ 1 are generalized eigenvectors

with eigenvalue λ 1 , e µ 1 +1 , · · · , e µ 1 + µ 2 are generalized eigenvectors with eigenvalue λ 2 , so on and so forth, and e µ 1 + ··· + µ r -1 +1 , · · · , e µ 1 + ··· + µ r are generalized eigenvectors with eigenvalue λ r .

By Lemma (A.2), (1 -λ L 1 ) µ 1 h ( ) i t is a linear transformation of x , t · · · , x t -µ 1 +1 for i = 1 , · · · , µ 1 ; (1 -λ L 2 ) µ 2 h ( ) i t is a linear transformation of x , t · · · , x t -µ 2 +1 for i = µ 1 +1 , · · · , µ 1 + µ 2 ; so on and so forth; (1 -λ L r ) µ r h ( ) i t is a linear transformation of x , t · · · , x t -µ r +1 for i = µ 1 + · · · + µ r -1 +1 , · · · , n .

glyph[negationslash]

glyph[negationslash]

We then apply lag operator polynomial Π j = i (1 -λ L j ) µ j to both sides of each equation. The lag polynomial in the LHS becomes (1 -λ L 1 ) µ 1 · · · (1 -λ L r ) µ r = χ † A ( L . ) For the RHS, since Π j = i (1 -λ L j ) µ j is of degree n -µ i , it lags the RHS by at most n -µ i additional steps, and the RHS becomes a linear transformation of x , t · · · , x t -n +1 .

Thus, for each i , χ † A ( L h ) ( ) i t is a linear transformation of x , t · · · , x t -n +1 .

The outputs of the LDS are defined as y t = Ch t + Dx t + ξ t = ∑ n i =1 h ( ) i t Ce i + Dx t + ξ t . By linearity, and since χ † A ( L ) is of degree n , both ∑ n i =1 h ( ) i t Ce i and χ † A ( L Dx ) t are linear transformations of x , t · · · , x t -n . We can write any such linear transformation as Γ( L x ) t for some m -byk matrix Γ( L ) of polynomials of degree n -1 . Thus, as desired,

$$\chi _ { A } ^ { \dagger } ( L ) y _ { t } = & \chi _ { A } ^ { \dagger } ( L ) \xi _ { t } + \Gamma ( L ) x _ { t }.$$

Assuming that there are no common factors in χ † A and Γ , χ † A is then the lag operator polynomial that represents the autoregressive part of y t . This assumption is the same as saying that y t cannot be expressed as a lower-order ARMA process. The reciprocal polynomial has the same coefficients in reverse order as the original polynomial. According to the lag operator polynomial on the LHS, 1 -ϕ L 1 -ϕ L 2 2 -··· -ϕ L n n = χ † A ( L ) , and L n -ϕ L 1 n -1 -···-ϕ n = χ A ( L ) , so the i -th order autoregressive parameter ϕ i is the negative value of the ( n -i ) -th order coefficient in the characteristic polynomial χ A .

## A.3 The hidden input case as a corollary

The statement about LDS without external inputs in Theorem 4.1 comes as a corollary to Theorem A.1, with a short proof here.

Proof. Define y ′ t = Ch t + Dx t to be the output without noise, i.e. y t = y ′ t + ξ t . By Theorem A.1, χ † A ( L y ) ′ t = Γ( L x ) t . Since we assume the hidden inputs x t are i.i.d.

Gaussians, y ′ t is then generated by an ARMA( n, n -1 ) process with autoregressive polynomial χ † A ( L ) .

The output noise ξ t itself can be seen as an ARMA( 0 0 , ) process. By Lemma A.1, ARMA( n, n -1 ) + ARMA( 0 0 , ) = ARMA( n +0 max( , n +0 , n -1 + 0) ) = ARMA( n, n ). Hence the outputs y t are generated by an ARMA( n, n ) process as claimed in Theorem 4.1. It is easy to see in the proof of Lemma A.1 that the autoregressive parameters do not change when adding a white noise [Granger and Morris, 1976].

## B Proof for eigenvalue approximation theorems

Here we restate Theorem 4.2 and Theorem 4.3 together, and prove it in three steps for 1) the general case, 2) the simple eigenvalue case, and 3) the explicit condition number bounds for the simple eigenvalue case.

Theorem B.1. Suppose y t are the outputs from an n -dimensional latent linear dynamical system with parameters Θ = ( A,B,C,D ) and eigenvalues λ , 1 · · · , λ n . Let ˆ Φ = (ˆ ϕ , 1 · · · , ϕ ˆ n ) be the estimated autoregressive parameters with error ‖ ˆ Φ -Φ ‖ = glyph[epsilon1] , and let r , 1 · · · , r n be the roots of the polynomial 1 -ˆ ϕ z 1 -··· -ˆ ϕ z n n .

Assuming the LDS is observable, the roots converge to the true eigenvalues with convergence rate O ( glyph[epsilon1] 1 /n ) . If all eigenvalues of A are simple (i.e. multiplicity 1), then the convergence rate is O ( glyph[epsilon1] ) . If A is symmetric, Lyapunov stable (spectral radius at most 1), and only has simple eigenvalues, then glyph[negationslash]

$$\cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad\cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \ \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cd. \quad \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cd. \quad \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad\cd. \quad \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot. \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cd� \quad \quad \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cd� \quad \quad \cdot \quad \cdot \quad \cd� \quad \cdot \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd. \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd. \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd1 \quad. \quad. \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \d� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \d. \quad. \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad. \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad. \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad. \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \cd� \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad. \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd� \quad \cd�$$

## B.1 General (1/n)-exponent bound

This is a known perturbation bound on polynomial root finding due to Ostrowski [Beauzamy, 1999].

Lemma B.1. Let Φ( ) = z z n + ϕ z 1 n -1 + · · · + ϕ n -1 z + ϕ n and Ψ( ) = z z n + ψ z 1 n -1 + · · · + ψ n -1 z + ψ n be two polynomials of degree n . If ‖ Φ -Ψ ‖ 2 &lt; glyph[epsilon1] , then the roots ( r k ) of Φ and roots (˜ ) r k of Ψ under suitable order satisfy

$$| r _ { k } - \tilde { r } _ { k } | \leq 4 C p \epsilon ^ { 1 / n },$$

where C = max 1 0 , ≤ ≤ k n {| ϕ n | 1 /n , | ψ n | 1 /n } .

The general O ( glyph[epsilon1] 1 /n ) convergence rate in Theorem 4.2 follows directly from Lemma B.1 and Theorem 4.1.

## B.2 Bound for simple eigenvalues

The 1 n -exponent in the above bound might seem not very ideal, but without additional assumptions the

1 n -exponent is tight. As an example, the polynomial x 2 -glyph[epsilon1] has roots x ± √ glyph[epsilon1] . This is a general phenomenon that a root with multiplicity m could split into m roots at rate O glyph[epsilon1] ( m ) , and is related to the regular splitting property [Hryniv and Lancaster, 1999, Lancaster et al., 2003] in matrix eigenvalue perturbation theory.

Under the additional assumption that all the eigenvalues are simple (no multiplicity), we can prove a better bound using the following idea with companion matrix: Small perturbation in autoregressive parameters results in small perturbation in companion matrix, and small perturbation in companion matrix results in small perturbation in eigenvalues.

Matrix eigenvalue perturbation theory The perturbation bound on eigenvalues is a well-studied problem [Greenbaum et al., 2019]. The regular splitting property states that, for an eigenvalue λ 0 with partial multiplicities m , 1 · · · , m k , an O glyph[epsilon1] ( ) perturbation to the matrix could split the eigenvalue into M = m 1 + · · · + m k distinct eigenvalues λ ij ( glyph[epsilon1] ) for i = 1 , · · · , k and j = 1 , · · · , m i , and each eigenvalue λ ij ( glyph[epsilon1] ) is moved from the original position by O glyph[epsilon1] ( 1 /m i ) .

For semi-simple eigenvalues, geometric multiplicity equals algebraic multiplicity. Since geometric multiplicity is the number of partial multiplicities while algebraic multiplicity is the sum of partial multiplicities, for semi-simple eigenvalues all partial multiplicities m i = 1 . Therefore, the regular splitting property corresponds to the asymptotic relation in equation 6. It is known that regular splitting holds for any semi-simple eigenvalue even for non-Hermitian matrices.

Lemma B.2 (Theorem 6 in [Lancaster et al., 2003]) . Let L λ, glyph[epsilon1] ( ) be an analytic matrix function with semisimple eigenvalue λ 0 at glyph[epsilon1] = 0 of multiplicity M . Then there are exactly M eigenvalues λ i ( glyph[epsilon1] ) of L λ, glyph[epsilon1] ( ) for which λ i ( glyph[epsilon1] ) → λ 0 as glyph[epsilon1] → 0 , and for these eigenvalues

$$\lambda _ { i } ( \epsilon ) = \lambda _ { 0 } + \lambda _ { i } ^ { \prime } \epsilon + o ( \epsilon ). \quad \quad ( 6 ) \quad \text{ori}$$

Companion Matrix Matrix perturbation theory tell us how perturbations on matrices change eigenvalues, while we are interested in how perturbations on polynomial coefficients change roots. To apply matrix perturbation theory on polynomials, we introduce the companion matrix , also known as the controllable canonical form in control theory.

Definition B.1. For a monic polynomial Φ( u ) = z n + ϕ z 1 n -1 + · · · + ϕ n -1 z + ϕ n , the companion matrix of the polynomial is the square matrix

$$C ( \Phi ) = \begin{bmatrix} 0 & 0 & \dots & 0 & - \varphi _ { n } \\ 1 & 0 & \dots & 0 & - \varphi _ { n - 1 } \\ 0 & 1 & \dots & 0 & - \varphi _ { n - 2 } \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & - \varphi _ { 1 } \end{bmatrix}.$$

The matrix C (Φ) is the companion in the sense that its characteristic polynomial is equal to Φ .

In relation to a pure autoregressive AR( p ) model, the companion matrix corresponds to the transition matrix in the linear dynamical system when we encode the values form the past p lags as a p -dimensional state

$$h _ { t } = \left [ y _ { t - p + 1 } \ \cdots$$

If y t = ϕ y 1 t -1 + · · · + ϕ y p t -p , then h t =

$$& \ s p l i t _ { - } \quad \mathbb { m } { y }$$

## Proof of Theorem 4.2 for simple eigenvalues

Proof. Let y t be the outputs of a linear dynamical system S with only simple eigenvalues, and let Φ = ( ϕ , 1 · · · , ϕ n ) be the ARMAX autoregressive parameters for y t . Let C (Φ) be the companion matrix of the polynomial z n -ϕ z 1 n -1 -ϕ z 2 n -2 - · · · -ϕ n . The companion matrix is the transition matrix of the LDS described in equation 7. Since this LDS the same autoregressive parameters and hidden state dimension as the original LDS, by Corollary 4.1 the companion matrix has the same characteristic polynomial as the original LDS, and thus also has simple (and hence also semi-simple) eigenvalues. The O glyph[epsilon1] ( ) convergence rate then follows from Lemma B.2 and Theorem 5.1, as the error on ARMAX parameter estimation can be seen as perturbation on the companion matrix.

A note on the companion matrix One might hope that we could have a more generalized result using Lemma B.2 for all systems with semi-simple eigenvalues instead of restricting to matrices with simple eigenvalues. Unfortunately, even if the original linear dynamical system has only semi-simple eigenvalues, in general the companion matrix is not semi-simple unless the original linear dynamical system is simple. This is because the

companion matrix always has its minimal polynomial equal to its characteristic polynomial, and hence has geometric multiplicity 1 for all eigenvalues. This also points to the fact that even though the companion matrix has the form of the controllable canonical form, in general it is not necessarily similar to the transition matrix in the original LDS.

## B.3 Explicit bound for condition number

In this subsection, we write out explicitly the condition number for simple eigenvalues in the asymptotic relation λ glyph[epsilon1] ( ) = λ 0 + κglyph[epsilon1] + ( ) o glyph[epsilon1] , to show how it varies according to the spectrum. Here we use the notation κ C, λ ( ) to note the condition number for eigenvalue λ in companion matrix C .

Lemma B.3. For a companion matrix C with simple eigenvalues λ , 1 · · · , λ n , the eigenvalues λ , ′ 1 · · · , λ ′ n of the perturbed matrix by C + δC satisfy

$$| \lambda _ { j } - \lambda _ { j } ^ { \prime } | \le$$

and the condition number κ C, λ ( j ) is bounded by glyph[negationslash]

glyph[negationslash]

$$\begin{array} {$$

where ρ C ( ) is the spectral radius, i.e. largest absolute value of its eigenvalues.

In particular, when ρ C ( ) ≤ 1 , i.e. when the matrix is Lyapunov stable, glyph[negationslash]

$$& \lambda _ { y \lambda _ { p } } ^ { y \lambda _ { p }$$

Proof. For each simple eigenvalue λ of the companion matrix C with column eigenvector v and row eigenvector w ∗ , the condition number of the eigenvalue is

$$\kappa ( C, \lambda ) = \frac { \| w \| _ { 2 } \| v \| _ { 2 } } { | w ^ { * } v | }. \quad \quad ( 1 1 ) \quad \ \overline { L }.$$

This is derived from differentiating the eigenvalue equation Cv = vλ , and multiplying the differentiated equation by w ∗ , which results in

$$w ^ { * } ( \delta C ) v + w ^ { * } C ( \delta v ) = \lambda w ^ { * } ( \delta v ) + w ^ { * } v ( \delta \lambda ). \\ \delta \lambda = \frac { w ^ { * } ( \delta C ) v } { w ^ { * } v }.$$

$$\delta \lambda = \frac { w ^ { * } ( \delta C ) v } { w ^ { * } v }.$$

Therefore,

$$| \delta \lambda | \leq \frac { \| w \| _ { 2 } \| v \| _ { 2 } } { | w ^ { * } v | } \| \delta C \| _ { 2 } = \kappa ( C, \lambda ) \| \delta C \| _ { 2 }. \quad ( 1 2 )$$

The companion matrix can be diagonalized as C = V -1 diag ( λ , 1 · · · , λ n ) V , the rows of the Vandermonde matrix V are the row eigenvectors of C , while the columns of V -1 are the column eigenvectors of C . Since the the j -th row V j, ∗ and the j -th column V -1 ∗ ,j have inner product 1 by definition of matrix inverse, the condition number is given by

$$\kappa ( C, \lambda _ { j } ) = \| V _ { j, * } \| _ { 2 } \ \| V _ { *, j } ^ { - 1 } \| _ { 2 }. \quad \ \ ( 1 3 )$$

Formula for inverse of Vandermonde matrix The Vandermonde matrix is defined as

$$V = \begin{bmatrix} 1 & \lambda _ { 1 } & \lambda _ { 1 } ^ { 2 } & \cdots & \lambda _ { 1 } ^ { p - 1 } \\ 1 & \lambda _ { 2 } & \lambda _ { 2 } ^ { 2 } & \cdots & \lambda _ { 2 } ^ { p - 1 } \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ 1 & \lambda _ { p } & \lambda _ { p } ^ { 2 } & \cdots & \lambda _ { p } ^ { p - 1 } \end{bmatrix}. \quad ( 1 4 )$$

The inverse of the Vandermonde matrix V is given by [El-Mikkawy, 2003] using elementary symmetric polynomial.

$$( V ^ { - 1 } ) _ { i, j } = \frac { ( - 1 ) ^ { i + j } S _ { p - i, j } } { \prod _ { k < j } ( \lambda _ { j } - \lambda _ { k } ) \prod _ { k > j } ( \lambda _ { k } - \lambda _ { j } ) }, \quad ( 1 5 )$$

where S p -i,j = S p -i ( λ , 1 · · · , λ j -1 , λ j +1 , · · · , λ p ) .

Pulling out the common denominator, the j -th column vector of V -1 is

$$\frac { ( - 1 ) ^ { j } } { \prod _ { k < j } ( \lambda _ { j } - \lambda _ { k } ) \prod _ { k > j } ( \lambda _ { k } - \lambda _ { j } ) } \begin{bmatrix} ( - 1 ) S _ { p - 1 } \\ ( - 1 ) ^ { 2 } S _ { p - 2 } \\ \vdots \\ ( - 1 ) ^ { p - 1 } S _ { 1 } \\ ( - 1 ) ^ { p } \end{bmatrix},$$

where the elementary symmetric polynomials are over variables λ , 1 · · · , λ j -1 , λ j +1 , · · · , λ p .

For example, if p = 4 , then the 3rd column (up to scaling) would be

$$\begin{matrix} ( 1 1 ) & & - 1 \\ & & \frac { - 1 } { ( \lambda _ { 3 } - \lambda _ { 1 } ) ( \lambda _ { 3 } - \lambda _ { 2 } ) ( \lambda _ { 4 } - \lambda _ { 3 } ) } \begin{bmatrix} - \lambda _ { 1 } \lambda _ { 2 } \lambda _ { 4 } \\ \lambda _ { 1 } \lambda _ { 2 } + \lambda _ { 1 } \lambda _ { 4 } + \lambda _ { 2 } \lambda _ { 4 } \\ - \lambda _ { 1 } - \lambda _ { 2 } - \lambda _ { 4 } \\ 1 \end{bmatrix}. \end{matrix}$$

Bounding the condition number As discussed before, the condition number for eigenvalue λ j is

$$\kappa ( C, \lambda _ { j } ) = \| V _ { j, * } \| _ { 2 } \ \| V _ { *, j } ^ { - 1 } \| _ { 2 }.$$

where V j, ∗ is the j -th row of the Vandermonde matrix V and V -1 ∗ ,j is the j -th column of V -1 .

$$\text{By definition } V _ { j, * } = \left [ 1 \ \lambda _ { j } \ \lambda _ { j } ^ { 2 } \ \cdots \ \lambda _ { j } ^ { p - 1 } \right ], \, \text{so} \\ \| V _ { j, * } \| _ { 2 } = \left ( \sum _ { i = 0 } ^ { p - 1 } \lambda _ { j } ^ { 2 i } \right ) ^ { 1 / 2 }.$$

Using the above explicit expression for V -1 , ‖ V -1 ∗ ,j ‖ 2 =

glyph[negationslash]

$$\frac { 1 } { \prod _ { k \neq j } | \lambda _ { j } - \lambda _ { k } | } \left ( \sum _ { i = 0 } ^ { p - 1 } S _ { i } ^ { 2 } ( \lambda _ { 1 }, \cdots, \lambda _ { j - 1 }, \lambda _ { j + 1 }, \cdots, \lambda _ { p } ) \right ) ^ { 1 / 2 } \, \dot { \ } d e t a r e \, \text{expl} \\ \cdot$$

Therefore, glyph[negationslash]

$$\Gamma \text{mere,} & & \text{or use} \\ \kappa ( C, \lambda _ { j } ) = \frac { 1 } { \prod _ { k \neq j } | \lambda _ { j } - \lambda _ { k } | } \left ( \sum _ { i = 0 } ^ { p - 1 } \lambda _ { j } ^ { 2 i } \right ) ^ { 1 / 2 } & & \text{Theo} \\ & & \left ( \sum _ { i = 0 } ^ { p - 1 } S _ { i } ^ { 2 } ( \lambda _ { 1 }, \cdots, \lambda _ { j - 1 }, \lambda _ { j + 1 }, \cdots, \lambda _ { p } ) \right ) ^ { 1 / 2 }. } & & \text{multi} \\ & & \text{size} \, \text{p} \\ & & \text{then}.$$

Note that both parts under ( · · · ) 1 / 2 are greater than or equal to 1 , so we can bound it below by glyph[negationslash]

$$\kappa ( C, \lambda _ { j } ) \geq \frac { 1 } { \prod _ { k \neq j } | \lambda _ { j } - \lambda _ { k } | }.$$

We could also bound the two parts above. The first part can be bounded by

$$\left ( \sum _ { i = 0 } ^ { p - 1 } \lambda _ { j } ^ { 2 i } \right ) ^ { 1 / 2 } \leq \sqrt { p } \max ( 1, | \lambda _ { j } | ) ^ { ( p - 1 ) }. \quad ( 1 7 ) \quad \ \Big | \ \ 1$$

While for the second part, since

$$| S _ { i } ( \lambda _ { 1 }, \cdots, \lambda _ { j - 1 }, \lambda _ { j + 1 }, \cdots, \lambda _ { p } ) | \leq \binom { p - 1 } { i } | \lambda | _ { \max } ^ { i }, \quad \Big |$$

we have that

$$\sum _ { i = 0 } ^ { p - 1 } S _ { i } ^ { 2 } ( \lambda _ { 1 }, \cdots, \lambda _ { j - 1 }, \lambda _ { j + 1 }, \cdots, \lambda _ { p } ) & & \quad \Big | \ e \\ \leq \sum _ { i = 0 } ^ { p - 1 } \binom { p - 1 } { i } | \lambda | _ { \max } ^ { 2 i } = ( 1 + | \lambda | _ { \max } ^ { 2 } ) ^ { p - 1 }. & & \quad \text{end} \\ \mathcal { C } \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \end{array}$$

Combining equation 17 and 18 for the upper bound, and putting it together with the lower bound, glyph[negationslash]

glyph[negationslash]

$$\frac { 1 } { \prod _ { k \neq j } | \lambda _ { j } - \lambda _ { k } | } & \leq \kappa ( C, \lambda _ { j } ) \leq \int \, \text{word} \\ \frac { \sqrt { p } } { \prod _ { k \neq j } | \lambda _ { j } - \lambda _ { k } | } \, ( \max ( 1, | \lambda _ { j } | ) ) ^ { p - 1 } \, ( 1 + \rho ( C ) ^ { 2 } ) ^ { \frac { p - 1 } { 2 } }, \\ \cdot \, \cdot \, \cdot \, & \quad -$$

as desired.

.

Theorem 4.3 follows from Lemma B.3, because the estimation error on the autoregressive parameters can be seen as the perturbation on the companion matrix, and the companion matrix has the same eigenvalues as the original LDS.

## C Iterated regression for ARMAX

Algorithm We generalize Algorithm 1 to accommodate for exogenous inputs. Since the exogenous inputs are explicitly observed, including exogenous inputs in the regression does not change the consistent property of the estimator.

Theorem A.1 shows that different output channels from the same LDS have the same autoregressive parameters in ARMAX models. Therefore, we could leverage multidimensional outputs by estimating the autoregressive parameters in each channel separately and average them.

Algorithm 2: Regularized iterated regression for AR parameter estimation in ARMAX

Input: A time series { y t } T t =1 where y t ∈ R m , exogenous input series { x t } T t =1 where x t ∈ R k , and guessed hidden state dimension n .

for d = 1 , · · · , m do

Let y ( d ) t be the projection of y t to the d -th dimension;

Initialize error term estimates ˆ glyph[epsilon1] t = 0 glyph[vector] ∈ R m for t = 1 , · · · , T ;

for

i

= 0

,

· · ·

, n do

Perform glyph[lscript] 2 -regularized least squares regression on y t against lagged terms of y t , x t , and ˆ glyph[epsilon1] t to solve for coefficients ˆ ϕ j ∈ R , ˆ θ j ∈ R , and ˆ γ j ∈ R k in the linear equation y ( d ) t = c + ∑ n j =1 ˆ ϕ y j ( d ) t -j + ∑ n -1 j =1 ˆ γ x j t -j + ∑ i j =1 ˆ ˆ θ glyph[epsilon1] j t -j , with glyph[lscript] 2 -regularization only on ˆ θ j ; Update ˆ glyph[epsilon1] t to be the residuals from the most recent regression;

end

$$| \text{ Record } \hat { \Phi } ^ { ( d ) } = ( \hat { \varphi } _ { 1 }, \cdots, \hat { \varphi } _ { n } ) ;$$

end

Return the average estimate ˆ = Φ 1 d ( ˆ Φ (1) + · · · + ˆ Φ ( m ) ) .

Again as before the i -th iteration of the regression only uses error terms from the past i lags. In other words, the initial iteration is an ARMAX( n, 0 , n -1 ) regression, the first iteration is an ARMAX( n, 1 , n -1 ) regression, and so forth.

Time complexity The iterated regression in each dimension involves n + 1 steps of least squares regression each on at most n k ( + 2) variables. Therefore, the total time complexity of Algorithm 2 is O nm ( (( nk ) 2 T +( nk ) 3 )) = O mn k T ( 3 2 + mn k 4 3 ) , where T is the sequence length, n the hidden state dimension, m the output dimension, and k the input dimension.

## D Additional simulation details

## D.1 Synthetic data generation

First, we generate K cluster centers by generating LDSs with random matrices A,B,C of standard i.i.d. Gaussians. We assume that the output y t only depends on the hidden state h t but not the input x t , i.e. the matrix D is zero. When generating the random LDSs, we require that the spectral radius ρ A ( ) ≤ 1 , i.e. all eigenvalues of A have absolute values at most 1, and regenerate a new random matrix if the spectral radius is above 1. Our method also applies to the case of arbitrary spectral radius, this requirement is for the purpose of preventing numeric overflow in generated sequence. We also require that the glyph[lscript] 2 distance d (Θ 1 , Θ ) = 2 ‖ λ A ( 1 ) -λ A ( 2 ) ‖ 2 between cluster centers are at least 0.2 apart.

Then, we generate 100 LDSs by randomly assigning them to the clusters. To obtain a LDS with assigned cluster center Θ = ( A , B , C c c c ) , we generate A ′ by adding a i.i.d. Gaussians to each entry of A c , while B ′ and C ′ are new random matrices of i.i.d. standard Gaussians. The standard deviation of the i.i.d. Gaussians for A ′ -A c is chosen such that the average distance to cluster centers is less than half of the inter-cluster distance between centers.

For each LDS, we generate a sequence by drawing hidden inputs x t ∼ N (0 , 1) and put noise ξ t ∼ N (0 , 0 01 . 2 ) on the outputs.

## D.2 Empirical correlation between AR distance and LDS distance.

Theorem 4.2 shows that LDSs with similar AR parameters also have similar eigenvalues. The converse of Theorem 4.2 is also true: dynamical systems with small eigenvalue distance have small autoregressive parameter distance, which follows from perturbation bounds for characteristic polynomials [Ipsen and Rehman, 2008]. Figure 2 shows simulation results where the AR parameter distance and the LDS eigenvalue distance are highly correlated.

Figure 2: The eigenvalue glyph[lscript] 2 distance and the autoregressive parameter glyph[lscript] 2 distance for 100 random linear dynamical systems with eigenvalues drawn uniformly randomly from [ -1 1] , . The two distance measures are highly correlated.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01039v3-figure-d9b1ff70-9e91-401d-8f03-bbefae146b5b.jpg)
Chloe Ching-Yun Hsu Michaela Hardt Moritz Hardt
Linear Dynamics: Clustering without identification
Linear dynamical systems are a fundamental and powerful parametric model
class. However, identifying the parameters of a linear dynamical system is a
venerable task, permitting provably efficient solutions only in special cases.
This work shows that the eigenspectrum of unknown linear dynamics can be
identified without full system identification. We analyze a computationally
efficient and provably convergent algorithm to estimate the eigenvalues of the
state-transition matrix in a linear dynamical system.
  When applied to time series clustering, our algorithm can efficiently cluster
multi-dimensional time series with temporal offsets and varying lengths, under
the assumption that the time series are generated from linear dynamical
systems. Evaluating our algorithm on both synthetic data and real
electrocardiogram (ECG) signals, we see improvements in clustering quality over
existing baselines.