Linear Dynamics: Clustering without identification

TL;DR

This paper introduces a method to identify the eigenspectrum of linear dynamical systems without full system identification, enabling efficient clustering of time series with temporal offsets and varying lengths.

Contribution

The work presents a novel, computationally efficient algorithm for estimating eigenvalues of linear systems, improving clustering of time series generated from such systems.

Findings

Algorithm accurately estimates eigenvalues without full system identification.

Improves clustering quality on synthetic and real ECG data.

Handles temporal offsets and varying lengths in time series.

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 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.

Tables

Table 1. 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.

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

Table 2. 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.10 (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)

Equations

$$h_t=A{h}_{t-1}+B{x}_t+{\zeta }_t$$

$$y_t=C{h}_t+D{x}_t+{\xi }_t.$$

$$y_t=c+{\Sigma }_{i=1}^p{\phi }_i{y}_{t-i}+{ϵ}_t,$$

$$y_t=c+{ϵ}_t+{\Sigma }_{i=1}^q{\theta }_i{ϵ}_{t-i},$$

$$y_t=c+{ϵ}_t+{\Sigma }_{i=1}^p{\phi }_i{y}_{t-i}+{\Sigma }_{i=1}^q{\theta }_i{ϵ}_{t-i}.$$

$$y_t=c+{ϵ}_t+{\Sigma }_{i=1}^p{\phi }_i{y}_{t-i}+{\Sigma }_{i=1}^q{\theta }_i{ϵ}_{t-i}+{\Sigma }_{i=0}^r{\gamma }_i{x}_{t-i},$$

$$\Phi (L)y_t=c+{ϵ}_t,$$

$$y_t=c+\Psi (L){ϵ}_t.$$

$$\Phi (L)y_t=c+\Psi (L){ϵ}_t+\Gamma (L)x_t$$

$$C(\Phi )=\left[\begin{array}{ccccc}0& 0& \dots & 0& -{\phi }_n\\ 1& 0& \dots & 0& -{\phi }_{n-1}\\ 0& 1& \dots & 0& -{\phi }_{n-2}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \dots & 1& -{\phi }_1\end{array}\right].$$

$$h_t={\left[\begin{array}{cccc}{y}_{t-p+1}& \cdots & y_{t-1}& y_t\end{array}\right]}^T.$$

$$\begin{array}{rl}& \left[\begin{array}{c}{y}_{t-p+1}\\ y_{t-p+2}\\ \cdots \\ y_{t-1}\\ y_t\end{array}\right]=\left[\begin{array}{ccccc}0& 1& 0& \dots & 0\\ 0& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & 1\\ {\phi }_p& {\phi }_{p-1}& {\phi }_{p-2}& \dots & {\phi }_1\end{array}\right]\left[\begin{array}{c}{y}_{t-p}\\ y_{t-p+1}\\ \cdots \\ y_{t-2}\\ y_{t-1}\end{array}\right]\\ & =C(-\Phi {)}^T{h}_{t-1}.\end{array}$$

$${\lambda }_i(ϵ)={\lambda }_0+{\lambda }_i^{\prime }ϵ+o(ϵ).$$

$$\frac{1}{{\prod }_{k\ne j}|{\lambda }_j-{\lambda }_k|}\le \kappa \le$$

$$\frac{\sqrt{n}}{{\prod }_{k\ne j}|{\lambda }_j-{\lambda }_k|}{\left(\max (1,|{\lambda }_j|)\right)}^{n-1}(1+\rho (A{)}^2{)}^{\frac{n-1}{2}},$$

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

$${\chi }_A^{†}(L)y_t={\chi }_A^{†}(L){\xi }_t+\Gamma (L)x_t,$$

$$(1-\lambda L{)}^{\mu }{h}_t^{(i)}=\text{linear transformation of }{x}_t,\cdots ,x_{t-\mu +1}.$$

$$(1-\lambda L)h_t^{(i)}$$

$$=⟨h_t^{(i)},e_i⟩-\lambda L⟨h_t^{(i)},e_i⟩$$

$$=⟨A{h}_{t-1}+B{x}_t,e_i⟩-⟨h_{t-1},\lambda e_i⟩$$

$$=⟨h_{t-1}^{(i)},(A^{\ast }-\lambda I)e_i⟩+⟨B{x}_t,e_i⟩.$$

$$(1-\lambda L{)}^2{h}_t^{(i)}=⟨h_{t-2},(A^{\ast }-\lambda I{)}^2{e}_i⟩$$

$$+⟨B{x}_{t-1},(A^{\ast }-\lambda I)e_i⟩+(1-\lambda L)⟨B{x}_t,e_i⟩,$$

$$(1-\lambda L{)}^k{h}_t^{(i)}-⟨h_{t-k},(A^{\ast }-\lambda I{)}^k{e}_i⟩=$$

$$\sum _{j=0}^{k-1}(1-\lambda L{)}^{k-1-j}{L}^j⟨B{x}_t,(A^{\ast }-\lambda I{)}^j{e}_i⟩,$$

$${\chi }_A^{†}(L)y_t=$$

$$|r_i-{\lambda }_i|\le \frac{\sqrt{n{2}^{n-1}}}{{\Pi }_{k\ne j}|{\lambda }_j-{\lambda }_k|}ϵ+\mathcal{O}({ϵ}^2).$$

$$|r_k-{\tilde{r}}_k|\le 4Cp{ϵ}^{1/n},$$

$${\lambda }_i(ϵ)={\lambda }_0+{\lambda }_i^{\prime }ϵ+o(ϵ).$$

$$C(\Phi )=\left[\begin{array}{ccccc}0& 0& \dots & 0& -{\phi }_n\\ 1& 0& \dots & 0& -{\phi }_{n-1}\\ 0& 1& \dots & 0& -{\phi }_{n-2}\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \dots & 1& -{\phi }_1\end{array}\right].$$

$$h_t={\left[\begin{array}{cccc}{y}_{t-p+1}& \cdots & y_{t-1}& y_t\end{array}\right]}^T.$$

$$\begin{array}{rl}& \left[\begin{array}{c}{y}_{t-p+1}\\ y_{t-p+2}\\ \cdots \\ y_{t-1}\\ y_t\end{array}\right]=\left[\begin{array}{ccccc}0& 1& 0& \dots & 0\\ 0& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & 1\\ {\phi }_p& {\phi }_{p-1}& {\phi }_{p-2}& \dots & {\phi }_1\end{array}\right]\left[\begin{array}{c}{y}_{t-p}\\ y_{t-p+1}\\ \cdots \\ y_{t-2}\\ y_{t-1}\end{array}\right]\\ & =C(-\Phi {)}^T{h}_{t-1}.\end{array}$$

$$|{\lambda }_j-{\lambda }_j^{\prime }|\le \kappa (C,{\lambda }_j)\Vert \delta C{\Vert }_2+o(\Vert \delta C{\Vert }_2^2),$$

Full text

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.

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 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 $\Theta=(A,B,C,D)$ receives inputs $x_{1},\cdots,x_{T}\in\mathbb{R}^{k}$, has hidden states $h_{0},\cdots,h_{T}\in\mathbb{R}^{n}$, and generates outputs $y_{1},\cdots,y_{T}\in\mathbb{R}^{m}$ according to the following time-invariant recursive equations:

[TABLE]

Assumptions. We assume that the stochastic noise $\zeta_{t}$ and $\xi_{t}$ are diagonal Gaussians. We also assume the system is observable, i.e. $C,CA,CA^{2},\cdots,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 $\mathbb{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 $\ell_{2}$ distance of the spectrum of the transition matrix $A$, i.e. $d(\Theta_{1},\Theta_{2})=\lVert\lambda(A_{1})-\lambda(A_{2})\rVert_{2}$, where $\lambda(A_{1})$ and $\lambda(A_{2})$ are the spectrum of $A_{1}$ and $A_{2}$ in sorted order. This distance definition satisfies non-negativity, 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.

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 $\mu$ with corresponding eigenvalue $\lambda$ if $(\lambda I-A)^{\mu}v=0$ and $(\lambda I-A)^{\mu-1}v\neq 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,

[TABLE]

where $\varphi_{1},\cdots,\varphi_{p}$ are autoregressive parameters, $c$ is a constant, and $\epsilon_{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,

[TABLE]

where $\theta_{1},\cdots,\theta_{q}$ are moving-average parameters, $c$ is a constant, and the errors $\epsilon_{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,

[TABLE]

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

[TABLE]

where $\{x_{t}\}$ is a known external time series, possibly multidimensional. In case $x_{t}$ is a vector, the parameters $\gamma_{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^{3}y_{t}=y_{t-3}$, and $(a_{2}L^{2}+a_{1}L+a_{0})y_{t}=a_{2}y_{t-2}+a_{1}y_{t-1}+a_{0}y_{t}$. The lag polynomials can be multiplied or inverted. An AR($p$) model can be characterized by

[TABLE]

where $\Phi(L)=1-\varphi_{1}L-\cdots-\varphi_{p}L^{p}$ is a polynomial of the lag operator $L$ of degree $p$. For example, any AR(2) model can be described as $(1-\varphi_{1}L-\varphi_{2}L^{2})y_{t}=c+\epsilon_{t}$.

Similarly, an MA($q$) can be characterized by a polynomial $\Psi(L)=\theta_{q}L^{q}+\cdots+\theta_{1}L+1$ of degree $q$,

[TABLE]

For example, for an MA(2) model the equation would be $y_{t}=c+(\theta_{2}L^{2}+\theta_{1}L+1)\epsilon_{t}$.

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

[TABLE]

where $\Phi,\Psi$, and $\Gamma$ are polynomials of degree $p,q$ and $r$. When the exogenous time series $x_{t}$ is multidimensional, $\Gamma(L)$ is a vector of degree-$r$ 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}\in\mathbb{R}^{m}$ be the outputs from a linear dynamical system with parameters $\Theta=(A,B,C,D)$, hidden dimension $n$, and inputs $x_{t}\in\mathbb{R}^{k}$. Each dimension of $y_{t}$ can be generated by an ARMAX($n,n,n-1$) model, whose autoregressive parameters $\varphi_{1},\cdots,\varphi_{n}$ can recover the characteristic polynomial of $A$ by $\chi_{A}(\lambda)=\lambda^{n}-\varphi_{1}\lambda^{n-1}-\cdots-\varphi_{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 $\lambda$ and rank $\mu$, the time series obtained from applying the lag operator polynomial $(1-\lambda L)^{\mu}$ to $\langle h_{t},e\rangle$ can be expressed as a linear combination of the past $k$ inputs $x_{t},\cdots,x_{t-k+1}$. Since $A$ is real-valued, $A$ and its adjoint $A^{*}$ share the same characteristic polynomial $\chi_{A}$.

We then consider the lag operator polynomial $\chi_{A}^{\dagger}(L)=L^{n}\chi_{A}(L^{-1})$, and show that the time series obtained from applying $\chi_{A}^{\dagger}(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 $\chi_{A}^{\dagger}$ 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 dynamical 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}<n_{2}$ can have the same autoregressive parameters if and only if $\chi_{A_{1}}(\lambda)=\chi_{A_{2}}(\lambda)\lambda^{n_{2}-n_{1}}$ and $\varphi_{n_{1}+1}=\cdots=\varphi_{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 over-parameterized 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 $\Theta=(A,B,C,D)$, eigenvalues $\lambda_{1},\cdots,\lambda_{n}$, and hidden inputs. Let $\hat{\Phi}=(\hat{\varphi}_{1},\cdots,\hat{\varphi}_{n})$ be the estimated autoregressive parameters for $\{y_{t}\}$ with error $\lVert\hat{\Phi}-\Phi\rVert=\epsilon$, and let $r_{1},\cdots,r_{n}$ be the roots of the polynomial $1-\hat{\varphi}_{1}z-\cdots-\hat{\varphi}_{n}z^{n}$.

Assuming the LDS is observable, the roots converge to the true eigenvalues with convergence rate $\mathcal{O}(\epsilon^{1/n})$. If all eigenvalues of $A$ are simple (no multiplicity), then the convergence rate is $\mathcal{O}(\epsilon)$.

Without additional assumptions, the $\frac{1}{n}$-exponent in the above general bound is tight. As an example, $z^{2}-\epsilon$ has roots $z\pm\sqrt{\epsilon}$. The general phenomenon that a root with multiplicity $m$ can split into $m$ roots at rate $O(\epsilon^{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(\epsilon)$ bound instead of $O(\epsilon^{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 $\Phi(u)=z^{n}+\varphi_{1}z^{n-1}+\cdots+\varphi_{n-1}z+\varphi_{n}$, the companion matrix, also known as the controllable canonical form in control theory, is the square matrix

[TABLE]

The matrix $C(\Phi)$ is the companion in the sense that it has $\Phi$ 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

[TABLE]

If $y_{t}=\varphi_{1}y_{t-1}+\cdots+\varphi_{p}y_{t-p}$, then $h_{t}=$

[TABLE]

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(\lambda,\epsilon)$ be an analytic matrix function with semi-simple eigenvalue $\lambda_{0}$ at $\epsilon=0$ of multiplicity $M$. Then there are exactly $M$ eigenvalues $\lambda_{i}(\epsilon)$ of $L(\lambda,\epsilon)$ for which $\lambda_{i}(\epsilon)\rightarrow\lambda_{0}$ as $\epsilon\rightarrow 0$, and for these eigenvalues

[TABLE]

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}-\lambda_{j}|\leq\kappa\epsilon+o(\epsilon^{2})$, then the condition number $\kappa$ is bounded by

[TABLE]

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

In particular, when $\rho(C)\leq 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}-\lambda_{j}|\leq\frac{\sqrt{n}(\sqrt{2})^{n-1}}{\prod_{k\neq j}|\lambda_{j}-\lambda_{k}|}\epsilon+o(\epsilon^{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($p$) models, 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 $\epsilon_{t}+\sum_{i=1}^{q}\theta_{i}\epsilon_{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.

Time complexity. The iterated regression involves $n+1$ steps of least squares regression each on at most $2n+1$ variables. Therefore, the total time complexity of Algorithm 1 is $O(n^{3}T+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 $\hat{\Phi}=(\hat{\varphi}_{1},\cdots,\hat{\varphi}_{n})$ from iterated regression converges in probability to the true parameters with rate

[TABLE]

or more explicitly, convergence in probability means that for all $\epsilon$, $\lim_{T\rightarrow\inf}\Pr(T^{1/2}|\hat{\Phi}-\Phi|>\epsilon)=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^{\prime},B^{\prime},C^{\prime},D^{\prime})$ are equivalent if $A^{\prime}=P^{-1}AP,B^{\prime}=P^{-1}B,C^{\prime}=CP,D^{\prime}=DP$, and $h_{t}^{\prime}=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 state-transition 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 multi-dimensional 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].

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^{3}T)$. 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-of-concept, 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 two-channel 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 $\ell_{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.

Appendix 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\leq p+q$ and $y\leq\max(p+n,q+m)$.

In shorthand notation, $\textrm{ARMA}(p,m)+\textrm{ARMA}(q,n)=\textrm{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 $\textrm{ARMAX}(p,m,r)+\textrm{ARMAX}(q,n,s)=\textrm{ARMAX}(p+q,\max(p+n,q+m),\max(r,s))$.

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 $\mu$ with corresponding eigenvalue $\lambda$ if $(\lambda I-A)^{\mu}v=0$ and $(\lambda I-A)^{\mu-1}v\neq 0$.

Relating the canonical basis to the characteristic polynomial, the characteristic polynomial can be completely factored into linear factors $\chi_{A}(\lambda)=(\lambda-\lambda_{1})^{\mu_{1}}(\lambda-\lambda_{2})^{\mu_{2}}\cdots(\lambda-\lambda_{r})^{\mu_{r}}$ over $\mathbb{C}$. The complex roots $\lambda_{1},\cdots,\lambda_{r}$ are eigenvalues of $A$. For each eigenvalue $\lambda_{i}$, there exist $\mu_{i}$ linearly independent generalized eigenvectors $v$ such that $(\lambda_{i}I-A)^{\mu_{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 $\Theta=(A,B,C,D)$, hidden dimension $n$, inputs $x_{t}\in\mathbb{R}^{k}$, and outputs $y_{t}\in\mathbb{R}^{m}$, the outputs $y_{t}$ satisfy

[TABLE]

where $L$ is the lag operator, $\chi_{A}^{\dagger}(L)=L^{n}\chi_{A}(L^{-1})$ is the reciprocal polynomial of the characteristic polynomial of $A$, and $\Gamma(L)$ is an $m$-by-$k$ 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 $\Theta=(A,B,C,D)$, hidden states $h_{t}\in\mathbb{R}^{n}$, inputs $x_{t}\in\mathbb{R}^{k}$, and outputs $y_{t}\in\mathbb{R}^{m}$ as defined in (1). For any generalized eigenvector $e_{i}$ of $A^{*}$ with eigenvector $\lambda$ and rank $\mu$, the lag operator polynomial $(1-\lambda L)^{\mu}$ applied to time series $h_{t}^{(i)}:=\langle h_{t},e_{i}\rangle$ results in

[TABLE]

Proof.

To expand the LHS, first observe that

[TABLE]

We can apply $(1-\lambda L)$ again similarly to obtain

[TABLE]

and in general we can show inductively that

[TABLE]

where the RHS is a linear transformation of $x_{t},\cdots,x_{t-k+1}$.

Since $(\lambda I-A^{*})^{\mu}e_{i}=0$ by definition of generalized eigenvectors, $\langle h_{t-\mu},(A^{*}-\lambda I)^{\mu}e_{i}\rangle=0$, and hence $(1-\lambda L)^{\mu}h_{t}^{(i)}$ itself is a linear transformation of $x_{t},\cdots,x_{t-\mu+1}$. ∎

Proof for Theorem A.1

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

Proof.

Let $\lambda_{1},\cdots,\lambda_{r}$ be the eigenvalues of $A$ with multiplicity $\mu_{1},\cdots,\mu_{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}\}_{i=1}^{n}$ for $A^{*}$, where $e_{1},\cdots,e_{\mu_{1}}$ are generalized eigenvectors with eigenvalue $\lambda_{1}$, $e_{\mu_{1}+1},\cdots,e_{\mu_{1}+\mu_{2}}$ are generalized eigenvectors with eigenvalue $\lambda_{2}$, so on and so forth, and $e_{\mu_{1}+\cdots+\mu_{r-1}+1},\cdots,e_{\mu_{1}+\cdots+\mu_{r}}$ are generalized eigenvectors with eigenvalue $\lambda_{r}$.

By Lemma (A.2), $(1-\lambda_{1}L)^{\mu_{1}}h_{t}^{(i)}$ is a linear transformation of $x_{t},\cdots,x_{t-\mu_{1}+1}$ for $i=1,\cdots,\mu_{1}$; $(1-\lambda_{2}L)^{\mu_{2}}h_{t}^{(i)}$ is a linear transformation of $x_{t},\cdots,x_{t-\mu_{2}+1}$ for $i=\mu_{1}+1,\cdots,\mu_{1}+\mu_{2}$; so on and so forth; $(1-\lambda_{r}L)^{\mu_{r}}h_{t}^{(i)}$ is a linear transformation of $x_{t},\cdots,x_{t-\mu_{r}+1}$ for $i=\mu_{1}+\cdots+\mu_{r-1}+1,\cdots,n$.

We then apply lag operator polynomial $\Pi_{j\neq i}(1-\lambda_{j}L)^{\mu_{j}}$ to both sides of each equation. The lag polynomial in the LHS becomes $(1-\lambda_{1}L)^{\mu_{1}}\cdots(1-\lambda_{r}L)^{\mu_{r}}=\chi_{A}^{\dagger}(L).$ For the RHS, since $\Pi_{j\neq i}(1-\lambda_{j}L)^{\mu_{j}}$ is of degree $n-\mu_{i}$, it lags the RHS by at most $n-\mu_{i}$ additional steps, and the RHS becomes a linear transformation of $x_{t},\cdots,x_{t-n+1}$.

Thus, for each $i$, $\chi_{A}^{\dagger}(L)h_{t}^{(i)}$ is a linear transformation of $x_{t},\cdots,x_{t-n+1}$.

The outputs of the LDS are defined as $y_{t}=Ch_{t}+Dx_{t}+\xi_{t}=\sum_{i=1}^{n}h_{t}^{(i)}Ce_{i}+Dx_{t}+\xi_{t}$. By linearity, and since $\chi_{A}^{\dagger}(L)$ is of degree $n$, both $\sum_{i=1}^{n}h_{t}^{(i)}Ce_{i}$ and $\chi_{A}^{\dagger}(L)Dx_{t}$ are linear transformations of $x_{t},\cdots,x_{t-n}$. We can write any such linear transformation as $\Gamma(L)x_{t}$ for some $m$-by-$k$ matrix $\Gamma(L)$ of polynomials of degree $n-1$. Thus, as desired,

[TABLE]

Assuming that there are no common factors in $\chi_{A}^{\dagger}$ and $\Gamma$, $\chi_{A}^{\dagger}$ 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-\varphi_{1}L-\varphi_{2}L^{2}-\cdots-\varphi_{n}L^{n}=\chi_{A}^{\dagger}(L)$, and $L^{n}-\varphi_{1}L^{n-1}-\cdots-\varphi_{n}=\chi_{A}(L)$, so the $i$-th order autoregressive parameter $\varphi_{i}$ is the negative value of the $(n-i)$-th order coefficient in the characteristic polynomial $\chi_{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}^{\prime}=Ch_{t}+Dx_{t}$ to be the output without noise, i.e. $y_{t}=y_{t}^{\prime}+\xi_{t}$. By Theorem A.1, $\chi_{A}^{\dagger}(L)y_{t}^{\prime}=\Gamma(L)x_{t}$. Since we assume the hidden inputs $x_{t}$ are i.i.d. Gaussians, $y_{t}^{\prime}$ is then generated by an ARMA($n,n-1$) process with autoregressive polynomial $\chi_{A}^{\dagger}(L)$.

The output noise $\xi_{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]. ∎

Appendix 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 $\Theta=(A,B,C,D)$ and eigenvalues $\lambda_{1},\cdots,\lambda_{n}$. Let $\hat{\Phi}=(\hat{\varphi}_{1},\cdots,\hat{\varphi}_{n})$ be the estimated autoregressive parameters with error $\lVert\hat{\Phi}-\Phi\rVert=\epsilon$, and let $r_{1},\cdots,r_{n}$ be the roots of the polynomial $1-\hat{\varphi}_{1}z-\cdots-\hat{\varphi}_{n}z^{n}$.

Assuming the LDS is observable, the roots converge to the true eigenvalues with convergence rate $\mathcal{O}(\epsilon^{1/n})$. If all eigenvalues of $A$ are simple (i.e. multiplicity 1), then the convergence rate is $\mathcal{O}(\epsilon)$. If $A$ is symmetric, Lyapunov stable (spectral radius at most 1), and only has simple eigenvalues, then

[TABLE]

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 $\Phi(z)=z^{n}+\varphi_{1}z^{n-1}+\cdots+\varphi_{n-1}z+\varphi_{n}$ and $\Psi(z)=z^{n}+\psi_{1}z^{n-1}+\cdots+\psi_{n-1}z+\psi_{n}$ be two polynomials of degree $n$. If $\lVert\Phi-\Psi\rVert_{2}<\epsilon$, then the roots $(r_{k})$ of $\Phi$ and roots $(\tilde{r}_{k})$ of $\Psi$ under suitable order satisfy

[TABLE]

where $C=\max_{1,0\leq k\leq n}\{|\varphi_{n}|^{1/n},|\psi_{n}|^{1/n}\}$.

The general $\mathcal{O}(\epsilon^{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 $\frac{1}{n}$-exponent in the above bound might seem not very ideal, but without additional assumptions the $\frac{1}{n}$-exponent is tight. As an example, the polynomial $x^{2}-\epsilon$ has roots $x\pm\sqrt{\epsilon}$. This is a general phenomenon that a root with multiplicity $m$ could split into $m$ roots at rate $O(\epsilon^{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 $\lambda_{0}$ with partial multiplicities $m_{1},\cdots,m_{k}$, an $O(\epsilon)$ perturbation to the matrix could split the eigenvalue into $M=m_{1}+\cdots+m_{k}$ distinct eigenvalues $\lambda_{ij}(\epsilon)$ for $i=1,\cdots,k$ and $j=1,\cdots,m_{i}$, and each eigenvalue $\lambda_{ij}(\epsilon)$ is moved from the original position by $O(\epsilon^{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(\lambda,\epsilon)$ be an analytic matrix function with semi-simple eigenvalue $\lambda_{0}$ at $\epsilon=0$ of multiplicity $M$. Then there are exactly $M$ eigenvalues $\lambda_{i}(\epsilon)$ of $L(\lambda,\epsilon)$ for which $\lambda_{i}(\epsilon)\rightarrow\lambda_{0}$ as $\epsilon\rightarrow 0$, and for these eigenvalues

[TABLE]

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 $\Phi(u)=z^{n}+\varphi_{1}z^{n-1}+\cdots+\varphi_{n-1}z+\varphi_{n}$, the companion matrix of the polynomial is the square matrix

[TABLE]

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

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

[TABLE]

If $y_{t}=\varphi_{1}y_{t-1}+\cdots+\varphi_{p}y_{t-p}$, then $h_{t}=$

[TABLE]

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 $\Phi=(\varphi_{1},\cdots,\varphi_{n})$ be the ARMAX autoregressive parameters for $y_{t}$. Let $C(\Phi)$ be the companion matrix of the polynomial $z^{n}-\varphi_{1}z^{n-1}-\varphi_{2}z^{n-2}-\cdots-\varphi_{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(\epsilon)$ 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 $\lambda(\epsilon)=\lambda_{0}+\kappa\epsilon+o(\epsilon)$, to show how it varies according to the spectrum. Here we use the notation $\kappa(C,\lambda)$ to note the condition number for eigenvalue $\lambda$ in companion matrix $C$.

Lemma B.3**.**

For a companion matrix $C$ with simple eigenvalues $\lambda_{1},\cdots,\lambda_{n}$, the eigenvalues $\lambda_{1}^{\prime},\cdots,\lambda_{n}^{\prime}$ of the perturbed matrix by $C+\delta C$ satisfy

[TABLE]

and the condition number $\kappa(C,\lambda_{j})$ is bounded by

[TABLE]

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

In particular, when $\rho(C)\leq 1$, i.e. when the matrix is Lyapunov stable,

[TABLE]

Proof.

For each simple eigenvalue $\lambda$ of the companion matrix $C$ with column eigenvector $v$ and row eigenvector $w^{*}$, the condition number of the eigenvalue is

[TABLE]

This is derived from differentiating the eigenvalue equation $Cv=v\lambda$, and multiplying the differentiated equation by $w*$, which results in

[TABLE]

[TABLE]

Therefore,

[TABLE]

The companion matrix can be diagonalized as $C=V^{-1}\text{diag}(\lambda_{1},\cdots,\lambda_{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

[TABLE]

Formula for inverse of Vandermonde matrix

The Vandermonde matrix is defined as

[TABLE]

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

[TABLE]

where $S_{p-i,j}=S_{p-i}(\lambda_{1},\cdots,\lambda_{j-1},\lambda_{j+1},\cdots,\lambda_{p}).$

Pulling out the common denominator, the $j$-th column vector of $V^{-1}$ is

[TABLE]

where the elementary symmetric polynomials are over variables $\lambda_{1},\cdots,\lambda_{j-1},\lambda_{j+1},\cdots,\lambda_{p}$.

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

[TABLE]

Bounding the condition number

As discussed before, the condition number for eigenvalue $\lambda_{j}$ is

[TABLE]

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}$.

By definition $V_{j,*}={\begin{bmatrix}1&\lambda_{j}&\lambda_{j}^{2}&\cdots&\lambda_{j}^{p-1}\end{bmatrix}}$, so

[TABLE]

Using the above explicit expression for $V^{-1}$, $\lVert V^{-1}_{*,j}\rVert_{2}=$

[TABLE]

.

Therefore,

[TABLE]

Note that both parts under $(\cdots)^{1/2}$ are greater than or equal to $1$, so we can bound it below by

[TABLE]

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

[TABLE]

While for the second part, since

[TABLE]

we have that

[TABLE]

Combining equation 17 and 18 for the upper bound, and putting it together with the lower bound,

[TABLE]

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.

Appendix 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.

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^{3}k^{2}T+mn^{4}k^{3})$, where $T$ is the sequence length, $n$ the hidden state dimension, $m$ the output dimension, and $k$ the input dimension.

Appendix 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 $\rho(A)\leq 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 $\ell_{2}$ distance $d(\Theta_{1},\Theta_{2})=\lVert\lambda(A_{1})-\lambda(A_{2})\rVert_{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 $\Theta=(A_{c},B_{c},C_{c})$, we generate $A^{\prime}$ by adding a i.i.d. Gaussians to each entry of $A_{c}$, while $B^{\prime}$ and $C^{\prime}$ are new random matrices of i.i.d. standard Gaussians. The standard deviation of the i.i.d. Gaussians for $A^{\prime}-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}\sim N(0,1)$ and put noise $\xi_{t}\sim 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.