EconBase
← Back to paper

Finite Time Identification in Unstable Linear Systems

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

48,972 characters · 8 sections · 39 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Finite Time Identification in Unstable Linear Systems

frontmatter, , \begin{keyword} Unstable Systems; Linear Dynamics; Finite Time Identification; Stabilization; Autoregressive Process; Non-Asymptotic Estimation \end{keyword} \begin{abstract} Identification of the parameters of stable linear dynamical systems is a well-studied problem in the literature, both in the low and high-dimensional settings. However, there are hardly any results for the unstable case, especially regarding {\em finite time bounds}. For this setting, classical results on least-squares estimation of the dynamics parameters are not applicable and therefore new concepts and technical approaches need to be developed to address the issue. Unstable linear systems arise in key real applications in control theory, econometrics, and finance. This study establishes finite time bounds for the identification error of the least-squares estimates for a fairly large class of heavy-tailed noise distributions, and transition matrices of such systems. The results relate the time length (samples) required for estimation to a function of the problem dimension and key characteristics of the true underlying transition matrix and the noise distribution. To establish them, appropriate concentration inequalities for random matrices and for sequences of martingale differences are leveraged. \end{abstract}

Introduction

Identification of the transition matrix in linear dynamical systems has been extensively studied in the literature for the stable case lutkepohl2005new,ljung1999system,soderstrom1989system. Further, new work has also addressed this topic under a high-dimensional scaling, with additional assumptions on sparsity of the parameters imposed on it basu2015regularized,zorzi2016ar,zorzi2017sparse. However, in settings where the underlying dynamics are {\em not stable}, this problem has {\em not} been adequately studied. A key issue that arises in this case is that the magnitude of the state vector explodes with high probability, exponentially over time lai1985asymptotic. Nevertheless, identification of the dynamics in the non-stable case is of interest due to a number of applications that give rise to such dynamics. In addition to adaptive control soderstrom2012discrete,kailath2000linear,kumar2015stochastic,bertsekas1995dynamic, these applications include a class of identification problems involving asset bubbles and high inflation episodes pesaran2005small,pesaran2010predictability,pesaran2002market,alogoskoufis1991phillips,garcia1991analysis,stock1996evidence,stock1998comparison,giacomini2006tests,nielsen2008explosive,engsted2006explosive,juselius2002high,nielsen2010analysis.

\textcolor{black}{Most existing work on the topic provides {\em asymptotic} results on the convergence lai1985asymptotic, as well as the limit distribution buchmann2007asymptotic,buchmann2013unified of the model parameters.} Specifically, early work investigated the limit distribution of the state vector under a set of restrictive assumptions on the dynamics matrix anderson1959asymptotic. Ensuing work dealt with the accuracy of identification in infinite time, for a class of structured transition matrices lai1983asymptotic. Further extensions to more general classes were established by Nielsen nielsen2005strong,nielsen2006order. Finally, additional asymptotic results together with the important concept of {\em irregularity} of the transition matrix which leads to inconsistency, are presented in the literature nielsen2009singular. However, finite time (i.e. non-asymptotic) results are not currently available.

In this work, we consider a linear dynamical system $x(t) \in \mathbb{R}^p, t=0,1,\cdots$ that evolves according to the following Vector Autoregressive (VAR) model

eqnarray[eqnarray omitted — 62 chars of source]

starting from an arbitrary initial state $x(0)$, \textcolor{black}{which can be either deterministic or stochastic}. Note that systems of longer but finite memory can also be written in the above form soderstrom2012discrete,kailath2000linear. We examine the general case where the system is not necessarily stable. The key contributions are: (i) establishing finite time identification bounds for the $\ell_2$ error of the least-squares estimates of the transition matrix $A_{0}$, (ii) under a fairly general heavy tailed noise (disturbance) process $\left\{w(t)\right\}_{t=1}^\infty$. In addition, the results due to the presence of a heavy-tailed noise term are of independent interest for the stable case as well. The novel results established provide insights on how the time length required for identification scales both with the dimension of the system, as well as with the characteristics of the transition matrix and the noise process.

In order to establish results for accurate finite time identification of $A_{0}$, one needs to address the following set of technical issues. Note that as long as $A_{0}$ has eigenvalues outside of the unit circle in the complex plane, the behavior of the Gram matrix of the state vector is governed by a random matrix. However, when $A_{0}$ has eigenvalues both inside and outside of the unit circle, the smallest eigenvalue of the Gram matrix scales linearly over time, while its largest eigenvalue grows exponentially, which in turn leads to the failure of the classical approaches to establish accurate identification. These issues are addressed in Subsections (ref) and (ref), respectively. In the proofs, we leverage selected concentration inequalities for random matrices tropp2012user, as well as an anti-concentration property of martingale difference sequences lai1983note.

The problem of fast accurate identification in unstable systems has a number of interesting applications. For example, in stochastic control, this includes the canonical problems of both stabilization, as well as design of an efficient adaptive policy for linear systems. First, since the dynamics are governed by unknown transition matrices, the control action can destabilize the system. Moreover, the user first needs to have an approximation of the dynamics, to be able to design a suitable control policy. Therefore, {\em accurate} identification of the dynamics of the transition matrices is necessary, even if they happen to lead to instability of the underlying system. More importantly, the identification result needs to be provided within a relative {\em short} time period for the user to be able to design the adaptive policy accordingly. More details are discussed in Example (ref).

\textcolor{black}{Applications of this setting in econometrics and finance also create the need to obtain finite time theoretic results. For example, in macroeconomics, the outstanding performance of the linear models marked them as a benchmark of forecasting the market pesaran2005small,stock1998comparison,giacomini2006tests. Their applications to the analysis of inflationary episodes in a number of OECD\footnote{Organization for Economic Co-operation and Development} countries pesaran2005small, as well as US stock prices engsted2006explosive,lin2017regularized are available in the literature. The former study establishes the structural non-stationarity of the process, where the latter verifies the explosive behavior of speculative bubbles. In particular, if a technology market is capable of important innovations with uncertain outcomes, it has been argued pesaran2010predictability that a bubble is very likely to emerge.}

\textcolor{black}{Another application involving unstable dynamics deals with episodes of hyperinflation. For example, Juselius and Mladenovic juselius2002high consider the case of (former) Yugoslavia and use data on various economic indicators to gain insights into the dynamics of the late 1990s episode. The analysis identifies wages, price level expectations, and currency depreciation as the key factors. In follow-up work, infinite time analysis techniques were used nielsen2010analysis, but as emphasized in the original work juselius2002high “hyperinflation episodes almost by definition are {\em short}." Therefore, the small sample size available can easily lead to incorrect inference, while finite time guarantees are informative about the sample size needed to make precise statements about the effects of different macroeconomic factors. Another hyperinflation episode from Germany in the early 1920's is studied by Nielsen nielsen2008explosive.}

Recently, the problem of forecasting non-stationary mixing kuznetsov2017generalization,kuznetsov2014generalization, and non-mixing kuznetsov2015learning time series has received attention, assuming the loss function employed is bounded. Unstable VAR models are a special, yet interesting, case of non-stationary time series. However, the problem of estimation/identification is not still addressed in the existing literature. Moreover, the results on forecasting are not applicable to the identification problem, since the least-squares loss function used in that study is not bounded. On the other hand, the obtained results on identification are applicable to forecasting.

The remainder of the paper is organized as follows. In Section (ref) we provide a rigorous formulation of the problem, introduce the identification procedure, and outline examples that require accurate identification but the system can not assumed to be stable. The contributions are discussed in Section (ref), where we study different scenarios. First, we provide identification results on (non-stationary) stable linear systems in Subsection (ref), followed by the explosive case (Subsection (ref)). Finally, we study the accurate identification of the dynamics for general systems in Subsection (ref).

Notations

The following {notation} is used throughout this paper. For a matrix $A \in \mathbb{C}^{p \times q}$, $A'$ denotes its transpose. When $p=q$, the smallest (respectively largest) eigenvalue of $A$ (in magnitude) is denoted by $\lambda_{\min} (A)$ (respectively $\lambda_{\max}(A)$). For $\gamma \in \mathbb{R}, \gamma > 0, x \in \mathbb{C}^q$, define the norm ${\left\vert\kern-0.25ex\left\vert x \right\vert\kern-0.25ex\right\vert}_{\gamma} = \left( \sum\limits_{i=1}^{q} \left| x_i \right|^\gamma \right)^{1/\gamma}$. For $\gamma = \infty$, define the norm ${\left\vert\kern-0.25ex\left\vert x \right\vert\kern-0.25ex\right\vert}_{\infty} = \max \limits_{1 \leq i \leq q} |x_i|$.

\textcolor{black}{We also use the following notation for the operator norm of matrices. For $\beta, \gamma \in \left(0,\infty\right], A \in \mathbb{C}^{p \times q}$ let,

equation*[equation* omitted — 383 chars of source]

Whenever $\gamma = \beta$, we simply write ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert A \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{\beta}$. To show the dimension of manifold $M$ over the field $F$, we use $\mathrm{dim}_{F} \left(M\right)$. The sigma-field generated by random vectors $X_1,\cdots,X_n$ is denoted by $\sigma \left( X_1,\cdots,X_n \right)$. Finally, the symbol $\vee$ denotes the maximum of two or more quantities.}

Problem Formulation and Preliminaries

The system $\left\{x(t)\right\}_{t=0}^\infty$ evolves according to (ref), while the unknown transition matrix $A_{0} \in \mathbb{R}^{p \times p}$ is not assumed to be stable, i.e. the eigenvalues of $A_{0}$ do not necessarily lie inside the unit circle. Further, $\left\{w(t)\right\}_{t=1}^\infty$ is the sequence of independent mean-zero noise vectors with covariance matrix $C$, i.e. $\mathbb{E} \left[w(t)\right]=0$, and $\mathbb{E} \left[w(t)w(t)'\right]=C$.

remarkThe results established also hold if the noise vectors are martingale difference sequences. Further, the generalization to heteroscedastic noise, where the covariance matrix $C$ is time varying, is rather straightforward.

\textcolor{black}{The objective is to identify $A_{0}$, using the least-squares estimator. One observes the state vector during a finite time interval, $\left\{x(t)\right\}_{t=0}^n$, and defines the sum-of-squares loss function

eqnarray*[eqnarray* omitted — 172 chars of source]

Then, $A_{0}$ is estimated by $\hat{A}^{(n)}$, which is the minimizer of the above sum-of-squares; $\mathcal{L}_{n} \left(\hat{A}^{(n)}\right) = \min\limits_{A_{} \in \mathbb{R}^{p \times p }} \mathcal{L}_{n} \left(A_{}\right)$.}

The {main} contribution of this paper is to establish that with high probability, accurate identification of the true transition matrix is achieved, excluding a pathological case. Formally, \textcolor{black}{for arbitrary accuracy $\epsilon>0$ and failure probability $\delta>0$,} $\hat{A}^{(n)}$ is with probability at least $1-\delta$ within an $\epsilon$-neighborhood of $A_{0}$, where apart from a logarithmic factor, the time length $n$ scales quadratically with $\epsilon^{-1}$, and logarithmically with $\delta^{-1}$. \textcolor{black}{In other words, for a fixed accuracy $\epsilon>0$, the probability that the identification error ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \hat{A}^{(n)}-A_{0} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{2}$ exceeds $\epsilon$, decays {\em exponentially} as $n$ grows.}

The following example elaborates on the problem of finite time identification for unstable dynamical systems in control theory.

example[Stabilization in adaptive control] Consider the linear stochastic system $\left[ A_{x}, A_{u}\right]$, where the state evolution is governed by the following dynamics: \begin{eqnarray*} x(t+1)=A_{x}x(t)+A_{u}u(t)+w(t+1). \end{eqnarray*}

In the previous equation, the vector $x(t) \in \mathbb{R}^p$ represents the state of the system, and $u(t) \in \mathbb{R}^r$ is the control action taken by the user. The {\em unknown} transition matrix $A_{x} \in \mathbb{R}^{p \times p}$ determines the evolution of the system, and the {\em unknown} input matrix $A_{u} \in \mathbb{R}^{p \times r}$ shows the effect of the control policy on the state of the system.

Due to the simplicity of the structure, the main interest is in linear feedbacks of the form $u(t)=Lx(t)$, where $L \in \mathbb{R}^{r \times p}$ is the feedback matrix. Further, in addition to preserving the linear nature of the system (which prevents the analysis from becoming mathematically intractable), linear feedbacks correspond to important objectives for a class of optimal control problems soderstrom2012discrete,kumar2015stochastic, including minimization of quadratic costs bertsekas1995dynamic,brunner2018stochastic. So, the linear dynamics are essentially determined by the closed-loop transition matrix $A_{0}=A_{x}+A_{u}L$.

\textcolor{black}{The system $\Theta_0=\left[ A_{x}, A_{u} \right]$ is assumed to be stabilizable, implying there exists a stabilizer $L_0$ such that the closed-loop matrix $A_{x}+A_{u}L_0$ is stable; $\left| \lambda_{\max} \left( A_{x}+A_{u}L_0 \right)\right|<1$. Finding such a stabilizer requires precise approximation of the true dynamics $\Theta_0$ bertsekas1995dynamic, as shown in the following example. Consider a system of dimension $p=3, r=2$, which is stabilizable, since exact knowledge of $\Theta_0$ yields $\left| \lambda_{\max} \left( A_{0} \right)\right|= 0.22$. Fig. (ref) depicts the scatter plot of the largest eigenvalue of the closed-loop transition matrix versus the relative magnitude of an Additive White Gaussian Noise (AWGN). A stabilizing linear feedback $L$ is applied to the system as if the dynamics parameter is $\Theta_0 + \Delta $ instead of $\Theta_0$, where entries of $\Delta$ are independent Gaussian measurement errors. It can be seen that a measurement error as small as $5\%$ in the identification of the system dynamics can lead to instability.}

\textcolor{black}{Fig. (ref) graphs the largest eigenvalue of the closed-loop matrix versus a perturbation in a single entry of $\Theta_0$. In fact, for different entries of $\Theta_0$, the linear feedback $L$ is designed as if the operator approximates a single entry incorrectly. Formally, for $\epsilon \geq 0$, only the $(i,j)$-th entry of $\Theta_0$ is approximated with error $\epsilon$, while all other entries are exactly provided to the operator. Fig. (ref) corresponds to the relationship between $\left| \lambda_{\max} \left( A_{0} \right)\right|$ and $\epsilon$, for different entries $(i,j)$. Therefore, stabilization is very sensitive to the perturbation, as an error of $3\%$ in relative magnitude in a single element of the system will totally destabilize the system. In many applications, especially if the system under consideration is not man-made, such precise information is not available. Hence, the matrix $A_{0}$ can not be assumed to be a priori stable.}

figure[figure omitted — 353 chars of source]
figure[figure omitted — 286 chars of source]

In addition, in order to design a desired policy (either steering the system to a specific state faradonbeh2016optimality or minimizing a cost function bertsekas1995dynamic), such an approximation is necessary. To obtain it, learning accurately the dynamics of an unstable system is needed. Importantly, such learning needs to conclude in finite time, because afterwards, the user needs to control the system to achieve the corresponding objective, determined by the application.

In order to establish high probability guarantees for accurate identification of the closed-loop matrix, we apply the results from Theorem (ref) given in the next Section. A random linear feedback, denoted by $L$, suffices to satisfy the assumptions of Theorem (ref) for the closed-loop matrix $A_{x}+ A_{u} L$. In fact, it suffices for $L$ to be a continuously distributed random matrix. This in turn, leads to accurate identification of $\left[ A_{x}, A_{u} \right]$, applying multiple random linear feedbacks, drawn independently. Note that direct identification of $\left[ A_{x}, A_{u} \right]$ is infeasible, since by observing the state sequence $\left\{ x(t) \right\}_{t=0}^\infty$, the best result one can provide is “closed-loop identification" soderstrom2012discrete,kumar1990convergence. \textcolor{black}{Specifically, for a given closed-loop transition matrix $A_{0}$, the set of parameters guiding the system's dynamics $A_{x},A_{u}$ which satisfy $A_{0}=A_{x}+A_{u}L$ is not unique if one exactly knows the feedback matrix $L$. This set is indeed a subspace of dimension $pr$ in the space $\mathbb{R}^{p \times (p+r)}$ the matrix $\left[A_{x},A_{u}\right]$ belongs to.}

To analyze the finite time behavior of the aforementioned identification procedure, the following is assumed for the tail-behavior of every coordinate of the noise vector.

assump[Sub-Weibull noise distribution] There exist positive constants $b, d$, and $\alpha$, such that for all $t=1,2,\cdots; i=1, \cdots, p; y >0$, \begin{eqnarray*} \mathbb{P} \left(\left|w_i(t)\right| > y\right) \leq b \: \mathrm{exp} \left(-\frac{y^\alpha}{d}\right). \end{eqnarray*}

\textcolor{black}{In case of random initial state $x(0)$, we assume it also follows a sub-Weibull distribution.} Intuitively, smaller values of the exponent $\alpha$ correspond to heavier tails for the noise distribution, and vice versa. Note that assuming a sub-Weibull distribution for the noise coordinates is more general than the sub-Gaussian (or sub-exponential) assumption routinely made in the literature abbasi2011online, where $\alpha \geq 2$ ($\alpha \geq 1$). In fact, when $\alpha <1$, the noise coordinates $w_i(t)$ do not need to have a moment generating function.

\textcolor{black}{Note that for establishing consistency of infinite time identification procedures, the noise vectors need to satisfy a moment condition, e.g. $\mathbb{E} \left[{\left\vert\kern-0.25ex\left\vert w(t) \right\vert\kern-0.25ex\right\vert}_{2}^{2+\alpha}\right]<\infty$, for some $\alpha>0$ lai1985asymptotic,lai1983asymptotic. On the other hand, finite time identification analysis results are usually obtained under an assumption of a light-tail (or even uniformly bounded) noise distribution; e.g. Gaussian process tropp2012user,abbasi2011online. Thus, the above assumption on sub-Weibull noise, that includes a family of heavy-tailed noise processes, provides a fairly general framework to narrow down the theoretical gap between asymptotic and non-asymptotic approaches. }

The noise coordinates can be either discrete or continuous random variables, and are not assumed to have a probability density function. To proceed, we define a property of the population covariance matrix of the system under study. It is easy to see that the following property is necessary and sufficient for accurate estimation of dynamics parameters.

deffn[Reachability] The pair $\left[A_{0},C\right]$ is called reachable if \begin{eqnarray*} \mathrm{rank}\left(\left[C^{1/2}, A_{0}C^{1/2}, \cdots, A_{0}^{p-1}C^{1/2} \right]\right)=p. \end{eqnarray*}

Clearly, reachability is equivalent to $\left| \lambda_{\min} \left( K(C) \right)\right|>0 $, where $K(C)=\sum\limits_{i=0}^{p-1}A_{0}^i C {A_{0}'}^i$. Specifically, if $C$ is positive definite, then $\left[A_{0},C\right]$ is reachable for all $A_{0} \in \mathbb{R}^{p\times p}$. Reachability is conceptually equivalent to the population covariance matrix of the system being positive definite. More precisely, since the noise vectors are independent, the covariance matrix of $x(t)$ is given by $\sum\limits_{i=0}^{t-1}A_{0}^i C {A_{0}'}^i$; i.e. reachability is in fact stating that for $t \geq p$, every coordinate of $x(t)$ has non-degenerate randomness.

\textcolor{black}{Further, reachability is particularly helpful if the actual evolution of the system is guided by VAR$\left( k \right)$ dynamics, for some $k>1$. In this case, the next step is determined by the $k$ previous lags: for $t \geq k$, the state sequence $\tilde{x}(t) \in \mathbb{R}^m$ evolves according to

eqnarray*[eqnarray* omitted — 93 chars of source]

for some initial vectors $\tilde{x}(0), \cdots, \tilde{x}(k-1) \in \mathbb{R}^m$, and transition matrices $A_{1}, \cdots, A_{k} \in \mathbb{R}^{m \times m}$, assuming $A_{k} \neq 0$. Arranging blocks of $\tilde{x}(t)$ accordingly, $x(t) = \left[ \tilde{x}(t+k-1)' , \cdots , \tilde{x}(t)' \right] ' \in \mathbb{R}^{km}$, the state evolution can be written in the form of (ref), for $A_{0}=

bmatrix[bmatrix omitted — 61 chars of source]

\in \mathbb{R}^{km \times km}$. Then, as long as the covariance matrix of $\tilde{w}(t)$ is full rank, reachability holds.}

Main results

Next, we establish the key identification results that characterize the time (samples) required, so that with high probability the $A_{0}$ least-squares estimate is accurate within a certain degree. First, we study the identification for stable systems where all eigenvalues of $A_{0}$ are inside the unit circle, i.e. $\left| \lambda_{\max} \left( A_{0} \right)\right|<1$. Subsequently, the explosive case where all eigenvalues of the transition matrix $A_{0}$ lie outside of the unit circle, i.e. $\left| \lambda_{\min} \left( A_{0} \right)\right|>1$, is examined. Finally, finite time identification results are presented for the general case which is the combination of these two regimes.

Some straightforward algebra shows that the least-squares estimator can be written as

eqnarray*[eqnarray* omitted — 77 chars of source]

where $V_n=\sum\limits_{t=0}^{n-1} x(t)x(t)'$ denotes the empirical covariance matrix of the state process \textcolor{black}{(once normalized by $n$)}, which is assumed to be non-singular.

The latter result implies that the behavior of $V_n$ needs to be carefully studied and this constitutes a major part of the following two subsections.

Stable systems

The stable case has been extensively studied before, customarily under the stronger assumption of sub-Gaussian noise abbasi2011online. Next, we generalize the results to sub-Weibull noise vectors defined in Assumption (ref). Further, these results will be used for the general case in Subsection (ref).

\textcolor{black}{For a stable transition matrix $A_{0} \in \mathbb{R}^{p \times p}$, we define the constant $\eta_{}\left(A_{0}\right)$, that is critical in specifying various constants that appear in the main results. Its definition is based on the Jordan decomposition of square matrices.}

\textcolor{black}{First, for $\lambda \in \mathbb{C}$, define the size $m$ Jordan matrix of $\lambda$ as follows.

eqnarray*[eqnarray* omitted — 310 chars of source]

Then, the Jordan decomposition of $A_{0}$ is given by $A_{0}=P^{-1}\Lambda P$, where $\Lambda$ is block diagonal, $\Lambda= \mathrm{diag}\left(\Lambda_1,\cdots, \Lambda_k\right)$, with $\Lambda_i \in \mathbb{C}^{m_i \times m_i}, i=1, \cdots, k$ being a Jordan matrix of $\lambda_i$.

deffnFor a stable matrix $A_{0}$, suppose that $A_{0}=P^{-1} \Lambda P$ is the Jordan decomposition as described above. For $t=1,2, \cdots$, let \begin{eqnarray*} \eta_{t}\left(\Lambda_i\right)=\inf\limits_{\rho \geq \left|\lambda_i\right|} t^{m_i-1} \rho^t \sum\limits_{j=0}^{m_i-1} \frac{\rho^{-j}}{j!}, \end{eqnarray*} and $\eta_{t}\left(\Lambda\right)= \max \limits_{1 \leq i \leq k} \eta_{t}\left(\Lambda_i\right)$. Then, letting $\eta_{0}\left(\Lambda\right)=1$, define \begin{eqnarray*} \eta_\left(A_{0}\right) = {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert P^{-1} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{\infty \to 2} {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert P \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{\infty} \sum\limits_{t=0}^{\infty} \eta_{t}\left(\Lambda\right). \end{eqnarray*}

Clearly, denoting the largest algebraic multiplicity of the eigenvalues of $A_{0}$ (which is the same to the largest block-size in the Jordan form) by $\mu\left(A_{0}\right)=\max\limits_{1 \leq i \leq k} m_i$, we have

eqnarray[eqnarray omitted — 207 chars of source]

}

In the stable regime, the state process has a stationary limit distribution. In this case, the empirical covariance matrix has an approximately deterministic behavior, which is described by its asymptotic distribution. Specifically, as time grows, $V_n$ appropriately normalized, can be approximated by $\kappa \left(C\right)$, where $\kappa \left(C\right)= \sum\limits_{i=0}^{\infty} A_{0}^i C {A_{0}'}^i$ denotes the asymptotic covariance matrix.

\textcolor{black}{The following lemma provides a finite lower bound for the time length (number of samples), based on the identification error $\epsilon$, and the failure probability $\delta$. For this purpose, using the parameters $b, d$, and $\alpha$ specified in Assumption (ref), we define the following constant. Henceforth, one can let $\alpha \to \infty$, if the noise vectors $w(1),w(2), \cdots$ are uniformly bounded.

eqnarray*[eqnarray* omitted — 470 chars of source]

}

lemmAssuming $\left| \lambda_{\max} \left( A_{0} \right)\right|<1$, let $c_1 $ be as defined above. Then, for arbitrary $\epsilon, \delta>0$ if \begin{eqnarray*} \frac{n}{\left(\log n\right)^{4/\alpha}} \geq \frac{c_1}{\epsilon^2} \left(-\log \delta\right)^{1+4/\alpha}, \end{eqnarray*} then \begin{eqnarray*} \mathbb{P} \left(\left| \lambda_{\max} \left( \frac{1}{n}V_{n+1}- \kappa \left(C\right) \right)\right| > \epsilon\right) \leq \delta. \end{eqnarray*}

A direct consequence of Lemma (ref) is the following corollary, which shows that high probability accurate identification can be ensured, if reachability, as defined in Definition (ref), is assumed. Note that reachability implies that $\kappa \left(C\right)$ is positive definite. \textcolor{black}{Using $c_1$ defined above, we define $c_2= 4 c_1 \left( {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert A_{0} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{2}^2 \vee 1 \right) \left| \lambda_{\min} \left( K(C) \right)\right|^{-2}+2$.}

corolSuppose that $\left| \lambda_{\max} \left( A_{0} \right)\right|<1$, and $\left[A_{0},C\right]$ is reachable. Then, for $c_2$ above, and for all $\epsilon,\delta>0$, \begin{eqnarray*} \frac{n}{\left(\log n\right)^{4/\alpha}} \geq \frac{c_2}{\epsilon^2} \left(-\log \delta\right)^{1+4/\alpha}, \end{eqnarray*} implies \begin{eqnarray*} \mathbb{P} \left({\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \hat{A}^{(n)}-A_{0} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{2}>\epsilon\right) < \delta. \end{eqnarray*}

Explosive systems

In the explosive case, the empirical covariance matrix $V_n$ grows exponentially with respect to $n$. In addition, unlike the stable case, $V_n$ appropriately normalized, can be approximated by a random matrix. Therefore, the eigenvalues of the normalized empirical covariance matrix are stochastic as well. In order to find deterministic bounds for the eigenvalues of $V_n$, new quantities, denoted by $\phi \left(A_{0}\right), \psi \left(A_{0},\delta\right)$, need to be defined.

Subsequently, after providing formal definitions of these quantities, we present in Lemma (ref) bounds for the eigenvalues. Then, a sufficient and necessary property of $A_{0}$ for accurate identification is introduced, followed by Propositions (ref), (ref), which establish the positiveness of $\phi \left(\cdot\right), \psi \left(\cdot,\cdot\right)$. This subsection concludes with Corollary (ref) that deals with identification in explosive systems.

First, for explosive $A_{0}$, we define the nonnegative functions $\phi \left(A_{0}\right), \psi \left(A_{0},\delta\right)$ as follows. Assuming $\left| \lambda_{\min} \left( A_{0} \right)\right|>1$, let $A_{0}=P^{-1}\Lambda P$ be the Jordan decomposition. Letting

eqnarray*[eqnarray* omitted — 120 chars of source]

for $\delta>0$ define

eqnarray*[eqnarray* omitted — 190 chars of source]

Note that according to this definition, all coordinates of the vector $Pz(\infty)$ are in magnitude at least $\psi \left(A_{0},\delta\right)$, with probability at least $1-\delta$. Next, define

eqnarray*[eqnarray* omitted — 386 chars of source]

where for an arbitrary matrix $M \in \mathbb{C}^{m \times k}$, $\left[M\right]_{\min}$ is the smallest magnitude of the nonzero entries of $M$:

eqnarray*[eqnarray* omitted — 129 chars of source]

In fact, $\phi \left(A_{0}\right)$ represents the deterministic portion of the smallest eigenvalue of the random matrix $F_\infty$ which approximates the normalized $V_n$. It only depends on $A_{0}$, while $\psi \left(A_{0},\delta\right)$ represents the stochastic portion which depends on both $A_{0}$ and the distribution of the noise sequence $\left\{w(t)\right\}_{t=1}^\infty$. Intuitively, $\phi \left(A_{0}\right)$ denotes the minimum nontrivial distance between the polynomials of $A_{0}^{-1}$ and the origin, and $\psi \left(A_{0},\delta\right)$ denotes the high probability minimum distance of the vector $Pz(\infty)$ from the origin. These minimum distances show up, because for $v \in \mathbb{R}^p$, $v'F_\infty v$ is determined by the product of a polynomial of $A_{0}^{-1}$ (with coefficients determined by $v$), and $Pz(\infty)$. More details are provided in the proof of Lemma (ref).

Now, the behavior of the normalized empirical covariance matrix can be controlled as follows:

lemmSuppose that $\left| \lambda_{\min} \left( A_{0} \right)\right|>1$; then, there is a constant $\xi \left(A_{0}\right) < \infty$ such that for all $n,\delta$, \begin{eqnarray*} \mathbb{P} \left(\left| \lambda_{\max} \left( A_{0}^{-n} V_{n+1} {A_{0}'}^{-n} \right)\right| > \xi \left(A_{0}\right) \left(-\log\delta\right)^{2/\alpha}\right) \leq \delta. \end{eqnarray*} \textcolor{black}{Further, there is a constant $n_1 < \infty$, such that for arbitrary $\epsilon,\delta>0$ if \begin{eqnarray} n \geq \frac{3 \left( \alpha+2 \right)}{\alpha \log \left| \lambda_{\min} \left( A_{0} \right)\right|} \log \left(\frac{- \log \delta}{\epsilon} \right) \vee n_1, \end{eqnarray} then with probability at least $1-4\delta$ it holds that \begin{eqnarray} \left| \lambda_{\min} \left( A_{0}^{-n} V_{n+1} {A_{0}'}^{-n} \right)\right| \geq \phi \left(A_{0}\right)^2\psi \left(A_{0},\delta\right)^2 - \epsilon. \end{eqnarray}}
remarkThe inequality (ref) is of interest for the following two reasons. First, the accuracy $\epsilon$ decays exponentially fast when $n$ grows. Second, the failure probability $\delta$ decays {\em doubly} exponentially fast with $n$.

This surprising strong behavior is intuitively caused by the exponential growth of $x(t)$. Broadly speaking, the growing signal (i.e. $x(t)$) to noise (i.e. $w(t)$) ratio leads to the super fast decay of $\epsilon$ and $\delta$. Note that commonly in identification problems, the decay rates of $\epsilon, \delta$ are square-root, and exponential, respectively.

If $\phi \left(A_{0}\right) \psi \left(A_{0},\delta\right)=0$, obviously (ref) holds. Thus, the main interest is in the case where $\phi \left(A_{0}\right) \psi \left(A_{0},\delta\right) \neq 0$, which we will show that holds under certain conditions, and is necessary to ensure accurate identification. In fact, the first case is of no interest, since it can be shown that $V_n$ will be singular, and thus identification of $A_{0}$ fails, even if the time period becomes infinitely large nielsen2009singular. For the second case, the transition matrix $A_{0}$ needs to be regular, according to the following definition. Regularity (of course in addition to reachability), leads to accurate identification as shown in Corollary (ref).

deffn[Regularity] $A_{} \in \mathbb{R}^{p \times p}$ is called regular if for any explosive eigenvalue of $A_{}$, denoted by $\lambda$, the geometric multiplicity of $\lambda$ is one.

Regularity essentially implies that the eigenspace corresponding to $\lambda$ is one dimensional. There are also equivalent formulations for regularity. Indeed, $A_{}$ is regular, if and only if for any explosive eigenvalue $\lambda$, in the Jordan decomposition of $A_{}$ there is only one block corresponding to $\lambda$. In other words, no matter how large the algebraic multiplicity of $\lambda$ is, its geometric multiplicity is one. Another equivalent formulation is the following one. $A_{}$ is regular if and only if

eqnarray*[eqnarray* omitted — 70 chars of source]

for all $\lambda \in \mathbb{C}$ such that $\left|\lambda \right|>1$. For example, let $P_1,P_2 \in \mathbb{C}^{2 \times 2}$ be arbitrary invertible matrices, and

eqnarray*[eqnarray* omitted — 171 chars of source]

where $\rho \in \mathbb{C}, \left|\rho\right|>1$. Then, $A_{1}$ is regular, where $A_{2}$ is not.

propoAssuming $\left| \lambda_{\min} \left( A_{0} \right)\right|>1$, regularity of $A_{0}$ is equivalent to $\phi \left(A_{0}\right)>0$.

The next proposition shows that positiveness of $\psi \left(A_{0},\delta\right)$ is implied by reachability. Proposition (ref) also reveals a linear scaling of $\psi \left(A_{0},\delta\right)$ with respect to $\delta$, when the noise is a continuous random variable.

propoAssume $\left| \lambda_{\min} \left( A_{0} \right)\right|>1$, and $\left[A_{0},C\right]$ is reachable. We then have $\psi \left(A_{0},\delta\right)> 0$. Moreover, if there is $i \geq p$, such that $w(i-p+1), \cdots, w(i)$ have bounded probability density functions (pdf) over certain subspaces of $\mathbb{R}^p$, then, \begin{eqnarray*} \psi \left(A_{0},\delta\right) \geq \psi \left(A_{0}\right) \delta, \end{eqnarray*} for some constant $\psi \left(A_{0}\right) >0$. If the bounded pdfs mentioned above correspond to the normal distribution, then \begin{eqnarray*} \psi \left(A_{0}\right) \geq \left(\frac{\pi \left| \lambda_{\min} \left( K(C) \right)\right|}{2\left| \lambda_{\max} \left( {A_{0}}^{i}{A_{0}'}^{i} \right)\right|}\right)^{1/2} p^{-1} \left(\min\limits_{1 \leq i \leq p}{\left\vert\kern-0.25ex\left\vert P_i \right\vert\kern-0.25ex\right\vert}_{2}\right) . \end{eqnarray*}

Now, we are ready to state the key result for the time length required to achieve accurate estimation for an explosive transition matrix.

corolSuppose that $\left| \lambda_{\min} \left( A_{0} \right)\right|>1$, $A_{0}$ is regular, and $\left[A_{0},C\right]$ is reachable. There exists a constant $n_2 <\infty$, such that for all $\epsilon,\delta>0$, \textcolor{black}{\begin{eqnarray} n \geq \frac{3 \left( \alpha + 4 \right)}{\alpha \log \left| \lambda_{\min} \left( A_{0} \right)\right|} \log \left(\frac{- \log \delta}{\epsilon \psi \left(A_{0},\delta\right)} \right) \vee n_2 \end{eqnarray}} implies \begin{eqnarray*} \mathbb{P} \left({\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \hat{A}^{(n)}-A_{0} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{2}>\epsilon\right) < 4\delta. \end{eqnarray*}

The time length specified in (ref) is similar to that of Lemma (ref) in terms of the accuracy $\epsilon$, while the dependence in $\delta$ is different. In fact, compared to (ref), the decay rate of $\delta$ is of the common exponential order as $n$ grows (assuming the linear scaling of $\psi \left(A_{0},\delta\right)$ with respect to $\delta$).

remarkAnother interesting property of explosive systems is that $n_1,n_2$ scale logarithmically with respect to the dimension $p$.

The dependency of $n_1,n_2$ on $A_{0}$, as well as $b, d$, and $\alpha$ specified in Assumption (ref) are outlined explicitly in the corresponding proofs. Moreover, the constants $n_1,n_2$ are in fact universal; for $\rho_1>0$, a single $n_1$ depending on $\rho_1$ implies (ref) for all matrices $A_{0}$ satisfying $\left| \lambda_{\min} \left( A_{0} \right)\right| \geq 1+\rho_1$.

In addition, let $\lambda_1 \left(A_{}\right), \cdots, \lambda_{k \left(A_{}\right)} \left(A_{}\right)$ be the distinct eigenvalues of $A_{}$. Then, there is a single universal constant $n_2$ depending on $\rho_1,\rho_2$, such that (ref) implies the desired estimation result of Corollary (ref), for all regular explosive transition matrices $A_{0}$ satisfying $1 + \rho_1 \leq \left| \lambda_{\min} \left( A_{0} \right)\right|$, and $0< \rho_2 \leq \min\limits_{1 \leq i < j \leq k \left( A_{0} \right)} \left| \lambda_i\left( A_{0} \right) - \lambda_j\left( A_{0} \right) \right|$.

General systems

The previous results enable us to establish the key result of the paper. Theorem (ref) establishes the accuracy of identification, when the regular matrix $A_{0}$ has no eigenvalue on the unit circle. \textcolor{black}{As the following well known fact states, this assumption includes almost all matrices mumford1999red.}

factThe set of all $p \times p$ real matrices with at least one eigenvalue on the unit circle has Lebesgue measure zero. Moreover, almost all matrices are regular.

\textcolor{black}{However, note that transition matrices with unit eigenvalues occur in applications, including resonating mechanical systems fossen2011parametric, the study of macroeconomic indicators engsted2006explosive,nelson1982trends and the timeline of bubbles during the crisis in the mid-late 2000s phillips2011dating. Therefore, addressing the identification problem for unit root transition matrices, even though they constitute a measure zero set, is an interesting direction for future work.}

Excluding two pathological cases of square matrices with at least one eigenvalue on the unit circle, and irregular matrices, the estimation of the transition matrix for a general unstable system is with high probability arbitrarily accurate, as determined in the following theorem. A well known fact states that there is an invertible matrix $M \in \mathbb{R}^{p \times p}$, such that $\tilde{A}=MA_{0}M^{-1} \in \mathbb{R}^{p \times p}$ is a block diagonal matrix,

eqnarray*[eqnarray* omitted — 88 chars of source]

where for $i=1,2$, we have $A_{i} \in \mathbb{R}^{p_i \times p_i}$, $p_1+p_2=p$, and

eqnarray*[eqnarray* omitted — 122 chars of source]

Technically, $p_1$ ($p_2$) is sum of the algebraic multiplicities of the stable (explosive) eigenvalues of the true unknown matrix $A_{0}$. Conceptually, it determines the dimension of a certain subspace of $\mathbb{R}^p$, on which the linear transformation $A_{0}$ is stable (explosive). Note that since $M$ is not known in advance, the above split of the true transition matrix to a stable one and an explosive one cannot be used in the identification procedure.

thrmSuppose that $A_{0}$ is regular, has no unit eigenvalue, $\left[A_{0},C\right]$ is reachable, and $A_{2}$ is as above. Then, there exist constants $c_3, n_3 <\infty$, such that for all $\epsilon, \delta>0$, \textcolor{black}{\begin{eqnarray} \frac{n}{\left(\log n\right)^{4/\alpha}} \geq \frac{c_3}{\epsilon^2} \left( \left( -\log \delta \right)^{1+4/\alpha} - \log \psi \left(A_{2},\delta\right) \right ) \vee n_3, \:\:\:\:\:\: \end{eqnarray}} implies that \begin{eqnarray*} \mathbb{P} \left({\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \hat{A}^{(n)}-A_{0} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_{2}>\epsilon\right) < 6\delta. \end{eqnarray*}

\textcolor{black}{Regarding the time length above in (ref), the exact specification of the constants $c_3, n_3$ requires some additional definitions, provided in the proof of Theorem (ref). Broadly speaking, the behavior of $c_3$ (respectively $n_3$) is similar to that of $c_2$ (resp. $n_2$) used in Corollary (ref) (resp. (ref)). Note that in order to compute $c_3$ (resp. $n_3$), one has to use the stable (resp. explosive) matrix $A_{1}$ (resp. $A_{2}$). Further, since $A_{0}$ is regular, regularity of $A_{2}$ is automatically guaranteed. Therefore, Corollaries (ref) and (ref) can be used. Note that the reachability condition is inherited from the matrix $A_{0}$, as formally presented in Proposition (ref). Thus, Proposition (ref) implies that $-\log \psi \left(A_{2},\delta\right) < \infty$, and it is up to a constant less than $-\log \delta$, if the noise vectors have bounded probability density functions. Therefore, using Proposition (ref), for continuously distributed noise vectors with bounded pdfs one can substitute (ref) with

eqnarray*[eqnarray* omitted — 187 chars of source]

}

Concluding Remarks

We studied the problem of providing finite time bounds for the least-squares estimates of general linear dynamical systems, where the transition matrix does not necessarily need to be stable. The relationships between different parameters involved, including time length, accuracy of the identification, failure probability, the transition and noise matrices, and dimension are investigated. We prove that apart from a pathological case of zero Lebesgue measure, the identification is with high probability accurate, if the length of the time period scales similar to standard results in estimation theory, i.e. quadratic scaling with the inverse identification error and logarithmic scaling with the failure probability.

These finite time results for such a widely used model can be helpful to obtain analogous results for more complicated models exhibiting temporal dependence, such as nonlinear systems. Further, the techniques used in this work to analyze the accuracy when the systems under study are not necessarily stable, provide insight for settings where additional knowledge on the structure of the dynamics is available. In particular, potential extensions to a high-dimensional setting (assuming that the transition matrix is sparse), or other structured classes such as low-rank matrices, \textcolor{black}{as well as addressing practically interesting cases of null measure, are topics of interest and for future investigation.}