Finite Sample Analysis of Stochastic System Identification

TL;DR

This paper provides a finite sample analysis for stochastic system identification, establishing non-asymptotic error bounds that decrease with more data and are valid even for marginally stable systems.

Contribution

It introduces a non-asymptotic analysis framework for system identification using modern probabilistic tools, extending classical results to finite samples and marginally stable systems.

Findings

Estimation errors decrease at a rate of 1/√N with high probability.

Non-asymptotic bounds align with classical asymptotic results.

Results are valid even for marginally stable systems.

Abstract

In this paper, we analyze the finite sample complexity of stochastic system identification using modern tools from machine learning and statistics. An unknown discrete-time linear system evolves over time under Gaussian noise without external inputs. The objective is to recover the system parameters as well as the Kalman filter gain, given a single trajectory of output measurements over a finite horizon of length $N$. Based on a subspace identification algorithm and a finite number of $N$ output samples, we provide non-asymptotic high-probability upper bounds for the system parameter estimation errors. Our analysis uses recent results from random matrix theory, self-normalized martingales and SVD robustness, in order to show that with high probability the estimation errors decrease with a rate of $1/\sqrt{N}$. Our non-asymptotic bounds not only agree with classical asymptotic results, but are also valid even when the system is marginally stable.

Equations

$$x_{k+1}$$

$$y_k$$

$${\hat{x}}_{k+1}$$

$$y_k$$

$$K=AP{C}^{\ast }(CP{C}^{\ast }+R{)}^{-1}$$

$$P=AP{A}^{\ast }+Q-AP{C}^{\ast }(CP{C}^{\ast }+R{)}^{-1}CP{A}^{\ast }$$

$${\hat{x}}_{k+1}=\mathbb{E}\left[x_{k+1}|y_0,\dots ,y_k\right],{\hat{x}}_0=0.$$

$${\Gamma }_k=\mathbb{E}\left[{\hat{x}}_k{\hat{x}}_k^{\ast }\right]$$

$$\bar{R}\triangleq \mathbb{E}\left[e_k{e}_k^{\ast }\right]=CP{C}^{\ast }+R$$

$${\Vert \hat{A}-S^{-1}AS\Vert }_2\le {ϵ}_A\left(\delta ,N\right)$$

$${\Vert \hat{C}-CS\Vert }_2\le {ϵ}_C\left(\delta ,N\right)$$

$${\Vert \hat{K}-S^{-1}K\Vert }_2\le {ϵ}_K\left(\delta ,N\right),$$

$$Y_k^{+}$$

$$Y_{+}$$

$$Y_{-}$$

$$\hat{X}$$

$$\mathcal{O}_{k}\triangleq\left[\begin{array}[]{cccc}C^{*}&A^{*}C^{*}&\cdots&(A^{*})^{k-1}C^{*}\end{array}\right]^{*},$$

$$\mathcal{K}_{k}\triangleq\left[\begin{array}[]{cccc}(A-KC)^{k-1}K&\dots&(A-KC)K&K\end{array}\right]$$

$$G\triangleq {\mathcal{O}}_f{\mathcal{K}}_p.$$

$$\mathcal{T}_{s}\triangleq\left[\begin{array}[]{cccc}I_{m}&0&&0\\ CK&I_{m}&\cdots&0\\ \vdots&\vdots&&\vdots\\ CA^{s-2}K&CA^{s-3}K&\cdots&I_{m}\end{array}\right].$$

$$Y_k^{+}={\mathcal{O}}_f{\hat{x}}_k+{\mathcal{T}}_f{E}_k^{+}.$$

$${\hat{x}}_k=K{y}_{k-1}+\cdots +(A-KC{)}^{p-1}K{y}_{k-p}+(A-KC{)}^p{\hat{x}}_{k-p}.$$

$$Y_{+}=G{Y}_{-}+{\mathcal{O}}_f(A-KC{)}^p\hat{X}+{\mathcal{T}}_f{E}_{+},$$

$$\hat{G}=Y_{+}{Y}_{-}^{\ast }(Y_{-}{Y}_{-}^{\ast }{)}^{-1}.$$

$$\hat{G}=\left[\begin{array}[]{cc}\hat{U}_{1}&\hat{U}_{2}\end{array}\right]\left[\begin{array}[]{cc}\hat{\Sigma}_{1}&0\\ 0&\hat{\Sigma}_{2}\end{array}\right]\left[\begin{array}[]{c}\hat{V}_{1}^{*}\\ \hat{V}_{2}^{*}\end{array}\right],$$

$${\hat{\mathcal{O}}}_f={\hat{U}}_1{\hat{\Sigma }}_1^{1/2},{\hat{\mathcal{K}}}_p={\hat{\Sigma }}_1^{1/2}{\hat{V}}_1^{\ast }.$$

$$\hat{C}={\hat{\mathcal{O}}}_f\left(1:m,:\right),\hat{K}={\hat{\mathcal{K}}}_p\left(:,(p-1)m+1:pm\right),$$

$${\hat{\mathcal{O}}}_f^u\triangleq {\hat{\mathcal{O}}}_f\left(1:m(f-1),:\right),$$

$$\hat{A}={\left({\hat{\mathcal{O}}}_f^u\right)}^{†}{\hat{\mathcal{O}}}_p^l,$$

$${\Sigma }_E=\mathbb{E}\left[{\mathcal{T}}_p{E}_k^{-}(E_k^{-}{)}^{\ast }{\mathcal{T}}_p^{\ast }\right]={\mathcal{T}}_p\mathrm{diag}(\bar{R},\dots ,\bar{R}){\mathcal{T}}_p^{\ast }.$$

$${\sigma }_E\triangleq {\sigma }_{\min }\left({\Sigma }_E\right).$$

$${\delta }_N\triangleq {\left(2(N+p-1)m\right)}^{-{\log }^2\left(2pm\right)\log \left(2(N+p-1)m\right)}.$$

$${\Vert G-\hat{G}\Vert }_2\le$$

Full text

1 Introduction

Identifying predictive models from data has been a fundamental problem across several fields, from classical control theory to economics and modern machine learning. System identification, in particular, has a long history of studying this problem from a control theoretic perspective [1]. Identifying linear state-space models:

[TABLE]

from input-output data has been the focus of time-domain identification. In fact, some identification algorithms can not only learn the system matrices in (1) but also the Kalman filter required for state estimation [2].

Most identification methods for linear systems either follow the prediction error approach [3] or the subspace method [2, 4]. The prediction error approach is usually non-convex and directly searches over the system parameters $A,B,C,D$ by minimizing a prediction error cost. The subspace approach is a convex one; first, Hankel matrices of the system are estimated, then, the parameters are realized via steps involving singular value decomposition (SVD). Methods inspired by machine learning have also also been employed [5]. In this paper, we focus on the subspace identification approach–see [6] for an overview. The convergence properties of subspace algorithms have been studied before in [7, 8, 9, 10, 11, 12]; the analysis relies on the assumption of asymptotic stability (spectral radius $\rho(A)<1$) and is limited to asymptotic results. In [7, 8] it is shown that the identification error can decrease as fast as $\mathcal{O}\left(1/\sqrt{N}\right)$ up to logarithmic factors, as the number of output data $N$ grows to infinity. While asymptotic results have been established, a finite sample analysis of subspace algorithms remains an open problem [2]. Another open question is whether the asymptotic stability condition $\rho(A)<1$ can be relaxed to marginal stability $\rho(A)\leq 1$.

From a machine learning perspective, finite sample analysis has been a standard tool for comparing algorithms in the non-asymptotic regime. Early work in finite sample analysis of system identification can be found in [13, 14, 15]. A series of recent papers [16, 17, 18] studied the finite sample properties of system identification from a single trajectory, when the system state is fully observed ($C=I$). Finite sample results for partially observed systems ($C\neq I$), which is a more challenging problem, appeared recently in [19, 20, 21]. These papers provide a non-asymptotic convergence rate of $1/\sqrt{N}$ for the recovery of matrices $A,B,C,D$ up to a similarity transformation. The results rely on the assumption that that the system can be driven by external inputs, i.e. $B,D\neq 0$. In [19], the analysis of the classical Ho-Kalman realization algorithm was explored. In [20], it was shown that with a prefiltering step, consistency can be achieved even for marginally stable systems where $\rho\left(A\right)\leq 1$. In [21], identification of a system of unknown order is considered. Finite sample properties of system identification algorithms have also been used to design controllers [22]. The dual problem of Kalman filtering has not been studied yet in this context.

In this paper, we perform the first finite sample analysis of system (1) in the case $B,D=0$, when we have no inputs, also known as stochastic system identification (SSI) [2]. This problem is more challenging than the case $B,D\neq 0$, since the system can only be driven through noise and establishing persistence of excitation is harder. We provide the first non-asymptotic guarantees for the estimation of matrices $A,C$ as well as the Kalman filter gain of (1). Similar to [18, 19], the analysis is based on new tools from machine learning and statistics [23, 24, 25]. In summary the main contributions of this paper are:

• •

To the best of our knowledge, our paper provides the first finite sample upper bounds in the case of stochastic system identification, where we have no inputs and the system is only driven by noise. We also provide the first finite sample guarantees for the estimation error of the Kalman filter gain.

• •

We prove that the outputs of the system satisfy persistence of excitation in finite time with high probability. This result is fundamental for the analysis of most subspace identification algorithms, which use outputs as regressors.

• •

We show that we can achieve a learning rate of $\mathcal{O}(\sqrt{1/N})$ up to logarithmic factors for marginally stable systems. To the best of our knowledge, the classical subspace identification results do not offer guarantees in the case of marginal stability where $\rho(A)\leq 1$. For stable systems ($\rho(A)<1$), this rate is consistent with classical asymptotic results [7].

All proofs are included in the Appendix.

2 Problem formulation

Consider the standard state space representation (1) with $B,D=0$, where $x_{k}\in\mathbb{R}^{n}$ is the system state, $y_{k}\in\mathbb{R}^{m}$ is the output, $A\in\mathbb{R}^{n\times n}$ is the system matrix, $C\in\mathbb{R}^{m\times n}$ is the output matrix, $w_{k}\in\mathbb{R}^{n}$ is the process noise and $v_{k}\in\mathbb{R}^{m}$ is the measurement noise. The noises $w_{k},\,v_{k}$ are assumed to be i.i.d. zero mean Gaussian, with covariance matrices $Q$ and $R$ respectively, and independent of each other. The initial state $x_{0}$ is also assumed to be zero mean Gaussian, independent of the noises, with covariance $\Sigma_{0}$. Matrices $A$, $C$, $Q$, $R$, $\Sigma_{0}$ are initially unknown. However, the following assumption holds throughout the paper.

Assumption 1**.**

The order of the system $n$ is known. System $A$ is marginally stable: $\rho\left(A\right)\leq 1$, where $\rho$ denotes the spectral radius. The pair $(A,C)$ is observable, $(A,Q^{1/2})$ is controllable and $R$ is strictly positive definite. $\diamond$

We leave the case of unknown $n$ for future work 111The results of Section 4 do not depend on the order $n$ being known.. The assumption $\rho(A)\leq 1$ is more general than the stricter condition $\rho(A)<1$ found in previous works, see [7, 8, 9, 10, 11, 12]. The remaining conditions in Assumption 1 are standard for the stochastic system identification problem to be well posed and the Kalman filter to converge.

The steady-state Kalman filter of system (1) is :

[TABLE]

where $K\in\mathbb{R}^{n\times m}$ is the steady-state Kalman gain

[TABLE]

and $P$ is the positive definite solution of the Riccati equation:

[TABLE]

A byproduct of Assumption 1 is that the closed-loop matrix $A-KC$ has all the eigenvalues strictly inside the unit circle [26]. The state of the Kalman filter is the minimum mean square error prediction:

[TABLE]

We denote its covariance matrix by:

[TABLE]

The innovation error sequence $e_{k}$ has covariance

[TABLE]

Since the original errors are Gaussian i.i.d., by the orthogonality principle the innovation error sequence $e_{k}$ is also Gaussian and i.i.d. The later property is true since we also assumed that the Kalman filter is in steady-state.

Assumption 2**.**

We assume that $\Sigma_{0}=P$, so that the Kalman filter (2) has converged to its steady-state. $\diamond$

Since the Kalman filter converges exponentially fast to the steady-state gain, this assumption is reasonable in many situations; it is also standard [7, 12]. Nonetheless, we leave the case of general $\Sigma_{0}$ for future work.

In the classical stochastic subspace identification problem, the main goal is to identify the Kalman filter parameters $A,C,K$ from output samples $y_{0}\dots,y_{N}$–see for example Chapter 3 of [2]. The problem is ill-posed in general since the outputs are invariant under any similarity transformation $\bar{A}=S^{-1}AS$, $\bar{C}=CS$, $\bar{K}=S^{-1}K$. Thus, we can only estimate $A,C,K$ up to a similarity transformation.

In this paper, we will analyze the finite sample properties of a subspace identification algorithm, which is based on least squares.

Problem 1** (Finite Sample Analysis of SSI).**

Consider a finite number $N$ of output samples $y_{0},\dots,y_{N-1}$, which follow model (1) with $B,D=0$, and an algorithm $\mathcal{A}$, which returns estimates $\hat{A},\hat{C},\hat{K}$ of the true parameters. Given a confidence level $\delta$ provide upper bounds $\epsilon_{A}\left(\delta,N\right)$, $\epsilon_{C}\left(\delta,N\right)$, $\epsilon_{K}\left(\delta,N\right)$ and an invertible matrix $S$ such that with probability at least $1-\delta$:

[TABLE]

where $\left\|\cdot\right\|_{2}$ denotes the spectral norm. The bounds $\epsilon$ can also depend on the model parameters $n,A,C,R,Q$ as well as the identification algorithm used. $\diamond$

3 Subspace Identification Algorithm

The procedure of estimating the parameters $A,C,K$ is based on a least square approach, see for example [7, 12]. It involves two stages. First, we regress future outputs to past outputs to obtain a Hankel-like matrix, which is a product of an observability and a controllability matrix. Second, we perform a balanced realization step, similar to the Ho-Kalman algorithm, to obtain estimates for $A,C,K$.

Before describing the algorithm, we need some definitions. Let $p,f$, with $p,f\geq n$ be two design parameters that define the horizons of the past and the future respectively. Assume that the total number of output samples is $\bar{N}=N+p+f-1$. Then, the future outputs $Y^{+}_{k}\in\mathbb{R}^{mf}$ and past outputs $Y^{-}_{k}\in\mathbb{R}^{mp}$ at time $k\geq p$ are defined as follows:

[TABLE]

By stacking the outputs for all sample sequences, over all times $p\leq k\leq N+p-1$, we form the batch outputs:

[TABLE]

The past and future noises $E^{+}_{k},E^{-}_{k},E_{+},E_{-}$ are defined similarly. Finally, define the batch states:

[TABLE]

The (extended) observability matrix $\mathcal{O}_{k}\in\mathbb{R}^{mk\times n}$ and the reversed (extended) controllability matrix $\mathcal{K}_{k}\in\mathbb{R}^{n\times mk}$ associated to system (2) are defined as:

[TABLE]

[TABLE]

respectively. We denote the Hankel(-like) matrix $\mathcal{O}_{f}\mathcal{K}_{p}$ by:

[TABLE]

Finally, for any $s\geq 2$, define block-Toeplitz matrix:

[TABLE]

3.1 Regression for Hankel Matrix Estimation

First, we establish a linear regression between the future and past outputs. This step is common in most (stochastic) subspace identification algorithms. From (2), for every $k$:

[TABLE]

Meanwhile, from (2), the state prediction $\hat{x}_{k}$ can be expressed in terms of the past outputs:

[TABLE]

After some algebra, we derive the linear regression:

[TABLE]

where the regressors $Y_{-}$ and the residuals $E_{+}$ are independent column-wise. The term $\mathcal{O}_{k}(A-KC)^{p}\hat{X}$ introduces a bias due to the Kalman filter truncation, where we use only $p$ past outputs instead of all of them. Based on (19), we compute the least squares estimate

[TABLE]

The Hankel matrix $G$ can be interpreted as a (truncated) Kalman filter which predicts future outputs directly from past outputs, independently of the internal state-space representation [2]. In this sense, the estimate $\hat{G}$ is a “data-driven" Kalman filter. Notice that persistence of excitation of the outputs (invertibility of $Y_{-}Y^{*}_{-}$) is required in order to compute the least squares estimate $\hat{G}$.

3.2 Balanced Realization

This step determines a balanced realization of the state-space, which is only one of the possibly infinite state-space representations–see Section 6 for comparison with other subspace methods. First, we compute a rank-$n$ factorization of the full rank matrix $\hat{G}$. Let the SVD of $\hat{G}$ be:

[TABLE]

where $\hat{\Sigma}_{1}\in\mathbb{R}^{n\times n}$ contains the $n-$largest singular values. Then, a standard realization of $\mathcal{O}_{f}$, $\mathcal{K}_{p}$ is:

[TABLE]

This step assumed knowing the order $n$ of the system, see Assumption 1. In addition, matrix $\mathcal{K}_{p}$ should have full rank $n$. This is equivalent to the pair $(A,K)$ being controllable. Otherwise, $\mathcal{O}_{f}\mathcal{K}_{p}$ will have rank less than $n$ making it impossible to accurately estimate $\mathcal{O}_{f}$.

Assumption 3**.**

The pair $(A,K)$ is controllable. $\diamond$

The above assumption is standard–see for example [12].

Based on the estimated observability/controllability matrices, we can approximate the system parameters as follows:

[TABLE]

where the notation $\hat{\mathcal{O}}_{f}\left(1:m,:\right)$ means we pick the first $m$ rows and all columns. The notation for $\hat{\mathcal{K}}_{p}$ has similar interpretation. For simplicity, define

[TABLE]

which includes the $m(f-1)$ “upper" rows of matrix $\hat{\mathcal{O}}_{f}$. Similarly, we define the lower part $\hat{\mathcal{O}}^{l}_{f}$. For matrix $A$ we exploit the structure of the extended observability matrix and solve $\hat{\mathcal{O}}^{u}_{f}\hat{A}=\hat{\mathcal{O}}^{l}_{p}$ in the least squares sense by computing

[TABLE]

where $\dagger$ denotes the pseudoinverse.

Now that we have a stochastic system identification algorithm, our goal is to perform finite sample analysis, which is divided in two parts. First, in Section 4, we provide high probability upper bounds for the error $\|G-\hat{G}\|_{2}$ in the regression step. Then, in Section 5, we analyze the robustness of the balanced realization step.

4 Finite Sample Analysis of Regression

In this section, we provide the finite sample analysis of the linear regression step of the identification algorithm. We provide high-probability upper bounds for estimation error $\|G-\hat{G}\|_{2}$ of the Hankel-like matrix $G$. Before we state the main result, let us introduce some notation. Recall the definition of the innovation covariance matrix $\bar{R}$ in (6). We denote the past noises’ weighted covariance by:

[TABLE]

The least singular value of the above matrix is denoted by:

[TABLE]

In Lemma B.2 in the Appendix, we show that $\sigma_{E}\geq\sigma_{\min}\left(R\right)>0$.

Theorem 1** (Regression Step Analysis).**

Consider system (2) under the Assumptions 1, 2, 3. Let $\hat{G}$ be the estimate (20) of the subspace identification algorithm given an output trajectory $y_{0},\dots,y_{N+p+f-1}$ and let $G$ be as in (17). Fix a confidence $\delta>0$ and define:

[TABLE]

There exist $N_{0},N_{1},N_{2}$ such that if $N\geq N_{0},N_{1},N_{2}$, (see definitions (B.1), (B.5), (C.3) in the Appendix), then with probability at least $1-\delta_{N}-6\delta$:

[TABLE]

where

[TABLE]

over-approximates the condition number of $\mathbb{E}\left[Y_{-}Y_{-}^{*}\right]$ and

[TABLE]

are system-dependent constants. $\diamond$

The final result in (26) has been simplified for compactness–see (C.4) for the full expression in the Appendix.

Remark 1** (Result interpretation).**

From (19), (20) the estimation error consists of two terms:

[TABLE]

The first term in (26) corresponds to the cross-term error, while the second term corresponds to the Kalman filter truncation bias term. To obtain consistency for $\hat{G}$, we have to let the term $\left\|(A-KC)^{p}\right\|_{2}$ go to zero with $N$. Recall that the matrix $A-KC$ has spectral radius less than one, thus, the second term decreases exponentially with $p$. By selecting $p=c\log N$, for some $c$, we can force the Kalman truncation error term to decrease at least as fast as the first one, see for example [7]. In this sense, the dominant term is the first one, i.e. the cross-term. Notice, that $f$ can be kept bounded as long as it is larger than $n$. Notice that the norm $\left\|\mathcal{O}^{\dagger}_{p}\right\|$ remains bounded by $\left\|\mathcal{O}^{\dagger}_{n}\right\|$ as $p$ increases–see Lemma E.3 in Appendix. $\diamond$

Remark 2** (Statistical rates).**

For marginally stable systems ($\rho(A)\leq 1$) and $p=c\log N$, we have $\log\kappa_{N}=\mathcal{O}\left(\log N\right)$, since $\left\|\mathcal{O}_{p}\right\|_{2},\operatorname{\mathrm{tr}}\Gamma_{N}$ depend at most polynomially on $p,N$–see Corollary E.1 in the Appendix. In this case,  (26) results in a rate of:

[TABLE]

To the best of our knowledge, there have not been any bounds for the performance of subspace algorithms in the case of marginally stable systems.

In the case of asymptotically stable systems ($\rho(A)<1$), we have $\kappa_{N}=\mathcal{O}\left(p\right)$, since $\left\|\mathcal{O}_{p}\right\|_{2},\operatorname{\mathrm{tr}}\Gamma_{N},\left\|\mathcal{T}_{p}\right\|_{2}$ are now $\mathcal{O}(1)$–see Corollary E.1 in the Appendix. Hence, if $p=c\log N$, we obtain a rate of:

[TABLE]

As a result, our finite sample bound (26) is consistent with the asymptotic bound in equation (14) of [7] as $N$ grows to infinity. $\diamond$

In the absence of inputs ($B,D=0$), the noise both helps and obstructs identification. Larger noise leads to better excitation of the outputs, but also worsens the convergence of the least squares estimator. To see how our finite sample bounds capture that, observe that larger noise leads to bigger $\sigma_{E}$ but also bigger $\left\|\bar{R}\right\|_{2}$. This trade-off is captured by $\mathcal{C}_{1}$.

If $N$ is sufficiently large (condition $N\geq N_{0},N_{1}$), the outputs are guaranteed to be persistently exciting in finite time; more details can be found in Section 4.1 and the Appendix. Meanwhile, condition $N\geq N_{2}$ is not necessary; it just leads to a simplified expression for the bound of the Kalman filter truncation error–see Section 4.3 and Appendix. The definitions of $N_{0},N_{1},N_{2}$ can be found in (B.1), (B.5), (C.3). Their existence is guaranteed even if $p$ varies slowly with $N$, i.e. logarithmically.

Obtaining the bound on the error $\|G-\hat{G}\|_{2}$ in (26) of Theorem 1 requires the following three steps:

Proving persistence of excitation (PE) for the past outputs, i.e. invertibility of $Y_{-}Y_{-}^{*}$. 2. 2.

Establishing bounds for the cross-term error in (29). 3. 3.

Establishing bounds for the the Kalman truncation error in (29).

In the following subsections, we sketch the proof steps.

4.1 Persistence of Excitation in Finite Time

The next theorem shows that with high probability the past outputs and noises are persistently exciting in finite time. The result is of independent interest and is fundamental since most subspace algorithms use past outputs as regressors.

Theorem 2** (Persistence of Excitation).**

Consider the conditions of Theorem 1 and $N_{0}$, $N_{1}$ as in (B.1), (B.5). If $N\geq N_{0},N_{1}$, then with probability at least $1-\delta_{N}-2\delta$ both of the following events occur:

[TABLE]

where $\succeq$ denotes comparison in the positive semidefinite cone. Hence, with probability at least $1-\delta_{N}-2\delta$ the outputs satisfy the PE condition:

[TABLE]

where $\sigma_{E}>0$ is defined in (24). $\diamond$

The above result implies that if the past noises satisfy a PE condition, then PE for the outputs is also guaranteed; the noises are the only way to excite the system in the absence of control inputs. The prove PE, from (2), the past outputs satisfy:

[TABLE]

Thus, their correlations are

[TABLE]

We can first show PE for the noise correlations $\mathcal{T}_{p}E_{-}E^{*}_{-}\mathcal{T}^{*}_{p}$, i.e. show that the event $\mathcal{E}_{E}$ occurs with high probability when $N$ is sufficiently large (condition $N\geq N_{0}$). This behavior is due to the fact that $\mathbb{E}\left[\mathcal{T}_{p}E_{-}E^{*}_{-}\mathcal{T}^{*}_{p}\right]=N\Sigma_{E}$ and the sequence $E^{-}_{k}$ is component-wise i.i.d. To prove this step, we use Lemma C.2 from [19]–see Lemma B.1 in the Appendix.

Meanwhile, the cross terms $\hat{X}E^{*}$ are much smaller and their norm increases with a rate of at most $\mathcal{O}(\sqrt{N})$ up to logarithmic terms. This is since $\mathbb{E}\left[\hat{X}E^{*}\right]=0$ and the product $\hat{X}E^{*}$ has martingale structure (see Appendix and Theorem 3 below). Eventually, if the number of samples $N$ is large enough (condition $N\geq N_{1}$), the cross-terms will be dominated by the noise and state correlations with high probability, which establishes output PE.

4.2 Cross-term error

To bound the cross-term error, we express it as a product of $E_{+}Y^{*}_{-}\left(Y_{-}Y_{-}^{*}\right)^{-1/2}$ and $(Y_{-}Y_{-}^{*})^{-1/2}$, as in [18]. The second term of the product can be bounded by applying Theorem 2. The first term is self-normalized and has martingale structure component-wise. In particular, the product $Y_{-}E^{*}_{+}$ is equal to:

[TABLE]

where every sum above is a martingale. To bound it, we apply the next theorem, which generalizes Theorem 1 in [24] and Proposition 8.2 in [18].

Theorem 3** (Cross terms).**

Let $\left\{\mathcal{F}_{t}\right\}_{t=0}^{\infty}$ be a filtration. Let $\eta_{t}\in\mathbb{R}^{m}$, $t\geq 0$ be $\mathcal{F}_{t}$-measurable, independent of $\mathcal{F}_{t-1}$. Suppose also that $\eta_{t}$ has independent components $\eta_{t,i}$ $i=1,\dots,m$, which are $1-$sub-Gaussian:

[TABLE]

Let $X_{t}\in\mathbb{R}^{d}$, $t\geq 0$ be $\mathcal{F}_{t-1}-$measurable. Assume that $V$ is a $d\times d$ positive definite matrix. For any $t\geq 0$, define:

[TABLE]

where

[TABLE]

for some integer $r$. Then, for any $\delta>0$, with probability at least $1-\delta$, for all $t\geq 0$

[TABLE]

$\diamond$**

The above theorem along with a Markov upper bound on $Y_{-}Y_{-}^{*}$ (see Lemma B.3 in the Appendix) are used to bound $E_{+}Y^{*}_{-}\left(Y_{-}Y_{-}^{*}\right)^{-1/2}$.

4.3 Kalman truncation error

For the Kalman truncation error term, we need to bound the term $\hat{X}Y^{*}_{-}(Y_{-}Y^{*}_{-})^{-1}$, which is $\mathcal{O}\left(1\right)$. Using the identities $\mathcal{O}_{p}^{\dagger}\mathcal{O}_{p}\hat{X}=\hat{X}$, and $Y_{-}=\mathcal{O}_{p}\hat{X}+\mathcal{T}_{p}E_{-}$, we derive the following equality:

[TABLE]

From Theorem 2, we obtain $\left\|\mathcal{T}_{p}E_{-}E_{-}^{*}\mathcal{T}^{*}_{p}(Y_{-}Y^{*}_{-})^{-1}\right\|_{2}\leq 2$. The last term in (33) can be treated like the cross-term in Section 4.2, by applying Theorems 2, 3 and Lemma B.3. It decreases with a rate of $\mathcal{O}\left(1/\sqrt{N}\right)$ up to logarithmic terms, so it is much smaller than the other terms in (33). To keep the final bound simple, we select $N_{2}$ such that

[TABLE]

with high probability–see also (C.3) for the definition of $N_{2}$.

5 Robustness of Balanced Realization

In this section, we analyze the robustness of the balanced realization. In particular, we upper bound the estimation errors of matrices $A,C,K$ in terms of the estimation error $\|G-\hat{G}\|_{2}$ obtained by Theorem 1.

Assume that we knew $G$ exactly. Then, the SVD of the true $G$, would be:

[TABLE]

for some $\Sigma_{1}\in\mathbb{R}^{n\times n}$. Hence, if we knew $G$ exactly, the output of the balanced realization would be:

[TABLE]

The respective matrices $\bar{C},\bar{K},\bar{A}$ are defined similarly, based on $\bar{\mathcal{O}}_{f},\,\bar{\mathcal{K}}_{p}$, as described in Section 3. The system matrices $\bar{C},\bar{K},\bar{A}$ are equivalent to the original matrices $C,K,A$ up to a similarity transformation $\bar{C}=CS$, $\bar{K}=S^{-1}K$, $\bar{A}=S^{-1}AS$ for some invertible $S$. For simplicity, we will quantify the estimation errors in terms of the similar $\bar{A},\bar{C},\bar{K}$ instead of the original $A,C,K$.

The next result follows the steps of [19] and relies on Lemma 5.14 of [25] and Theorem 4.1 of [27]. Let $\sigma_{n}\left(\cdot\right)$ denote the $n-$th largest singular value.

Theorem 4** (Realization robustness).**

Consider the true Hankel-like matrix $G$ defined in (17) and the noisy estimate $\hat{G}$ defined in (20). Let $\hat{A},\hat{C},\hat{K},\hat{\mathcal{O}}_{f}$, $\hat{\mathcal{K}}_{p}$ be the output of the balanced realization algorithm based on $\hat{G}$. Let $\bar{A},\bar{C},\bar{K},\bar{\mathcal{O}}_{f}$, $\bar{\mathcal{K}}_{p}$ be the output of the balanced realization algorithm based on the true $G$. If $G$ has rank $n$ and the following robustness condition is satisfied:

[TABLE]

there exists an orthonormal matrix $T\in\mathbb{R}^{n\times n}$ such that:

[TABLE]

where

[TABLE]

The notation $\hat{\mathcal{O}}^{u}_{f},\bar{\mathcal{O}}^{u}_{f},$ refers to the upper part of the respective matrix (first $(f-1)m$ rows)–see Section 3.2. $\diamond$

Remark 3**.**

The result states that if the error of the regression step is small enough, then the realization is robust. The singular value $\sigma_{n}\left(G\right)$ can be quite small. Hence, the robustness condition (36) can be quite restrictive in practice. However, such a condition is a fundamental limitation of the SVD procedure. The gap between the singular values of $G$ and the singular values coming from the noise $G-\hat{G}$ should be large enough; this guarantees that the singular vectors related to small singular values of $G$ are separated from the singular vectors coming from the noise $G-\hat{G}$, which can be arbitrary. See also Wedin’s theorem [28]. Such robustness conditions have also appeared in model reduction theory [29]. $\diamond$

The term $\frac{\sqrt{\left\|G\right\|_{2}}+\sigma_{o}}{\sigma^{2}_{o}}$ which appears in the bound of $A$ is $\mathcal{O}\left(1\right)$. The value of $\sigma^{-1}_{o}$ is random since it depends on $\hat{\mathcal{O}}^{h}_{f}$. However, we could replace it by a deterministic bound. From:

[TABLE]

$\sigma_{o}$ will eventually be lower bounded by $\sigma_{n}\left(\bar{\mathcal{O}}^{h}_{f}\right)/2$ if the error $\left\|\hat{\mathcal{O}}_{f}-\bar{\mathcal{O}}_{f}T\right\|_{2}\leq\sigma_{n}\left(\bar{\mathcal{O}}^{h}_{f}\right)/2$ is small enough. The norm $\|G\|_{2}\leq\left\|\mathcal{O}_{f}\right\|_{2}\left\|\mathcal{K}_{p}\right\|_{2}$ is bounded as $p$ increases, since $A-KC$ is asymptotically stable and $f$ is finite.

Remark 4** (Total bounds).**

The final upper bounds for the estimation of the system parameters $A,C,K$, as stated in Problem 1, can be found by combining the finite sample guarantees of the regression step (Theorem 1) with the robustness analysis of the realization step (Theorem 4). All matrix estimation errors depend linearly on the Hankel matrix estimation error $\|G-\hat{G}\|_{2}$. As a result, all matrix errors have the same statistical rate as the error of $G$, i.e. their estimation error decreases at least as fast as $\mathcal{O}\left(1/\sqrt{N}\right)$ up to logarithmic factors. $\diamond$

6 Discussion and Future Work

One of the main differences between the subspace algorithm considered in this paper and other stochastic subspace identification algorithms is the SVD step. The other algorithms perform SVD on $W_{1}GW_{2}$ instead of $G$, where $W_{1},W_{2}$ are full rank weighting matrices, possibly data dependent [30, 3, 2]. From this point of view, the results of Section 4 (upper bound for $\|G-\hat{G}\|$ in Theorem 1 and persistence of excitation in Theorem 2) are fundamental for understanding the finite sample properties of other subspace identification algorithms. Here, we studied the case $W_{1}=I,W_{2}=I$, which is not the standard choice [12]. It is subject of future work to explore how the choice of $W_{1},W_{2}$ affects the realization step, especially the robustness condition of the SVD step.

We could also provide finite sample bounds for the estimation of the closed-loop matrix $A_{c}\triangleq A-KC$. This can be done in two ways. One option is to form matrix $\hat{A}-\hat{K}\hat{C}$ from the estimates $\hat{A},\hat{C},\hat{K}$. Alternatively, we could estimate $\hat{A_{c}}$ directly from $\hat{\mathcal{K}}_{p}$ in the same way that we estimated $\hat{A}$. However, we do not have finite sample guarantees for the stability of the closed-loop estimates $\hat{A}-\hat{K}\hat{C}$, $\hat{A_{c}}$. It would be interesting to address this in future work.

Another direction for future work is repeating the analysis when the Kalman filter has not reached steady-state, i.e. relax Assumption 2. Finally, in this work, we only considered upper bounds. It would be interesting to study lower bounds as well to evaluate the tightness of our upper bounds. In any case, from lower bounds for fully observed systems [17], the factor of $1/\sqrt{N}$ is tight.

Appendix A Proof of Theorem 3

Let us first state a result which follows from the arguments of Chapter 4 of [23]. See also Proposition 8.1 of [18].

Proposition A.1** ([23]).**

Consider a matrix $M\in\mathbb{R}^{m\times n}$. Let $S^{n-1}$ denote the unit sphere. Then for any $\epsilon>0$:

[TABLE]

$\diamond$**

Second, we state Theorem 1 of [24], which upper bounds self-normalized martingales.

Theorem A.1** (Theorem 1 in [24]).**

Let $\left\{\mathcal{F}_{t}\right\}_{t=0}^{\infty}$ be a filtration. Let $\eta_{t}\in\mathbb{R}$, $t\geq 0$ be real-valued, $\mathcal{F}_{t}$-measurable, and conditionally $1-$sub-Gaussian:

[TABLE]

Let $X_{t}\in\mathbb{R}^{d}$, $t\geq 0$ be vector valued and $\mathcal{F}_{t-1}-$measurable. Assume that $V$ is a $d\times d$ positive definite matrix. For any $t\geq 0$, define:

[TABLE]

Then, for any $\delta>0$, with probability at least $1-\delta$, for all $t\geq 0$

[TABLE]

$\diamond$**

Finally, we state a standard linear algebra result. We include the proof for completeness.

Lemma A.1** (Block Matrix Norm).**

Assume

[TABLE]

for matrices of appropriate dimensions. Then:

[TABLE]

$\diamond$**

Proof.

Consider a vector $x$ such that $Mx$ is defined. Then, from triangle inequality, the definition of matrix norm and Cauchy-Schwartz:

[TABLE]

where we used that $\left\|x\right\|^{2}_{2}=\left\|x_{1}\right\|_{2}^{2}+\dots+\left\|x_{r}\right\|^{2}$. ∎

Now, we can prove Theorem 3. Notice that:

[TABLE]

We can analyze each component $\bar{V}_{t}^{-1/2}\sum_{s=1}^{t}X_{s}\eta^{*}_{k+i}$ separately and apply a union bound afterwards, since by Lemma A.1:

[TABLE]

Now, fix a $0\leq i<r$ and let

[TABLE]

Consider an arbitrary element of the unit sphere $\xi\in S^{m-1}$. The scalar $\eta^{*}_{s+i}\xi$ is conditionally $1$-sub-Gaussian and satisfies the conditions of Theorem A.1. Thus, with probability at least $1-\frac{\delta}{r5^{m}}$:

[TABLE]

Now, we apply Proposition A.1 for $\epsilon=1/2$:

[TABLE]

Finally, by (A.3) and a union bound over all components:

[TABLE]

∎

Appendix B Persistence of Excitation

The main focus of this section is the proof of Theorem 2, which provides finite sample guarantees for the PE of the past noises and the past outputs. We also include upper bounds for the sample correlations of the past outputs and the states $\hat{x}_{k}$. Finally, we provide the definition of $N_{0},N_{1}$ that we hided in the main theorem statements.

The following result shows that with high probability, the past noises are persistently exciting. It follows from Lemma C.2 of [19], which in turn is based on results for random circulant matrices [31].

Lemma B.1** (Noise PE).**

Consider the conditions of Theorem 2 and the definition of $\delta_{N}$ in (25):

[TABLE]

There exists a universal constant $c$ (independent of system and algorithm parameters) such that if $N\geq 2cpm\log{1/\delta_{N}}$, then with probability at least $1-\delta_{N}$ the event:

[TABLE]

occurs, where $\Sigma_{E}$ is defined in (23). $\diamond$

Proof.

We can rewrite $E^{-}_{k}=\operatorname{diag}(\bar{R}^{1/2},\dots,\bar{R}^{1/2})U^{-}_{k}$, where $U^{-}_{k}$ is defined similarly to $E^{-}_{k}$ but has components with unit covariance. Now from Lemma C.2 of [19] applied on $U^{-}_{k}$ we obtain that with probability at least $1-\delta_{N}$:

[TABLE]

Multiplying by $\mathcal{T}_{p}\operatorname{diag}(\bar{R}^{1/2},\dots)$ from the left and $\operatorname{diag}(\bar{R}^{1/2},\dots)\mathcal{T}_{p}^{*}$ from the right gives the desired result. ∎

From the above lemma it follows that $N$ should be large enough to guarantee PE for the noises. In particular, $N$ should be larger than $N_{0}$, where

[TABLE]

Such a $N_{0}$ exists since the term $2cpm\log^{2}\left(2pm\right)\log^{2}\left(2(N+p-1)m\right)$ depends logarithmically on $N$.

To guarantee PE for the noises, we also need to show that the smallest singular value $\sigma_{E}=\sigma_{\min}\left(\Sigma_{E}\right)$ is positive.

Lemma B.2**.**

Let $\Sigma_{E}$ be as in (23). Then:

[TABLE]

$\diamond$**

Proof.

The past errors $\mathcal{T}_{p}E^{-}_{k}$ can be rewritten as:

[TABLE]

where $V^{-}_{k}$, $W^{-}_{k}$ are defined similarly to $Y^{-}_{k}$ and consist of the past measurement and process noises respectively. Matrix $\mathcal{T}$ is a block Toeplitz matrix; we omit its analytical expression. By independence of the measurement noise, the process noise and $x_{k-p}-\hat{x}_{k-p}$:

[TABLE]

This in turn implies $\sigma_{E}\geq\sigma_{\min}\left(R\right)$. By Assumption 1, $R\succ 0$ and $\sigma_{\min}\left(R\right)>0$. ∎

Lemma B.3** (Markov upper bounds).**

Consider system (2) and recall the definition of $\Gamma_{k}$, $\Sigma_{E}$ in (5), (23). We have the following upper bounds:

[TABLE]

Proof.

We only show (B.3). The proof of (B.2) is similar. From Markov’s inequality we obtain:

[TABLE]

What remains is to bound the expectation in the right-hand side. Notice that $\|Y^{-}_{k}\left(Y^{-}_{k}\right)^{*}\|_{2}=\operatorname{\mathrm{tr}}Y^{-}_{k}\left(Y^{-}_{k}\right)^{*}$, since $Y^{-}_{k}$ has unit rank. Hence from the triangle inequality:

[TABLE]

The right hand side of the above inequality is

[TABLE]

The result follows from $\operatorname{\mathrm{tr}}\left(\mathcal{O}_{p}\Gamma_{k-p}\mathcal{O}^{*}_{p}\right)\leq\left\|\mathcal{O}_{p}\right\|^{2}_{2}\operatorname{\mathrm{tr}}\Gamma_{k-p}$ along with

[TABLE]

since the sequence $\Gamma_{k}$ is monotone (see the following Lemma). ∎

The following result is standard and we include it for completeness.

Lemma B.4** (Monotonicity of $\Gamma_{k}$).**

Consider system (2), under Assumption 1. The sequence $\Gamma_{k}=\mathbb{E}\left[\hat{x}_{k}\hat{x}^{*}_{k}\right]$ is monotone: $\Gamma_{k}\succeq\Gamma_{k-1}$.

Proof.

Notice that since $\hat{x}_{0}=0$, we have $\Gamma_{0}=0$. Define $\bar{Q}=K\bar{R}K^{*}$. By the orthogonality principle, $\hat{x}_{k}$ and $e_{k}$ are uncorrelated. Hence

[TABLE]

For $k=1$ we obtain:

[TABLE]

But the operator $\mathcal{L}$ is monotone, which implies:

[TABLE]

The result $\Gamma_{k}\succeq\Gamma_{k-1}$ follows by induction. ∎

Proof of Theorem 2

Some arguments are similar to Section 9 of [18].

Step 1: Noise PE. Under the condition $N\geq N_{0}$, from Lemma B.1 the event $\mathcal{E}_{E}$ occurs with probability at least $1-\delta_{N}$.

Step 2: Cross terms are small. Next, we show that the cross terms $\hat{X}E^{*}_{-}$ are small. As in Lemma B.1, express $E_{-}=\operatorname{diag}(\bar{R}^{1/2},\dots,\bar{R}^{1/2})H_{-}$, where $H_{-}$ is defined similarly to $E_{-}$ but has unit variance components. Define:

[TABLE]

Notice that

[TABLE]

Hence, by Theorem 3 applied to $\bar{V}_{N}$, $\hat{X}H^{*}_{-}$ the event:

[TABLE]

occurs with probability at least $1-\delta$.

Next we upper bound term $\bar{V}_{N}$. From Lemma B.3, the event:

[TABLE]

occurs with probability at least $1-\delta$. This implies that:

[TABLE]

Combining the two events $\mathcal{E}_{1},\mathcal{E}_{2}$ and by a union bound, with probability at least $1-2\delta$ the event:

[TABLE]

occurs, where

[TABLE]

As a consequence, if $u\in\mathbb{R}^{mp}$, $\left\|u\right\|_{2}=1$ is an arbitrary unit vector:

[TABLE]

Step 3: Output PE Consider an arbitrary unit vector $u\in\mathbb{R}^{mp}$, $\left\|u\right\|_{2}=1$. Consider the events $\mathcal{E}_{E}$ and $\mathcal{E}_{XE}$ from steps 1,2. With probability $1-\delta_{N}-2\delta$, since $N\geq N_{0}$ the event $\mathcal{E}_{E}\cap\mathcal{E}_{XE}$ occurs. It remains to show that on $\mathcal{E}_{E}\cap\mathcal{E}_{XE}$ the outputs satisfy PE for sufficiently large $N$. Define

[TABLE]

From (32), (B.4) for $N\geq N_{0}$ on $\mathcal{E}_{E}\cap\mathcal{E}_{XE}$:

[TABLE]

with $\beta\geq{\sigma_{E}}/{2}$. Now let $N_{1}$ be such that:

[TABLE]

Since $\mathcal{C}_{XE}$ grows at most logarithmically with $N$, $N_{1}$ always exists. Now, since $N\geq N_{1}$ and $\beta\geq\sigma_{E}/2$:

[TABLE]

The above inequality follows from the following lemma. ∎

Lemma B.5** (Minimum of function).**

Let $\beta\geq b>0$, for some $b>0$ and consider the function:

[TABLE]

If $\gamma\leq 1,\frac{b}{2}$, then

[TABLE]

$\diamond$**

Proof.

By elementary calculus:

[TABLE]

Thus, if $\gamma\leq 1,b/2$, we have $f\left(\alpha,\beta\right)\geq 0$. ∎

Appendix C Proof of Theorem 1

Step 1: Since $N\geq N_{0},N_{1}$, from Theorem 2, with probability at least $1-\delta_{N}-2\delta$, the event $\mathcal{E}_{Y}\cap\mathcal{E}_{E}$ occurs.

Step 2: Next, we analyze the cross term. Define:

[TABLE]

Notice that on the event $\mathcal{E}_{Y}\cap\mathcal{E}_{E}$, we have that $\bar{Z}_{N}\preceq 2Z_{N}$ since $Z_{N}\succeq N\frac{\sigma_{E}}{4}I_{mp}$. Hence,

[TABLE]

Now the proof continues as in the case of cross-terms in the proof of Theorem 2. We apply Theorem 3 to $\bar{Z}_{N},S_{N}$ and use Lemma B.3 to upper bound $\bar{Z}_{N}$. Then, conditioned on $\mathcal{E}_{Y}\cap\mathcal{E}_{E}$, with probability $1-2\delta$:

[TABLE]

where

[TABLE]

Next, we bound $\left\|Z^{-1/2}_{N}\right\|$ separately on the event $\mathcal{E}_{Y}\cap\mathcal{E}_{E}$ by $\frac{2}{\sqrt{N\sigma_{E}}}$.

Finally, conditioned on $\mathcal{E}_{E}\cap\mathcal{E}_{Y}$, with probability at least $1-2\delta$:

[TABLE]

Step 3: We bound the Kalman truncation term. Recall that:

[TABLE]

On the event $\mathcal{E}_{E}\cap\mathcal{E}_{Y}$, we have:

[TABLE]

since $Y_{-}Y^{*}_{-}\succeq\frac{1}{2}\mathcal{T}_{p}E_{-}E_{-}^{*}\mathcal{T}^{*}_{p}$. Hence we obtain:

[TABLE]

From the discussion in Section 4, we only need to find $N_{2}$ such that for $N\geq N_{2}$ with high probability:

[TABLE]

Define

[TABLE]

Notice that on the event $\mathcal{E}_{E}\cap\mathcal{E}_{Y}$, we have $Y_{-}Y^{*}_{-}\succeq\frac{1}{2}\bar{B}_{N}$, which implies:

[TABLE]

Now we can treat the right-hand side in the same way as the cross-term above. By an application of Theorem 3 and Lemma B.3, we obtain that conditioned on $\mathcal{E}_{Y}\cap\mathcal{E}_{E}$, with probability $1-2\delta$:

[TABLE]

where

[TABLE]

Thus, we define:

[TABLE]

Such an $N_{2}$ exists since $\mathcal{C}_{N}$ grows at most logarithmically with $N$.

Step 4: Final expression From the previous step and a union bound, for $N\geq N_{0},N_{1},N_{2}$ with probability at least $1-\delta_{N}-6\delta$:

[TABLE]

Step 5: Simplification of final expression To simplify the final expression, we use

[TABLE]

and

[TABLE]

since $p\geq n+1\geq 2$.

Appendix D Proof of Theorem 4

The proof follows the steps of [19].

Step 1: Bounds for observability/controllability matrix

Denote the rank $n$ approximation of $\hat{G}$ by

[TABLE]

By definition, $\hat{G}_{n}$ is the matrix which minimizes $\left\|\hat{G}-M\right\|_{2}$, among all rank $n$ matrices $M$. Thus, by optimality:

[TABLE]

since $G$ has also rank $n$. As a result, we have:

[TABLE]

From (D.1) and the robustness condition (36) we have:

[TABLE]

Hence, we can now apply Theorem 5.14 of [25], which states that there exists an orthonormal matrix $T$ such that:

[TABLE]

where $\left\|\cdot\right\|_{F}$ denotes the Frobenius norm.

Since matrices $G,\hat{G}_{n}$ have rank-$n$, the sum $G-\hat{G}_{n}$ has rank at most $2n$. Thus, we can bound the Frobenius norm $\left\|G-\hat{G}_{n}\right\|_{F}$ in terms of spectral norm:

[TABLE]

where the second inequality follows from (D.1). For simplicity, we also use $\frac{2}{\sqrt{2}-1}\leq 5$. Thus, from (D.3) and the above inequalities:

[TABLE]

Finally, since the spectral norm is always smaller than the Frobenius one:

[TABLE]

As a corollary,

[TABLE]

Step 2: Bounds for system parameters

Bounds for $C,K$

Since $\hat{C}-\bar{C}T$ is a sub-matrix of $\hat{\mathcal{O}}_{f}-\bar{\mathcal{O}}_{f}T$, we immediately obtain:

[TABLE]

Similarly, for $K$ we have:

[TABLE]

Bounds for $A$

For simplicity, denote $\hat{M}\triangleq\hat{\mathcal{O}}^{h}_{f}$, $\bar{M}\triangleq\bar{\mathcal{O}}^{h}_{f}$ and $\hat{N}\triangleq\hat{\mathcal{O}}^{l}_{f}$, $\bar{N}\triangleq\bar{\mathcal{O}}^{l}_{f}$. Based on this notation:

[TABLE]

After some algebraic manipulations:

[TABLE]

First, notice that $\left\|\bar{N}\right\|_{2}\leq\left\|\bar{\mathcal{O}}_{f}\right\|_{2}=\sqrt{\left\|G\right\|}$, where the inequality follows from the fact that $\bar{N}$ is a submatrix of $\bar{\mathcal{O}}_{f}$; equality follows from the definition of $\bar{\mathcal{O}}_{f}=U_{1}\Sigma_{1}^{1/2}$. Second, $\left\|\hat{M}^{\dagger}\right\|_{2}=\frac{1}{\sigma_{n}\left(\hat{M}\right)}\leq\frac{1}{\sigma_{o}}$ and third

[TABLE]

since $\hat{N}-\bar{N}T$ is a submatrix of $\hat{\mathcal{O}}_{f}-\bar{\mathcal{O}}_{f}T$. Finally, from Theorem 4.1 of [28]

[TABLE]

Combining all previous bounds, we obtain

[TABLE]

Appendix E Bounds for system matrices

In this section, we formally prove bounds for $\mathcal{O}_{k},\mathcal{T}_{k},\Gamma_{k}$, which we implicitly used in Section 4.

First, we prove a standard result for the norm of (block) Toeplitz matrices.

Lemma E.1** (Toeplitz norm).**

Let $M\in\mathbb{R}^{m_{1}n\times m_{2}n}$, for some integers $n,m_{1},m_{2}$ be an (upper) block triangular Toeplitz matrix:

[TABLE]

where $M_{i}\in\mathbb{R}^{m_{1}\times m_{2}}$, $i=1,\dots,n$. Then:

[TABLE]

Proof.

A standard technique in the analysis of Toeplitz matrices is to write them in terms of linear combinations of powers of companion matrices (see [32], equation (0.9.7)). We can write $M$ as:

[TABLE]

where $\otimes$ is the Kronecker product, $I_{n}$ is the identity matrix of dimension $n$, and $J_{n}$ is the companion matrix:

[TABLE]

But we have $\left\|D\otimes F\right\|_{2}=\left\|D\right\|_{2}\left\|F\right\|_{2}$ for any matrices $D,F$ (see Theorem 4.2.15 in [33]). Also the companion matrix has norm one $\left\|J_{n}\right\|_{2}=1$ since $J_{n}J_{n}^{*}=\operatorname{diag}(1,\dots,1,0)$. The result follows from the triangle inequality. ∎

Next, we provide a bound for the sum of powers of $A$.

Lemma E.2**.**

Consider the series $S_{t}=\sum_{i=0}^{t}\left\|A^{i}\right\|_{2}$. We have the following two cases:

• •

If the system is asymptotically stable $\rho(A)<1$, then $\left\|S_{t}\right\|_{2}=\mathcal{O}\left(1\right)$

• •

If the system is marginally stable ($\rho(A)=1$), then $\left\|S_{t}\right\|_{2}=\mathcal{O}\left(t^{\kappa}\right)$, where $\kappa$ is the largest Jordan block of $A$ corresponding to a unit circle eigenvalue $\left|\lambda\right|=1$.

Proof.

Proof of first part. By Gelfand’s formula [32], for every $\epsilon>0$, there exists a $i_{0}=i_{0}(\epsilon)$ such that $\left\|A^{i}\right\|\leq(\rho(A)+\epsilon)^{i}$, for all $i\geq i_{0}$. Just pick $\epsilon$ such that $\rho\left(A\right)+\epsilon<1$. Then,

[TABLE]

Proof of second part. Assume that $A$ is equal to a $n\times n$ Jordan block corresponding to $\lambda=1$. The proof for the other cases is similar. Then we have that:

[TABLE]

By Lemma E.1, we obtain:

[TABLE]

where the second inequality is classical, see Exercise 0.0.5 in [23].

Finally, we have:

[TABLE]

∎

Corollary E.1**.**

The norms $\left\|\mathcal{O}_{k}\right\|_{2}$, $\left\|\mathcal{T}_{k}\right\|_{2}$, $\left\|\Gamma_{k}\right\|_{2}$ depend at most polynomially in $k$ if $\rho\left(A\right)\leq 1$ (they are $\mathcal{O}\left(\text{poly}(k)\right)$). If the system is asymptotically stable they are upper bounded for all $k$ (they are $\mathcal{O}\left(1\right)$).

Proof.

Observability matrix: Consider the sum $S_{t}=\sum_{i=0}^{t}\left\|A^{i}\right\|_{2}$ of the previous lemma. We have:

[TABLE]

The result follows by applying the previous lemma.

Block Toeplitz matrix Consider the sum $S_{t}=\sum_{i=0}^{t}\left\|A^{i}\right\|_{2}$ of the previous lemma. By Lemma E.1 we have:

[TABLE]

The result follows by applying the previous lemma.

Covariance matrix We have

[TABLE]

The result follows by using a similar argument as in the previous lemma for $\sum_{i=0}^{k-1}\left\|A^{i}\right\|^{2}$. ∎

Lemma E.3** (Lest non-zero singular values of $G$, $\mathcal{O}_{p}$ are increasing).**

The $n-$th singular value of $\mathcal{O}_{p}$ is increasing with $p$:

[TABLE]

The same is true for the $n-$th singular value of $G$:

[TABLE]

Proof.

We only prove the first part. The other proof is similar. Notice that $\mathcal{O}_{p_{1}}$ can be rewritten as:

[TABLE]

for some matrix $M$. Hence

[TABLE]

Thus $\sigma^{2}_{n}\left(\mathcal{O}_{p_{1}}\right)\geq\sigma^{2}_{n}\left(\mathcal{O}_{p_{2}}\right)$. ∎