Learning Sparse Dynamical Systems from a Single Sample Trajectory

TL;DR

This paper proposes a Lasso-based method for identifying sparse linear systems from a single trajectory, providing finite-time guarantees and improved sample complexity bounds, with applications demonstrated in power systems.

Contribution

Introduces a novel sparse system identification approach with sharp finite-time guarantees and polynomial sample complexity bounds, extending existing methods.

Findings

High-probability correct sparsity pattern recovery

Sample complexity scales polynomially with nonzero elements

Effective bounds on estimation error influenced by system properties

Abstract

This paper addresses the problem of identifying sparse linear time-invariant (LTI) systems from a single sample trajectory generated by the system dynamics. We introduce a Lasso-like estimator for the parameters of the system, taking into account their sparse nature. Assuming that the system is stable, or that it is equipped with an initial stabilizing controller, we provide sharp finite-time guarantees on the accurate recovery of both the sparsity structure and the parameter values of the system. In particular, we show that the proposed estimator can correctly identify the sparsity pattern of the system matrices with high probability, provided that the length of the sample trajectory exceeds a threshold. Furthermore, we show that this threshold scales polynomially in the number of nonzero elements in the system matrices, but logarithmically in the system dimensions --- this improves on existing sample complexity bounds for the sparse system identification problem. We further extend these results to obtain sharp bounds on the $\ell_{\infty}$-norm of the estimation error and show how different properties of the system---such as its stability level and \textit{mutual incoherency}---affect this bound. Finally, an extensive case study on power systems is presented to illustrate the performance of the proposed estimation method.

Equations

$$x(t+1)=Ax(t)+Bu(t)+w(t)$$

$$\displaystyle(\hat{A},\hat{B})=\arg\min_{A,B}\$$

$$+\lambda (\Vert A{\Vert }_1+\Vert B{\Vert }_1).$$

$$(\hat{A},\hat{B})=\arg \underset{A,B}{\min }$$

$$+\lambda (\Vert A{\Vert }_1+\Vert B{\Vert }_1).$$

$$Y=\left[\begin{array}{c}x(1{)}^{\top }\\ ⋮\\ x(T{)}^{\top }\end{array}\right],X=\left[\begin{array}{cc}x(0{)}^{\top }& u(0{)}^{\top }\\ ⋮& ⋮\\ x(T-1{)}^{\top }& u(T-1{)}^{\top }\end{array}\right],W=\left[\begin{array}{c}w(0{)}^{\top }\\ ⋮\\ w(T-1{)}^{\top }\end{array}\right].$$

$$\hat{\Psi }=\arg \underset{\Psi }{\min }\frac{1}{2T}\Vert Y-X\Psi {\Vert }_F^2+\lambda \Vert \Psi {\Vert }_1$$

$$\mathbb{E}\{\exp (tx)\}\le \exp \left(\frac{b^2{t}^2}{2}\right)$$

$$x(0)=\underset{\tilde{T}\to \infty }{\lim }\sum _{\tau =-\tilde{T}}^{-1}(A+B{K}_0{)}^{-\tau -1}(w(\tau )+Bv(\tau ))$$

$$(A+B{K}_0)Q^{\ast }(A+B{K}_0{)}^{\top }-Q^{\ast }+{\sigma }_w^{2}I+{\sigma }_v^{2}B{B}^{\top }=0$$

$$M^{\ast }=\left[\begin{array}{cc}{Q}^{\ast }& Q^{\ast }{K}_0^T\\ K_0{Q}^{\ast }& K_0{Q}^{\ast }{K}_0^T+{\sigma }_v^{2}I\end{array}\right]$$

$$\underset{1\le j\le n}{\max }\left\{\underset{i\in {\mathcal{A}}_j^c}{\max }\left\{{\Vert M_{i{\mathcal{A}}_j}^{\ast }(M_{{\mathcal{A}}_j{\mathcal{A}}_j}^{\ast }{)}^{-1}\Vert }_1\right\}\right\}\le 1-\gamma$$

$$\underset{1\le j\le n}{\min }{\lambda }_{\min }(M_{{\mathcal{A}}_j{\mathcal{A}}_j}^{\ast })\ge C_{\min }$$

$$\underset{1\le j\le n}{\max }{\left|\left|\left|(M_{{\mathcal{A}}_j{\mathcal{A}}_j}^{\ast }{)}^{-1}\right|\right|\right|}_{\infty }\le D_{\max }$$

$$\underset{1\le j\le n}{\min }\left\{\underset{i\in {\mathcal{A}}_j}{\max }\left\{|{\Psi }_{ij}^{\ast }|\right\}\right\}\ge {\Psi }_{\min }$$

$$\lambda =c_1\cdot \frac{C}{1-\rho }\cdot \frac{{\eta }^2}{\gamma }\sqrt{\frac{\log ((n+m)/\delta )}{T}}$$

$$T\ge c_2\cdot \frac{C^4}{(1-\rho {)}^4}\cdot \frac{D_{\max }^2}{{\gamma }^2{C}_{\min }^2{\Psi }_{\min }^2}\cdot k^2\log ((n+m)/\delta ),$$

$$\Vert \hat{\Psi }-{\Psi }^{\ast }{\Vert }_{\infty }\le c_3\cdot \frac{C}{1-\rho }\cdot \frac{D_{\max }{\eta }^2}{\gamma }\sqrt{\frac{\log ((n+m)/\delta )}{T}}$$

$$\displaystyle\min_{\Psi}\ \$$

$$\displaystyle\mathrm{s.t.}\ \$$

$$\mu (A)=\underset{1\le i<j\le t_2}{\max }\frac{|A_{:,i}^{\top }{A}_{:,j}|}{\Vert A_{:,i}{\Vert }_2\Vert A_{:,j}{\Vert }_2}$$

$$\Vert {\Psi }_{:,j}^{\ast }{\Vert }_0<\frac{1}{2}\left(1+\frac{1}{\mu (X)}\right)$$

$$\mu (X)$$

$$M_i{\ddot{\theta }}_i+D_i{\dot{\theta }}_i=P_{M_i}-P_{E_i}$$

$$P_{E_i}=\sum _{j\in {\mathcal{N}}_i}{B}_{ij}({\theta }_i-{\theta }_j)$$

$$x_i(t+1)=\left(A_{ii}{x}_i(t)+\sum _{j\in {\mathcal{N}}_i}{A}_{ij}{x}_j(t)\right)+B_{ii}{u}_i(t)+w_i(t)$$

$$A_{ii}=\left[\begin{array}{cc}1& dt\\ -\frac{{\sum }_{j\in {\mathcal{N}}_i}{B}_{ij}}{M_i}dt& 1-\frac{D_i}{M_i}dt\end{array}\right],A_{ij}=\left[\begin{array}{cc}0& 0\\ \frac{B_{ij}}{M_i}dt& 0\end{array}\right],B_{ii}=\left[\begin{array}{c}0\\ 1\end{array}\right]$$

$$L({\Psi }_{:,j})=\Vert Y-X{\Psi }_{:,j}{\Vert }_2^2$$

$${\hat{\Psi }}_{:,j}=\arg \min \frac{1}{2T}L({\Psi }_{:,j})+\lambda \Vert {\Psi }_{:,j}{\Vert }_1$$

$$\hat{\Psi }=\left[\begin{array}{cccc}{\hat{\Psi }}_{:,1}& {\hat{\Psi }}_{:,2}& \cdots & {\hat{\Psi }}_{:,n}\end{array}\right]$$

$$G=-\nabla L({\Psi }_{:,j}){|}_{{\Psi }_{:,j}={\Psi }_{:,j}^{\ast }}=\frac{1}{T}{X}^T{W}_{:,j},$$

$$M={\nabla }^{2}L({\Psi }_{:,j}){|}_{{\Psi }_{:,j}={\Psi }_{:,j}^{\ast }}=\frac{1}{T}{X}^{T}X$$

$$\Vert G{\Vert }_{\infty }\le \frac{\lambda \gamma }{3},$$

Full text

1 Introduction

Modern cyber-physical systems, such as power grids, autonomous transportation systems, and distributed computing and sensing networks, are characterized by being large scale, spatially distributed, and by having complex ever changing dynamics and interconnected topologies. The distributed optimal control literature addresses set-point tracking and regulation in the distributed setting by assuming known dynamics with a sparse interconnections. Indeed, the underlying sparsity structure of a distributed system is aggressively (and necessarily) exploited, with foundational results showing that both tractability [1] and scalability [2, 3, 4, 5] in controller synthesis are only possible when the underlying dynamical system is suitably sparse. However, in this large-scale, dynamic, and complex setting, it is unclear how to obtain the necessary models of the dynamical systems. To address this issue, we use data-driven approaches to identify both the interconnected topology and the dynamic behavior of these systems for which first-principle modeling becomes either intractable or impractical for such large-scale dynamic systems.

This then raises a more fundamental question: how can data-driven methods be appropriately integrated into safety-critical control loops? This question has been addressed in the context of learning [6, 7], and control of a small-scale and dense unknown systems, e.g., a single autonomous vehicle or robot [8, 9, 10, 11, 12]. These works make clear that if a learned model is to be integrated into a safety-critical control loop, then it is essential that the uncertainty associated with the learned model be explicitly quantified. This way, the learned model and the uncertainty bounds can be integrated with tools from robust control to provide strong guarantees of system performance and stability. This paper takes a first step towards extending these results to the large-scale distributed setting by providing a sample efficient and computationally tractable algorithm for the identification of sparse dynamical systems, as well as providing sharp estimates on the corresponding model uncertainty.

Main contributions: We show that large-scale sparse system models can be identified with a complexity scaling quadratically with the number of nonzero elements in the underlying dynamical system—for systems composed of a large number of subsystems that only interact with a small number of local neighbors, this computational saving can be significant. We further provide sharp bounds on the corresponding model uncertainty, paving the way for the use of these models in safety-critical control loops. Finally, in contrast to previous work, we show that such models can be extracted from a single trajectory of the system. In the context of large-scale systems, the system resets needed by methods relying on independent trajectories become prohibitively more expensive and impractical—indeed contrast resetting a robotic arm and a power distribution network, and the increase in difficulty becomes apparent. Note that we defer a detailed comparison of our results to prior work to Section 3.

Paper organization: In Section 2, we formally define the sparse system-identification task that we consider, and introduce our Lasso-like estimator based on a single system trajectory. Section 3 presents our main result, and compares and contrasts it with existing results in the literature. We also show that some of the technical assumptions that we make are necessary for a well-posed problem. We provide an empirical study of our method on a power system in Section 4. We end with conclusions in Section 5. The proofs are deferred to the appendix to streamline the presentation.

Notation: For a matrix $M$, the symbols ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|M\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$, ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|M\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\infty}$, $\|M\|_{F}$, $\|M\|_{1}$, and $\|M\|_{\infty}$ are used to denote its induced spectral, induced infinity, Frobenius, element-wise $\ell_{1}/\ell_{1}$, and element-wise $\ell_{\infty}/\ell_{\infty}$ norms, respectively. Furthermore, $\|M\|_{0}$ refers to the number of nonzero elements in $M$. The symbols $M_{:j}$ and $M_{j:}$ indicate the $j^{\text{th}}$ column and row of $M$, respectively. For a set $\mathcal{I}$, the symbol $|\mathcal{I}|$ denotes its cardinality. Given the index sets $\mathcal{U}$ and $\mathcal{V}$, define $M_{\mathcal{U}\mathcal{V}}$ as the $|\mathcal{U}|\times|\mathcal{V}|$ submatrix of $M$ obtained by removing the rows and columns with indices not belonging to $\mathcal{U}$ and $\mathcal{V}$. The symbols $c$ and $c_{i}$ play the role of universal constants throughout the paper. $\mathbb{E}\left\{x\right\}$ denotes the expected value of a random variable $x$. For an event $\mathcal{E}$, the notation $\mathbb{P}(\mathcal{E})$ refers to its probability of occurrence. The notation $x_{n}\overset{a.s.}{\rightarrow}x$ means that a sequence of random variables $x_{n}$ converges to $x$ almost surely.

2 Problem Statement

Consider the linear time-invariant (LTI) system

[TABLE]

where $A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times m}$ are the unknown state and input matrices, respectively. Furthermore, $x(t)\in\mathbb{R}^{n}$, $u(t)\in\mathbb{R}^{m}$, and $w(t)\in\mathbb{R}^{n}$ are the respective state, input, and disturbance vectors at time $t$.

The goal of this work is to estimate the underlying parameters of the dynamics, based on a limited number of sample trajectories, i.e., a sequence $\{(x^{(i)}(\tau),u^{(i)}(\tau))\}_{\tau=0}^{T}$ with $i=1,2,...,d$, where $d$ is the number of available sample trajectories and $T$ is the length of each sample trajectory. To simplify the notations, the superscript $i$ is dropped from the sample trajectories when $d=1$.

This paper is concerned with the identification of high dimensional but sparse system matrices $(A,B)$. Such high-dimensional sparse parameters arise in the context of large-scale distributed and multi-agent systems, where dynamic coupling arises due to local interactions between subsystems–it is this local interaction structure that results in correspondingly sparse system matrices. Examples of such systems include power grids, intelligent transportation systems, and distributed computation and sensing networks.

We now compare and contrast two approaches to collecting sample trajectories from a dynamical system (1):

Fixed $\mathbf{d}$ and variable $\mathbf{T}$: In this method, the number of sample trajectories $d$ is set to a fixed value (e.g., $d=1$) and instead, a sufficiently long time horizon (also referred to as learning time) $T$ is chosen to collect enough information about the dynamics. This approach is most suitable when the open-loop system is stable, or if a stabilizing controller is provided—note that this assumption of stability is necessary, as even a simple least-squares estimator may not be consistent if the system has unstable modes [6]. From a practical perspective, system instability may also impose limits on how large the learning time can be in order to ensure system safety, thereby restricting the amount of data that can be collected.

Fixed $\mathbf{T}$ and variable $\mathbf{d}$: In this approach, the learning time $T$ is fixed and instead, the number of sample trajectories is chosen to be sufficiently large. Notice that this method is not dependent on the system stability. However, one needs to reset the initial state of the system at the beginning of each sample trajectory, which may not be possible in practice, especially in the case of large-scale systems.

This work focuses on sparse system identification using a single trajectory, where it is assumed that the system is either stable, or equipped with an initial stabilizing controller, and our goal is to both identify the supports of the sparse system matrices $(A,B)$ and estimate their values, using a single sample trajectory. As mentioned in [8], in many applications, the existence of an initial stabilizing controller for the unknown system (1) is not restrictive. In fact, [9] and [12] respectively introduce offline and adaptive procedures for designing such an initial stabilizing controller.

Indeed, one can cast the sparse system identification task as a supervised learning problem, where the goal is to fit the linear model (1)—parameterized by $(A,B)$—to a limited number of measurements $\{(x(\tau),u(\tau))\}_{\tau=0}^{T}$. Motivated by this observation, one can consider the following $M$-estimator:

[TABLE]

where the first term corresponds to the maximum likelihood estimation of $(A,B)$ when the disturbance noise has a zero-mean Gaussian distribution, and the second term has the role of promoting sparsity in the estimated $(\hat{A},\hat{B})$.

Before proceeding, it is essential to note that there are fundamental limits on the performance of the introduced estimator. In particular, the above optimization problem may not have a unique solution for any length of the sample trajectory. To see this, suppose that $u(t)=K_{0}x(t)$ and $K_{0}$ is equal to the identity matrix. Then, the above optimization problem reduces to

[TABLE]

It is easy to see that, given any optimal solution $(\hat{A},\hat{B})$ to the above optimization, $(\tilde{A},\tilde{B})=(\alpha\hat{A},(1-\alpha)\hat{B})$ is also optimal for any $0\leq\alpha\leq 1$. To break this symmetry and to guarantee the identifiability of the parameters, it is essential to inject an input noise to the system at every time $t$. In particular, we assume that $u(t)=K_{0}x(t)+v(t)$, where $v(t)$ is a random vector with a user-defined distribution. As another example, if $A$ is stable and $K_{0}=0$, the need to introduce noise in the input is inevitable in order to identify the matrix $B$.

To further analyze the properties of the above estimator, one can write (1) in a compact form. Let $\Psi^{*}=\begin{bmatrix}A&B\end{bmatrix}^{\top}$ denote the true parameters of the system. Furthermore, define

[TABLE]

The system identification problem is then reduced to estimating the unknown parameter $\Psi^{*}$ given the design matrix $X$, and the observation matrix $Y$ that is corrupted with the noise matrix $W$. We can therefore rewrite optimization problem (2) compactly as

[TABLE]

which corresponds to the so-called Lasso estimator, initially popularized in statistics and machine learning to estimate the support parameter values of a sparse linear model [13]. The non-asymptotic properties of this estimator have been widely studied in the literature [14, 15, 16], all highlighting its sub-linear sample complexity under suitable technical conditions. In particular, they show that under the so-called mutual incoherency of the design matrix and the sparsity of the unknown parameters, the minimum number of observations for the accurate estimation of the Lasso scales logarithmically in the dimension of $\Psi$. Motivated by these results, one may speculate that the proposed estimator (2) benefits from a similar logarithmic sample complexity. However, the validity of the derived non-asymptotic estimation error bounds on the Lasso is contingent upon a number of assumptions on the independence between the design matrix $X$ and the noise matrix $W$ [14, 17]; such assumptions do not necessarily hold in the sparse system identification problem, partly due to the dependency between the states, the inputs and the disturbance noise. The problematic nature of this dependency becomes more evident by noting that the Lasso may not be consistent when the design and noise matrices are dependent [18].

This lack of independence in the design and noise matrices of the sparse system identification problem has been the main roadblock in deriving similar sub-linear sample complexity bounds for the sparse system identification problem and it leaves the following question unanswered:

Is the estimator (2) consistent, and if so, what is its sample complexity?

3 Main Results

Despite the fact that in general, the Lasso may not be a consistent estimator when the design and noise matrices are dependent, we exploit the underlying structure of the system identification problem to control this dependency and provide an affirmative answer to the posed question. In other words, we show that not only is the proposed estimator (2) consistent, but that it also enjoys a logarithmic sample complexity in the state and input dimensions, under appropriate conditions. To this goal, we first provide a number of definitions.

Definition 1**.**

A zero-mean (centered) random variable $x$ is sub-Gaussian with parameter $b$ if its moment generating function satisfies

[TABLE]

for every $t$.

For a centered sub-Gaussian random variable $x$ with parameter $b$, one can easily verify that $\mathbb{P}(|x|>t)\leq 2\exp\left(\frac{t^{2}}{2b^{2}}\right)$. The most commonly known examples of such random variables are Gaussian, Bernoulli, and any bounded random variable.

Definition 2**.**

Given a sub-Gaussian random variable $x$, its sub-Gaussian norm, denoted by $\|x\|_{\psi}$ is defined as the smallest $r>0$ such that the inequality $\mathbb{E}\{{x^{2}}/{r^{2}}\}\leq 2$ is satisfied.

It is well-known that the above two definitions are closely related. In particular, it can be verified that $\frac{1}{\sqrt{5}}b\leq\|x\|_{\psi}\leq\sqrt{\frac{8}{3}}b$ for a sub-Gaussian random variable with parameter $b$.111This is a standard result; see [19] and [20] for a simple proof. For a random vector $x$ with sub-Gaussian elements, $\|x\|_{\psi}$ is defined as $\max_{i}\{\|x_{i}\|_{\psi}\}$.

As mentioned before, we assume that the dynamical system is equipped with an initial static and stabilizing state-feedback controller $K_{0}$. More specifically, we assume that at any given time $t$, the input $u(t)$ is equal to $K_{0}x(t)+v(t)$, where $v(t)$ is a user-defined input noise with independent and centered sub-Gaussian elements whose non-zero variance is upper bounded by $\sigma_{v}^{2}$ (for stable systems, $K_{0}$ can be set to zero). Similarly, we assume that the disturbance noise at every time $t$ is a random vector with independent and centered sub-Gaussian elements whose variance is upper bounded by $\sigma_{u}^{2}$. Further, let $\eta>0$ be the smallest positive constant such that $\max\{\|w(t)\|_{\psi},\|v(t)\|_{\psi}\}\leq\eta$; such a constant is guaranteed to exist as $w$ and $v$ are assumed to be centered sub-Gaussian random variables.

Remark 1**.**

Most of the existing results on the sample complexity of the system identification problem assume a centered Gaussian distribution for the input noise [7, 21, 9]. Despite having desirable finite-time properties, these types of Gaussian inputs may jeopardize the safety of the dynamical system due to their unbounded range. Accordingly, in many control systems, the input is constrained to have a limited power. These types of constraints can be translated into $\ell_{\infty}$ or $\ell_{2}$ bounds on the input signal. Due to the fact that such bounded random signals are sub-Gaussian, our results are readily applied to system identification problems with input constraints.

Notice that for LTI systems, the uniform asymptotic stability of the closed-loop system is equivalent to its exponential stability. In other words, an LTI system is uniformly asymptotically stable if and only if there exist constants $C\geq 1$ and $0<\rho<1$ such that ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(A+BK_{0})^{\tau}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq C\rho^{\tau}$ for every time $\tau$. Without loss of generality, let $C\geq 1$ and $0\leq\rho<1$ be the smallest constants such that ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(A+BK_{0})^{\tau}B\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq C\rho^{\tau}$, ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|K_{0}(A+BK_{0})^{\tau}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq C\rho^{\tau}$ and ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|K_{0}(A+BK_{0})^{\tau}B\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq C\rho^{\tau}$ for every time $\tau$. Note that the existence of such $C\geq 1$ and $0<\rho<1$ is guaranteed due to the exponential stability of the closed-loop system.

Furthermore, we assume that the initial state $x(0)$ rests at its stationary distribution or, equivalently, the following equality holds:

[TABLE]

Note that, for exponentially stable systems, the state converges to its stationary distribution exponentially fast and therefore, the stationarity of $x(0)$ is a reasonable assumption. Furthermore, using the above equality, it is easy to see that $x(0)$ is a random vector whose elements are (dependent) centered sub-Gaussian random variables with bounded parameters. Moreover, one can verify that its covariance $\mathbb{E}\{x(0)x(0)^{\top}\}=Q^{*}$ satisfies the following Lyapunov equation:

[TABLE]

Accordingly, $Q^{*}$ can be used to derive the covariance matrix $M^{*}$ for the random vector $\begin{bmatrix}x(0)^{\top}&(K_{0}x(0)+v(0))^{\top}\end{bmatrix}^{\top}$:

[TABLE]

Define $\mathcal{A}_{j}=\{i:\Psi^{*}_{ij}\not=0\}$ and let $\mathcal{A}^{c}_{j}$ refer to its complement. Denote $k$ as the maximum number of nonzero elements in any column of ${\Psi^{*}}$.

Assumption 1**.**

The following inequalities are satisfied

• A1

(Mutual incoherence)

[TABLE]

• A2

(Bounded eigenvalue)

[TABLE]

• A3

(Bounded infinity norm)

[TABLE]

• A4

(Nonzero gap)

[TABLE]

for some constants $0<\gamma<1$, $1\geq C_{\min}>0$, $D_{\max}\geq 1$ and $1\geq\Psi_{\min}>0$.

Next, we present the main result of the paper.

Theorem 1**.**

Assume that $k\geq 2$ and

[TABLE]

where $c_{1}$ and $c_{2}$ are universal constants. Then, the following statements hold with probability of at least $1-\delta$:

(Correct sparsity recovery) (4) has a unique solution and recovers the true sparsity pattern of $\Psi^{*}$.

• 2.

($\ell_{\infty}$-norm error) We have

[TABLE]

where $c_{3}$ is a universal constant.

Remark 2**.**

As mentioned before, the injection of a random input noise is essential to guarantee the identifiability of the parameters. This is also reflected in the above theorem: in order to guarantee a finite sample complexity for the proposed estimator, it is crucial to have $C_{\min}>0$, which is only possible if $\sigma_{v}>0$.

A number of observations can be made based on Theorem 1. First, it implies that if $\gamma$, $C$, $D_{\max}$, $C_{\min}$, $\Psi_{\min}$, and $\rho$ do not scale with the system dimension, then $T=\Omega(k^{2}\log(n+m))$ is enough to guarantee the correct sparsity recovery and a small estimation error. Notice that for sparse systems, this quantity can be much smaller than the system dimension. Second, the sample complexity of the proposed estimator depends on $\frac{C}{1-\rho}$, which is a measure of the system stability. In particular, for highly stable systems, $\frac{C}{1-\rho}$ is small, resulting in an improved accuracy of the proposed estimator with smaller $T$. In contrast, when the system is close to its stability margin, $\frac{C}{1-\rho}$ will grow which negatively affects the estimation error as well as the lower bound on $T$. Another intuitive interpretation of $\frac{C}{1-\rho}$ is that it measures the amount of dependency between the states at different times: for highly stable systems where $\rho$ is small, $(x(t),u(t))$ is only weakly dependent on $(x(\tau),u(\tau))$ for $\tau=0,\dots,t-1$, thereby facilitating the estimation of the unknown parameters. We finally mention that this dependency is in contrast with the recent discoveries on the sample complexity of the least-squares estimator, which support the favorable effect of a large $\rho$ on the accuracy of the estimator [22]. We leave investigating whether this seemingly contradictory observation is an artifact of our methodology (e.g., mixing the initial state to the stationary distribution), or is fundamental to the sparse system identification problem, to future work.

Remark 3**.**

In order to further enhance the accuracy of the proposed estimator, one can perform a least-squares estimation restricted to the nonzero elements of the estimated parameter, after obtaining its sparsity pattern via the proposed method. Although, theoretically, this post-model-selection estimation method may not improve the estimation error rate, it will incur less bias [23]. We will show in our simulations that the effect of this post-processing step can be significant in the accuracy of the estimation.

3.1 Comparison to prior art

As mentioned before, another line of work focuses on unstructured system identification, where either the learning time $T$ or the number of sample trajectories $d$ is allowed to grow. In [9], the authors consider the sample complexity of the system identification problem with multiple sample trajectories via least-squares, where it is shown that the proposed estimator incurs a small error, provided that $d=\Omega(n+m)$. Revisiting (20) reveals that the proposed method outperforms the sample complexity of ordinary least-squares when $k$ is significantly smaller than $n+m$, i.e., exploiting prior knowledge of the system sparsity leads to a reduction in sample complexity. In [6, 22, 11, 12], the authors consider unstructured system identification from a single sample trajectory under different assumptions on system stability and/or the initial state of the system. However, similar to [9], none of these works take advantage of the underlying sparsity structures of the system matrices. As a result, they cannot correctly estimate the sparsity structure of $(A,B)$ and suffer from poor dependencies on the system dimensions in the large-scale and structure setting.

Subsequently, a Lasso-type estimator is proposed in [21] to further exploit the underlying sparsity pattern of $(A,B)$ with $d$ sample trajectories, each with a zero initial state. In particular, it is shown that $d=\Omega\left(\frac{\kappa({\Sigma})^{2}}{\gamma^{2}\Psi_{\min}^{2}}k\log(n+m)\right)$ is enough to ensure the correct sparsity recovery and a small estimation error with high probability, where $\kappa({\Sigma})$ is the condition number of the finite-time controllability matrix of the system. Comparing this quantity with (20), one can observe that the former has a better dependency on $k$. However, $\kappa({\Sigma})$ is highly dependent on the learning time $T$. In fact, it is easy to show that for unstable systems, $\kappa({\Sigma})$ may grow exponentially fast with respect to $T$. On the other hand, (20) is free of such dependency and instead, it is in terms of the stationary distributions of the state and input vectors.

Moreover, our work is a major extension to the results of [7], where the authors address a similar sparse system identification problem with a single sample trajectory. First, unlike the presented results, [7] only considers autonomous systems, i.e., systems (1) with $B$=0. Second, [7] only ensures the correct sparsity recovery of the true parameters. In contrast, we extend these results to obtain non-asymptotic bounds on the estimation error. As demonstrated in [9, 8], having these bounds is essential for the design of near-optimal and robustly stabilizing controllers. Third, [7] requires that the closed-loop system be contractive with respect to the spectral norm, i.e., that $|||(A+BK_{0})|||<1$, whereas we only require system stability. Notice that the former condition is much stronger, as in practice, stable systems are often not contractive in spectral norm. Finally, the validity of the non-asymptotic bounds introduced in [7] heavily relies on the Gaussian nature of the disturbance and input noises. As an extension to this result, our proposed method targets a larger class of uncertainties for the disturbance and input noises, thereby allowing for norm bounded disturbance and input signals.

3.2 Mutual incoherency

In this subsection, we analyze the mutual incoherence condition on the steady-state covariance matrix $M^{*}$. In particular, we explain why this assumption is not an artifact of the proposed method, but that it rather stems from a fundamental limitation of any sparsity-promoting technique for the system identification problem. We show that similar mutual incoherence assumptions are indeed necessary to recover the correct sparsity of system parameters by using a class of oracle estimators.

We assume that the oracle estimator can measure the disturbance matrix $W$ and that it can work with sample trajectories of an arbitrary length. With these assumptions, the oracle estimator solves the following optimization problem to estimate the parameters of the system:

[TABLE]

Clearly, this oracle estimator cannot be used in practice since 1) the disturbance matrix $W$ is unknown, 2) the learning time $T$ is finite, and 3) the corresponding optimization problem is non-convex and NP-hard in its worst case. Setting aside these restrictions for now, there are fundamental limits on the consistency of this estimator. To explain this, we introduce the mutual-coherence metric for a matrix (note the difference between this definition and Assumption A1). For a given matrix $A\in\mathbb{R}^{t_{1}\times t_{2}}$, its mutual-coherence $\mu(A)$ is defined as

[TABLE]

In other words, $\mu(A)$ measures the maximum correlation between distinct columns of $A$. Reminiscent of the classical results in the compressive sensing literature, it is well-known that the optimal solution $\Psi^{*}$ of (9) is unique if the following identifiability condition

[TABLE]

holds for $j=1,2,...,n$ (see, e.g., Theorem 2.5 in [24]). Furthermore, this bound is tight, implying that there exists an instance of the problem for which the violation of $\|\Psi^{*}_{:,j}\|_{0}<\frac{1}{2}\left(1+\frac{1}{\mu(X)}\right)$ for some $j$ results in the non-uniqueness of the optimal solution. On the other hand, according to Lemma 3 (to be introduced later) and the Borel-Cantelli lemma, $\frac{1}{T}X^{\top}X$ converges to $M^{*}$ almost surely, as $T\rightarrow\infty$. This implies that

[TABLE]

The above analysis reveals that the off-diagonal entries of $M^{*}$ play a crucial role in the identifiability of the true parameters: as these elements become smaller relative to the diagonal entries, the oracle estimator can correctly identify the structure of $\Psi$ for a wider range of sparsity levels. Similarly, our proposed mutual incoherence assumption is expected to be satisfied when the off-diagonals of $M^{*}$ have small magnitudes, relative to the diagonal entries. This implies that Assumption A1 is a natural condition to impose in order to ensure the correct sparsity recovery of $\Psi$. Furthermore, in practice, $M^{*}$ will be close to a diagonally dominant matrix with exponentially decaying off-diagonal entries, provided that the matrices $A$, $B$, and $K_{0}$ have sparse structures [25].

4 Numerical Experiments

As a case study, we consider the frequency control problem for power systems, where the goal is to control the governing frequency of the entire network, based on the so-called swing equations. Assume that there exist $N_{g}$ generators in the system. It is common to describe the per-unit swing equations using the well-known direct current (DC) approximation:

[TABLE]

where $\theta_{i}$ is the voltage angle at generator $i$, $P_{M_{i}}$ is the mechanical power input at generator $i$, and $P_{E_{i}}$ denotes the active power injection at the bus connected to generator $i$. Furthermore, $M_{i}$ and $D_{i}$ are the inertia and damping coefficients at generator $i$, respectively. Under the DC approximation, the relationship between active power injection and voltage is defined as follows:

[TABLE]

where $n$ is the number of generators in the network, $\mathcal{N}_{i}$ collects the neighbors of generator $i$, and $B_{ij}$ is the susceptance of the line $(i,j)$. After discretization with the sampling time $dt$, the system of swing equations is reduced to the following dynamical system:

[TABLE]

where $x_{i}=\begin{bmatrix}\theta_{i}&\dot{\theta_{i}}\end{bmatrix}^{\top}$, $u_{i}(t)=P_{M_{i}}$, and

[TABLE]

The goal is to identify the underlying dynamical system based on a single sample trajectory consisting of a sequence of mechanical power inputs and their effects on the angles and frequencies of different generators. To assess the performance of the proposed method, we generate several instances of the problem according to the following rules:

the generators are connected via a randomly generated tree with a maximum degree of $10$.

• -

the parameters $B_{ij}$, $M_{i}$, $D_{i}$ are uniformly chosen from $[0.5,1]$, $[1,2]$, $[0.5,1.5]$, respectively.

Furthermore, the sampling time $dt$ is set to $0.1$. We assume that the disturbance noise has a zero-mean Gaussian distribution with covariance $0.01I_{2\times 2}$. Notice that the magnitude of the noise is comparable to those of the nonzero elements in $A$ and $B$. Furthermore, the mechanical input is set to $u_{i}(t)=-0.1(\theta_{i}+\dot{\theta_{i}})+v_{i}(t)$, where $v_{i}(t)$ is a randomly generated input noise, distributed according to a zero-mean Gaussian distribution with variance $0.05$. Notice that the first term in the input signal is used to ensure the closed-loop stability.

The reported results are for a serial implementation in MATLAB R2017b, and the function lasso is used to solve (2). It is worthwhile to note that the running time can be further reduced via parallelization; this is trivially possible due to the decomposable nature of the problem. The mismatch error is defined as the total number of false positives and false negatives in the sparsity pattern of the estimated parameters $(\hat{A},\hat{B})$. Furthermore, relative learning time (RLT) is defined as the learning time normalized by the dimension of the system, and relative mismatch error (RME) is used to denote the mismatch error normalized by the total number of elements in $A$ and $B$. In all of our experiments, the regularization coefficient $\lambda$ is set to $\lambda=\sqrt{\frac{0.03\log(n+m)}{T}}$. Note that this value does not require any additional fine-tuning and is at most a constant factor away from (6).

Figure 1a illustrates the mismatch error (averaged over 10 different trials) with respect to the learning time $T$ and for different number of generators $N_{g}$ that are chosen from $\{100,200,400,800\}$. These correspond to the total system dimensions of $\{300,600,1200,2400\}$. Note that the largest instance has more than $3.84$ million unknown parameters. Not surprisingly, the learning time needed to achieve a small mismatch error increases as the dimension of the system grows. Conversely, a smaller value for RLT is needed to achieve infinitesimal RME for larger systems. In particular, when $N_{g}$ is equal to $100$, $200$, $400$, and $800$, the minimum RLT to guarantee RME $\leq 0.1\%$ is equal to $3.83$, $1.42$, $0.50$, and $0.16$, respectively.

As mentioned before, the accuracy of the proposed estimator can be improved by additionally applying the least-squares over the nonzero elements of $(\hat{A},\hat{B})$. Figure 1b illustrates the normalized 2-norm estimation error of this approach (abbreviated as LASSO+LS), compared to the proposed method without any post-processing step (abbreviated as LASSO), and the least-squares estimator (abbreviated as LS) when $N_{g}$ is set to $200$. It can be observed that both LASSO+LS and LS significantly outperform LS; in fact, LS is not even well-defined if the learning time is strictly less than the system dimensions. Furthermore, on average, the estimation error for LASSO+LS is $1.91$ times smaller than that of LASSO.

Finally, only 32 out of 360 generated instances did not satisfy the proposed mutual incoherence condition. However, this violation did not have a significant effect on the accuracy of the proposed estimator. To further investigate the frequency of the instances that satisfy this condition, we plot the histogram of the mutual incoherence parameter $\gamma$ for 2000 randomly generated instances with fixed $N_{g}=200$. It can be seen in Figure 1c that the mutual incoherence condition is violated only for $5.15\%$ of the instances.

5 Conclusions

The problem of sparse system identification of linear time-invariant (LTI) systems is considered in this work, where the goal is to estimate the sparse structure of the system matrices based on a single sample trajectory of the dynamics. A Lasso-type estimator is introduced to identify the parameters of the system, while promoting their sparsity via a $\ell_{1}$-regularization technique. By carefully examining the underlying properties of the system—such as its stability and mutual incoherency—we provide non-asymptotic bounds on the accuracy of the proposed estimator. In particular, we show that it correctly identifies the sparsity structure of the system matrices and enjoys a sharp upper bound on its estimation error, provided that the learning time exceeds a threshold. We further show that this threshold scales polynomially in the number of nonzero elements but logarithmically in the system dimensions.

Appendix A Proof of Theorem 1

In this section, we present the sketch of the proof for the main theorem. Define

[TABLE]

and

[TABLE]

for every $j\in\{1,2,...,n\}$. It is easy to verify that

[TABLE]

Furthermore, the Gradient and Hessian of $L(\cdot)$ are equal to

[TABLE]

Note that $G$ can be different for every $j$. However, we keep this dependency implicit in the notations to streamline the presentation. The following Lemma is at the core of our subsequent analysis:

Lemma 1** (Proposition 4.1 [7]).**

Suppose that the following conditions are satisfied:

[TABLE]

Then, (11) recovers the true sparsity pattern of $\Psi^{*}_{:,j}$.

The first step in proving Theorem 1 is to verify that the conditions of Lemma 1 hold with high probability. To this goal, first we write $x(t)$ and $u(t)$ in terms of $x(0)$, $w(\tau)$ and $v(\tau)$ for $\tau=0,1,\dots,t$:

[TABLE]

Instead of initiating the system at $x(0)$ with the stationary distribution, we will start at the time $-T_{0}$, with a modified initial state $x(-T_{0})=w(-T_{0}-1)+Bv(-T_{0}-1)$, where $w(-T_{0}-1)$ and $v(-T_{0}-1)$ have the same distributions as the disturbance and input noises, respectively. Since the system is stable, by taking $T_{0}\rightarrow\infty$ and invoking the Continuous Mapping Theorem, the matrices

[TABLE]

and

[TABLE]

converge in distribution to the same matrices when the system is initialized at a state with the stationary distribution. Therefore, without loss of generality, we will focus on the former. Based on this observation, one can write

[TABLE]

This implies that the elements in $G$ and $M$ can be written as quadratic functions of the disturbance and input noises in the form of $G_{i}=z^{\top}R_{G}z$ and $M_{ij}=z^{\top}R_{M}z$, where $z\in\mathbb{R}^{(n+m)(t+T_{0}+1)}$ is a random vector, defined as

[TABLE]

The following theorem will be used in our analysis to provide concentration bounds on $G$ and $M$.

Theorem 2** (Hanson-Wright inequality [26]).**

Let $x=\begin{bmatrix}x_{1}&x_{2}&\dots&x_{n}\end{bmatrix}$ be a random vector with independent zero-mean sub-Gaussian elements. Given a square and symmetric matrix $P$, the following inequality holds

[TABLE]

for every $t\geq 0$, where $c$ is a universal constant.

For a symmetric matrix $P$, we have $\|P\|^{2}_{F}=\sum_{k=1}^{n}\lambda^{2}_{k}$. Therefore, the above theorem implies that, for a sub-Gaussian random vector $z$ with independent elements, we have

[TABLE]

provided that $t\leq\left(\frac{\sum_{k}\lambda_{k}^{2}}{\max_{k}|\lambda_{k}|}\right)\|z\|^{2}_{\psi}$. The assumptions of Lemma 1 can be seen to hold directly as a consequence of the following two lemmas:

Lemma 2**.**

Let $i\in\{1,2,...,n+m\}$ and suppose that $\epsilon<\frac{3C\eta^{2}}{1-\rho}$. Then, there exists a universal constant $c_{4}$ such that

[TABLE]

Proof.

See Appendix B.1. ∎

Lemma 3**.**

Let $i,j\in\{1,2,...,n+m\}$ and suppose that $\epsilon\leq\frac{4C^{2}\eta^{2}}{(1-\rho)^{2}}$. Then, there exists a universal constant ${c}_{5}$ such that

[TABLE]

Proof.

See Appendix B.2. ∎

The following proposition shows that for a fixed column $j$, the proposed estimator (11) correctly recovers the sparsity pattern with high probability.

Proposition 1**.**

Assume that $k\geq 2$ and the following conditions are satisfied:

[TABLE]

for universal constants $c_{6},c_{7}\geq 0$. Then, (11) recovers the true sparsity pattern of $\Psi^{*}_{:,j}$ with probability of at least $1-\delta$.

Proof.

The Lemmas 2 and 3 can be used to prove statement. The details are provided in Appendix B.3. ∎

The next lemma provides a deterministic upper bound on the estimation error in terms of the deviations of $M$ and $G$ from their mean.

Lemma 4**.**

Assume that

[TABLE]

and (11) recovers the correct sparsity pattern of $\Psi^{*}_{:,j}$. Then, the following inequality holds for $E=\hat{\Psi}_{:,j}-\Psi^{*}_{:,j}$:

[TABLE]

Proof.

See Appendix B.4. ∎

The next lemma shows that the condition of Proposition 4 holds with high probability, provided that $T$ is large enough.

Proposition 2**.**

Assume that

[TABLE]

for some universal constant $c_{5}\geq 0$. Then, the following inequality holds with probability of at least $1-\delta$

[TABLE]

Proof.

Notice that $|\mathcal{A}_{j}|\leq k$. One can verify that

[TABLE]

provided that $\frac{\epsilon}{k}\leq\frac{4C^{2}\eta^{2}}{(1-\rho)^{2}}$. Setting $\epsilon=\frac{\min\{1,2\eta^{2}\}}{2D_{\max}}$ and recalling that $D_{\max},C\geq 1$, one can verify that $\frac{\epsilon}{k}\leq\frac{4C^{2}\eta^{2}}{(1-\rho)^{2}}$ is satisfied. Furthermore, by choosing $c_{8}=\frac{16}{c_{5}}$, one can certify that (16) is enough to ensure that the right hand side of the above inequality is upper bounded by $\delta$, thereby completing the proof. ∎

Proof of Theorem 1: First note that (4) can be decomposed into $n$ disjoint sub-problems over different columns of $\Psi$, each in the form of (11). Consider the following choices for $\lambda$ and $T$:

[TABLE]

where $c_{4}$, $c_{5}$, $c_{6}$, $c_{7}$, and $c_{6}$ are introduced in Lemmas 2, 3, and Propositions 1, 2. Based on the Proposition 1 and the above choices for $\lambda$ and $T$, (11) recovers the sparsity pattern of $\Psi^{*}_{:,j}$ for a given column index $j$ with probability of at least $1-\delta$. Furthermore, based on Proposition 2, the lower bound on $T$ guarantees that the inequality

[TABLE]

holds with probability of at least $1-\delta$. This, together with Proposition 4 results in

[TABLE]

with probability of at least $1-2\delta$. Now, it suffices to obtain concentration bounds for different terms of the above inequality. Based on (18) and Lemma 2, one can write

[TABLE]

This implies that, with the following choices

[TABLE]

for any $\zeta_{1}>1,\zeta_{2}>2$ that satisfy

[TABLE]

we have

[TABLE]

Note that the last term on the right hand side is due to a simple union bound on the events that (21) holds and (11) recovers the correct sparsity pattern of $\Psi^{*}_{:,j}$. Now, upon defining

[TABLE]

the inequalities in (27) are satisfied, provided that $T\geq\max\{\frac{1}{c_{4}},\frac{2}{c_{5}}\}\cdot\log(4k/\delta)$. Furthermore, combining (29) and (30) with (A) results in

[TABLE]

After plugging (29) and (30) into (26) and (25), the above inequality is reduced to

[TABLE]

with probability of at least $1-3\delta$. Due to (20), one can write

[TABLE]

Therefore,

[TABLE]

with probability of at least $1-3\delta$. Now, to conclude the proof, it suffices to perform a union bound on different columns of the solution with indices $1\leq j\leq n$. This results in

[TABLE]

with probability of at least $1-3n\delta$. Replacing $\delta$ with $\frac{\delta}{3n}$ in the above inequality concludes the proof.∎

Appendix B Proof of Auxiliary Lemmas

B.1 Proof of Lemma 2

To prove this lemma, we first introduce some notations. Define the matrix

[TABLE]

where $X(\tau)$ is a matrix valued time-dependent signal. Furthermore, define the symmetrized matrix $\tilde{R}_{1}(\cdot)=\left({R}_{1}(\cdot)+{R}_{1}(\cdot)^{T}\right)/2$. Finally, for a matrix $N$, define $[N]_{i\rightarrow j}$ as a matrix with the same size as $H$ and with all rows equal to zero except for the $j^{th}$ row which is equal to the $i^{th}$ row of $N$.

Lemma 5**.**

Let $\lambda_{k}$ be the $k^{th}$ eigenvalue of the matrix $R_{G}$ defined as

[TABLE]

Then, the following relations hold

[TABLE]

Proof.

Notice that

[TABLE]

Similar to the proof of Lemma A.3 in [7], one can verify that

[TABLE]

This completes the proof of the second statement. Finally, it is easy to see that the rank of $R_{G}$ is upper bounded by $2T$. This, together with the bound on the maximum eigenvalue completes the proof of the third statement. ∎

Define the matrix $P_{ji}\in\mathbb{R}^{n(T+T_{0}+1)\times m(T+T_{0}+1)}$ as

[TABLE]

where $E_{ji}\in\mathbb{R}^{n\times m}$ is a 0-1 matrix with 1 at its $(j,i)^{th}$ entry and 0 otherwise.

Lemma 6**.**

Let $\lambda_{k}$ be the $k^{th}$ eigenvalue of the matrix $\tilde{R}_{G}$ defined as

[TABLE]

Then, the following relations hold

[TABLE]

Proof.

The proof of the first statement follows directly from Lemma 5. Furthermore, it is easy to verify that the rank of $\tilde{R}_{G}$ is upper bounded by $4T$. This, together with the upper bound on the maximum eigenvalue completes the proof of the third statement. ∎

Proof of Lemma 2: One can easily verify that

if $i\in\{1,2,\dots,n\}$, then $G_{i}=\frac{1}{T}X_{:,i}^{T}W_{:,j}=\frac{1}{T}z^{T}R_{G}z$ where $z\in\mathbb{R}^{(n+m)(T+T_{0}+1)}$ is a random vector with independent zero-mean sub-Gaussian elements and $\|z\|_{\psi}\leq 1$.

• -

if $i\in\{n+1,\dots,n+m\}$, then $G_{i}=\frac{1}{T}X_{:,i}^{T}W_{:,j}=\frac{1}{T}z^{T}\tilde{R}_{G}z$ where $z\in\mathbb{R}^{(n+m)(T+T_{0}+1)}$ is a random vector with independent zero-mean sub-Gaussian elements and $\|z\|_{\psi}\leq 1$.

Furthermore, note that the diagonal entries of both $R_{G}$ and $\tilde{R}_{G}$ are zero and hence, $\mathbb{E}\left\{\frac{1}{T}z^{T}R_{G}z\right\}=\mathbb{E}\left\{\frac{1}{T}z^{T}\tilde{R}_{G}z\right\}=0$. This, together with Hanson-Wright inequality and Lemmas 5 and 6 completes the proof.∎

B.2 Proof of Lemma 3

Define the matrix

[TABLE]

and

[TABLE]

where the matrix $P_{j}\in\mathbb{R}^{T\times m(T+T_{0}+1)}$ has the form

[TABLE]

and $e_{j}\in\mathbb{R}^{1\times m}$ with $1$ at its $j^{th}$ entry and 0 otherwise. These notations will be used in the subsequent lemma.

Lemma 7**.**

Let $\{k_{1},k_{2},k_{3},k_{4}\}\in\{1,2,3,4\}^{4}$, where $k_{1}\not=k_{4}$ and $k_{2}\not=k_{3}$. Furthermore, let $\lambda_{k}$ be the $k^{th}$ eigenvalue of the following matrix

[TABLE]

Then, the following relations hold

[TABLE]

Proof.

To show the validity of the first statement, one can write

[TABLE]

Furthermore, similar to the proof of Lemma A.4 in [7], one can verify that

[TABLE]

Combining this with the above inequality completes the proof of the first statement. Finally, note that $R_{M}(k_{1},k_{2},k_{3},k_{4})$ can be written as

[TABLE]

which implies that its rank is upper bounded by $2T$. This, together with the upper bound on the maximum eigenvalue completes the proof. ∎

Lemma 8**.**

We have $\mathbb{E}(M)=M^{*}$.

Proof.

Define

[TABLE]

The theorem can be proven by showing

[TABLE]

In what follows, we show the validity of the first equality. The other equalities can be proven in a similar manner. We have

[TABLE]

Furthermore, notice that $x(0)$ has a stationary distribution and hence, $\mathbb{E}(x(0)x(0)^{T})=Q^{*}$. Furthermore,

[TABLE]

where the second inequality is due to (LABEL:lyap). Similarly, one can show that $\mathbb{E}(x(\tau)x(\tau)^{T})=Q^{*}$ for every $\tau\in\{2,3,\dots,T-1\}$ and hence,

[TABLE]

This completes the proof. ∎

Proof of Lemma 3: Due to Lemma 8 and upon taking $T_{0}\rightarrow\infty$, we have

[TABLE]

and hence, it suffices to obtain a bound for $\mathbb{P}\{|M_{ij}-\mathbb{E}(M_{ij})|>\epsilon\}$. We should consider four cases:

If $i,j\in\{1,2,\dots,n\}$, then $M_{ij}=\frac{1}{T}z^{T}R_{M}(1,2,1,2)z$, where $z\in\mathbb{R}^{(n+m)(T+T_{0}+1)}$ is a random vector with independent zero-mean sub-Gaussian elements and $\|z\|_{\psi}\leq 1$.

• -

If $i\in\{1,2,\dots,n\}$ and $j\in\{n+1,n+2,\dots,n+m\}$, then $M_{ij}=\frac{1}{T}z^{T}R_{M}(1,2,3,4)z$, where $z\in\mathbb{R}^{(n+m)(T+T_{0}+1)}$ is a random vector with independent zero-mean sub-Gaussian elements and $\|z\|_{\psi}\leq 1$.

• -

If $i\in\{n+1,n+2,\dots,n+m\}$ and $j\in\{1,2,\dots,n\}$, then $M_{ij}=\frac{1}{T}z^{T}R_{M}(3,4,1,2)z$, where $z\in\mathbb{R}^{(n+m)(T+T_{0}+1)}$ is a random vector with independent zero-mean sub-Gaussian elements and $\|z\|_{\psi}\leq 1$.

• -

If $i\in\{n+1,n+2,\dots,n+m\}$ and $j\in\{n+1,n+2,\dots,n+m\}$, then $M_{ij}=\frac{1}{T}z^{T}R^{(4)}_{M}(3,4,3,4)z$, where $z\in\mathbb{R}^{(n+m)(T+T_{0}+1)}$ is a random vector with independent zero-mean sub-Gaussian elements and $\|z\|_{\psi}\leq 1$.

Invoking the Hanson-Wright inequality and Lemma 7 for the aforementioned cases completes the proof.∎

B.3 The proof of Proposition 1

We need the following lemma:

Lemma 9**.**

We have

[TABLE]

Proof.

One can easily verify that

[TABLE]

and hence

[TABLE]

Therefore, with the assumption $\sigma_{w},\sigma_{w}\leq 1$ and the fact that $\sigma_{u},\sigma_{v}\leq\sqrt{5}\eta$ (the proof of which is simple and can be found, e.g., in [19]), one can write

[TABLE]

This completes the proof. ∎

Based on this lemma, we will take a similar approach to the proof of Theorem 3.1 in [7] to prove the correct sparsity recovery of the system matrices.

Proof of Proposition 1: To prove this proposition, we need to show that the conditions of Lemma 1 holds with high probability. To ensure that the first condition on $G$ implies the second one, it suffices to have

[TABLE]

Noting that $0<\gamma<1$, one can verify that the following bound on $\lambda$ is enough to guarantee that the above inequality holds:

[TABLE]

Furthermore, to ensure the last two conditions on $M$, it suffices to have

[TABLE]

Based on the above analysis, it suffices to have

[TABLE]

in order to ensure the exact recovery with probability of at least $1-\delta$. First, we derive conditions under which (67a) holds. Based on Lemma 2, one needs to ensure the following inequalities

[TABLE]

where (68c) is a technical condition that is required by Lemma 2. It can be easily verified that (68a) is satisfied with the choice of

[TABLE]

Based on the chosen value for $\lambda$ and in order to satisfy (68b), we should have the following lower bound on $T$

[TABLE]

Similarly, to ensure the validity of (68c), we should have

[TABLE]

Now, we will derive the conditions under which (67b) is satisfied using Lemma 3. To this goal, first we need to show that the following condition is satisfied:

[TABLE]

which is reduced to

[TABLE]

with the choice of $\epsilon=\frac{\gamma C_{\min}}{12\sqrt{k}}$. However, the above inequality implies that

[TABLE]

A sufficient condition for the correctness of the above inequality is to have $k\geq 2$. To see this, note that

[TABLE]

where the last inequality is due to Lemma 9. Therefore,

[TABLE]

which implies $k\geq 2$. Finally, to verify (67b) and according to Lemma 3, it suffices to have

[TABLE]

This implies that

[TABLE]

Based on the above analysis, the inequalities  (70), (71), and (78) impose lower bounds on $T$. Comparing these inequalities with (20), one can verify that the latter dominates all of them. This completes the proof.∎

B.4 Proof of Lemma 4

To prove this lemma, first we introduce the KKT conditions for (11).

Lemma 10** (KKT conditions).**

$\hat{\Psi}_{:,j}$* is an optimal solution for (11) if and only if it satisfies*

[TABLE]

for some $S\in\partial\|\hat{\Psi}_{:,j}\|_{1}$, where $\partial\|\hat{\Psi}_{:,j}\|_{1}$ is the sub-differential of $\|\cdot\|_{1}$ at $\hat{\Psi}_{:,j}$.

Proof.

The proof is trivial and is omitted for brevity. ∎

The following lemma is an immediate consequence of the KKT conditions.

Lemma 11**.**

Assuming that (11) recovers the correct sparsity pattern of $\Psi^{*}_{:,j}$, the following equalities hold for $E=\hat{\Psi}_{:,j}-\Psi^{*}_{:,j}$:

[TABLE]

Proof.

Due to the correct sparsity recovery, we have $E_{\mathcal{A}_{j}^{c}}=0$. This, together with the KKT conditions imply that

[TABLE]

Solving the above equation with respect to $E_{\mathcal{A}_{j}}$ will conclude the proof. ∎

Proof of Lemma 4: Based on Lemma 11, one can write

[TABLE]

In what follows, we will provide a bound for each term in the above inequality. For $Z_{2}$, one can write

[TABLE]

On the other hand, we have

[TABLE]

and therefore

[TABLE]

This leads to

[TABLE]

where the last inequality is due to the assumption (14). Combining the above inequality with (B.4) gives rise to

[TABLE]

Now we will bound $Z_{1}$. Similar to $Z_{2}$, we have

[TABLE]

Putting together (B.4) and (88) completes the proof.∎