Sample Complexity of Sparse System Identification Problem

TL;DR

This paper introduces a sparsity-promoting estimator for identifying high-dimensional sparse linear systems efficiently, requiring fewer samples than the system dimension and capable of exact sparsity recovery.

Contribution

It proposes a novel block-regularized estimator that achieves accurate sparse system identification with fewer samples and guarantees exact sparsity recovery in high-dimensional settings.

Findings

Estimator achieves small element-wise error with sufficient samples.

Sample complexity depends polynomially on block size and sparsity, logarithmically on system dimension.

Method guarantees exact recovery of the system's sparsity structure.

Abstract

In this paper, we study the system identification problem for sparse linear time-invariant systems. We propose a sparsity promoting block-regularized estimator to identify the dynamics of the system with only a limited number of input-state data samples. We characterize the properties of this estimator under high-dimensional scaling, where the growth rate of the system dimension is comparable to or even faster than that of the number of available sample trajectories. In particular, using contemporary results on high-dimensional statistics, we show that the proposed estimator results in a small element-wise error, provided that the number of sample trajectories is above a threshold. This threshold depends polynomially on the size of each block and the number of nonzero elements at different rows of input and state matrices, but only logarithmically on the system dimension. A by-product of this result is that the number of sample trajectories required for sparse system identification is significantly smaller than the dimension of the system. Furthermore, we show that, unlike the recently celebrated least-squares estimators for system identification problems, the method developed in this work is capable of \textit{exact recovery} of the underlying sparsity structure of the system with the aforementioned number of data samples. Extensive case studies on synthetically generated systems, physical mass-spring networks, and multi-agent systems are offered to demonstrate the effectiveness of the proposed method.

Equations

$$x[t+1]=Ax[t]+Bu[t]+w[t]$$

$$\underset{A,B}{\min }\sum _{i=1}^d\sum _{t=0}^{T-1}{\Vert x^{(i)}[t+1]-\left(A{x}^{(i)}[t]+B{u}^{(i)}[t]\right)\Vert }_2^2$$

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

$$Y=X\cdot \Theta +W$$

$${\Theta }_{\text{ls}}=\arg \underset{\Theta }{\min }\Vert Y-X\Theta {\Vert }_F^2$$

$${\Theta }_{\text{ls}}=(X^{\top }X{)}^{-1}{X}^{\top }Y={\Theta }^{\ast }+(X^{\top }X{)}^{-1}{X}^{\top }W$$

$$E={\Theta }_{\text{ls}}-{\Theta }^{\ast }=(X^{\top }X{)}^{-1}{X}^{\top }W$$

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

$$\hat{\Theta }=\arg \underset{\Theta }{\min }\frac{1}{2d}\Vert Y-X\Theta {\Vert }_F^2+{\lambda }_d\Vert \Theta {\Vert }_{\mathrm{block}}$$

$$x^{(i)}[T-1]=$$

$$+A^{T-2}{w}^{(i)}[0]+A^{T-3}{w}^{(i)}[1]+\cdots +w^{(i)}[T-2]$$

$$X_{i,:}^{\top }\sim N\left(0,\tilde{\Sigma }\right)$$

$$\tilde{\Sigma }=\left[\begin{array}{cc}{F}_T{\Sigma }_u{F}_T^{\top }+G_T{\Sigma }_w{G}_T^{\top }& 0\\ 0& {\Sigma }_u\end{array}\right]$$

$$F_T=[A^{T-2}B{A}^{T-3}B\cdots B]$$

$$G_T=[A^{T-2}{A}^{T-3}\cdots I]$$

$$n_{\max }=\underset{1\le i\le \bar{n}}{\max }{n}_i,$$

$$p_{\max }=\max \left\{n_{\max },m_{\max }\right\},$$

$${\sigma }_{\max }^2=\underset{1\le i\le n+m}{\max }{\tilde{\Sigma }}_{ii}$$

$$\underset{j=1,...,\bar{n}}{\max }\left\{\underset{i=1,...,|{\mathcal{A}}_j^c|}{\max }{\Vert {\tilde{\Sigma }}_{(i),{\mathcal{A}}_j}({\tilde{\Sigma }}_{{\mathcal{A}}_j,{\mathcal{A}}_j}{)}^{-1}\Vert }_1\right\}\le 1-\gamma$$

$${\Lambda }_{\min }\le {\lambda }_{\min }(\tilde{\Sigma })\le {\lambda }_{\max }(\tilde{\Sigma })\le {\Lambda }_{\max }$$

$$t_{\min }\le \underset{1\le j\le \bar{n}}{\min }\underset{i\in {\mathcal{A}}_j}{\min }{\Vert {{\Theta }^{\ast }}^{(i,j)}\Vert }_{\infty }$$

$$n_{\max }=O\left((\bar{n}+\bar{m}{)}^{{\alpha }_n}\right)$$

$$m_{\max }=O\left((\bar{n}+\bar{m}{)}^{{\alpha }_m}\right)$$

$${\lambda }_d=\Theta \left({\sigma }_{\max }\sqrt{\frac{\log (n+m)+\log (1/\delta )}{d}}\right)$$

$$d=\Omega \left(\kappa (\tilde{\Sigma }{)}^2{k}_{\max }\left(\log (n+m)+\log (1/\delta )\right)\right)$$

$$\displaystyle\|\hat{\Theta}\!-\!\Theta^{*}\|_{\infty}\!=\!O\Bigg{(}\kappa(\tilde{\Sigma})\left(1\!+\!\sqrt{\frac{k_{\max}(k_{\max}+\log({n}+{m})+\log(1/\delta))}{d}}\right)\sqrt{\frac{\log({n}+{m})+\log(1/\delta)}{d}}\Bigg{)}$$

$${\lambda }_d=\Theta \left({\sigma }_{\max }\sqrt{\frac{D\log (\bar{n}+\bar{m})+D^2\log (1/\delta )}{d}}\right)$$

$$d=\Omega \left(\kappa (\tilde{\Sigma }{)}^2{k}_{\max }\left(D\log (\bar{n}+\bar{m})+D^2\log (1/\delta )\right)\right)$$

$$\displaystyle g=\|\hat{\Theta}\!-\!\Theta^{*}\|_{\infty}\!=\!O\Bigg{(}\kappa(\tilde{\Sigma})\left(1\!+\!\sqrt{\frac{k_{\max}(k_{\max}n_{\max}\!+\!\log(\bar{n}+\bar{m})\!+\!\log(1/\delta))}{d}}\right)\sqrt{\frac{D\log(\bar{n}\!+\!\bar{m})\!+\!D^{2}\log(1/\delta)}{d}}\Bigg{)}$$

$${\lambda }_d=\Theta \left({\sigma }_{\max }{\log }^2(n+m)\sqrt{\frac{\log (1/\delta )}{d}}\right)$$

$$d=\Omega (\kappa (\tilde{\Sigma }{)}^2{k}_{\max }^2{\log }^4(n+m)\log (1/\delta ))$$

$$\Vert \hat{\Theta }-{\Theta }^{\ast }{\Vert }_{\infty }=O\left(\kappa (\tilde{\Sigma }){\log }^2(n+m)\sqrt{\frac{\log (1/\delta )}{d}}\right)$$

$$\sqrt{\frac{k_{\max }(k_{\max }{n}_{\max }+\log (\bar{n}+\bar{m})+\log (1/\delta ))}{d}}=O(1)$$

$${\lambda }_d=\Theta \left({\sigma }_{\max }(n+m{)}^{\left({\beta }_n+{\beta }_m\right)}\sqrt{\frac{\log (1/\delta )}{d}}\right)$$

$$d=\Omega (\kappa (\tilde{\Sigma }{)}^2{k}_{\max }^2(n+m{)}^{2\left({\beta }_n+{\beta }_m\right)}\log (1/\delta ))$$

Full text

I Introduction

With their ever-growing size and complexity, real-world dynamical systems are hard to model. Today’s systems are complex and large, often with a massive number of unknown parameters, which render them doomed to the so-called curse of dimensionality. Therefore, system operators should rely on simple and tractable estimation methods to identify the dynamics of the system via a limited number of recorded input-output interactions, and then design control policies to ensure the desired behavior of the entire system. The area of system identification is created to address this problem [1].

Despite its long history in control theory, most of the results on system identification deal with the asymptotic behavior of various estimation methods [1, 2, 3, 4]. Although these results shed light on the theoretical consistency of the existing methodologies, they are not applicable in the finite time/sample settings. In many applications, including neuroscience, transportation networks, and gene regulatory networks, the dimensionality of the system is overwhelmingly large, often surpassing the number of available input-output data [5, 6, 7]. Under such circumstances, the dynamics of the system should be estimated under the large dimension-small sample size regime and classical approaches for checking the asymptotic consistency of an estimator face major breakdowns. Simple examples of such failures are widespread in high-dimensional statistics. For instance, it is known that, given a number of independent and identically distributed (i.i.d.) samples with Gaussian distribution, none of the eigenvectors and eigenvalues of the sample covariance matrix are consistent estimators of their true counterparts if the sample size and the dimension of variables grow at the same rate [8]. As another example, it is well-known that the least-squares estimators, which are widely used in system identification problems, seize to exist uniquely when the sample size is smaller than the dimension of the system [9].

Recently, a special attention has been devoted to the sparse system identification problem, where the states and inputs are assumed to possess localized or low-order interactions. These methods include, but are not restricted to, selective $\ell_{1}$-regularized estimator [10], identification based on compressive sensing [11, 12], sparse estimation of polynomial system dynamics [13], kernel-based regularization [14], and low rank estimation in frequency domain [15]. On the other hand, with the unprecedented interest in data-driven control approaches, such as model-free reinforcement learning, robust control, and adaptive control [16, 17, 18], a question arises as to what the minimum number of input-output data samples should be to guarantee a small error in the estimated model. Answering this question has been the subject of many recent studies on the sample complexity of the system identification problem [19, 20, 21, 22, 23]. Most of these results are tailored to a specific type of dynamics, depend on the stability of the open-loop system, or do not exploit the a priori information about the structure of the system.

In this work, the objective is to employ modern results on high-dimensional statistics to reduce the sample complexity for one of the most fundamental classes of systems in control theory, namely linear time-invariant (LTI) systems with perfect state measurements. This type of dynamical systems forms the basis of many classical control problems, such as Linear Quadratic Regulator and Linear Quadratic Gaussian problems. Our results are built upon the fact that, in many practical large-scale systems, the states and inputs exhibit sparse interactions with one another, which in turn translates into a sparse representation of the state-space equations of the system. Driven by the existing non-asymptotic results on the classical Lasso problem, the main focus of this paper is on the block-regularized estimators for the system identification problem, where the goal is to promote sparsity on different blocks of input and state matrices. To this goal, the $\ell_{\infty}$-norms of the blocks are penalized instead of their $\ell_{1}$-norms. One motivation behind employing this type of estimator stems from topology extraction in consensus networks, especially in the multi-agent setting [24, 25]. In this problem, given a number of subsystems (agents) whose interactions are defined via an unknown sparse topology network, the objective is to estimate the state-space model governing the entire system based on a limited number of input-output sample trajectories. Since the subsystems have their own local state and input vectors with potentially different sizes, the parameters of the state-space model admit a block-sparse structure.

One objective is to show that the developed estimator for recovering the block-sparsity of the true dynamics is guaranteed to achieve infinitesimal estimation error with a small number of samples. In particular, we use an $\ell_{1}/\ell_{\infty}$-regularized least-squares estimator, i.e., a least-squares estimator accompanied by a $\ell_{\infty}$ regularizer on different blocks, and show that, with an appropriate scaling of the regularization coefficient, $\Omega(k_{\max}(D\log(\bar{n}+\bar{m})+D^{2}))$ sample trajectories are enough to guarantee a small estimation error with a high probability, where $k_{\max}$ is the maximum number of nonzero elements in the rows of input and state matrices, $D$ is the size of the largest block in these matrices, and $\bar{n}$ and $\bar{m}$ are the number of row blocks in the state and input matrices, respectively. This is a significant improvement over the recently derived sample complexity of $\Omega(n+m)$ for the least-squares estimator ($n$ and $m$ are the state and input dimensions, respectively), in the case where the system is sparse and the sizes of all blocks are small relative to the system dimension. While the traditional Lasso is heavily studied in the literature [26, 27], the high-dimensional behavior of the block-regularized estimator is less known when the dimensions of blocks are arbitrary. The paper [28] analyzes the high-dimensional consistency of this estimator when each block of the regression parameter is a row vector. Furthermore, that work assumes that the regression parameter consists of only one column of blocks. In an effort to make these results applicable to the block-sparse system identification problem, we significantly generalize the existing non-asymptotic properties to problems with an arbitrary number of blocks, each with general sizes.

Moreover, we derive upper bounds on the element-wise error of the proposed estimator. In particular, we prove that $\Omega(k_{\max}^{2}(D\log(\bar{n}+\bar{m})+D^{2}))$ sample trajectories is enough to ensure that the estimation error decreases at the rate $O(\sqrt{({D\log(\bar{n}+\bar{m})+D^{2}})/{d}})$, where $d$ is the number of available sample trajectories. We show that if the number of nonzero elements in the columns (in addition to the rows) of input and state matrices are upper bounded by $k_{\max}$, the operator norm of the estimation error of the proposed estimator is arbitrarily smaller than that of its un-regularized least-squares counterpart introduced in [22]. Another advantage of the proposed estimator over its least-squares analog is its exact recovery property. More specifically, we show that while the least-squares estimator is unable to identify the sparsity pattern of the input and state matrices for any finite number of samples, the proposed estimator recovers the true sparsity pattern of these matrices with a sublinear number of sample trajectories. It is worthwhile to mention that this work generalizes the results in [21], where the authors use a similar regularized estimator to learn the dynamics of a particular type of systems. However, [21] ignores the block structure of the system and assumes autonomy and inherent stability, all of which will be relaxed in this work. To demonstrate the efficacy of the developed regularized estimator, three case studies are offered on synthetically generated systems, mass-spring networks, and multi-agent systems.

This work is a significant extension of our previous conference papers on Lasso-type estimators for system identification [29] and non-asymptotic analysis of block-regularized linear regression problems [30]. In particular, by combining the properties of the block-regularized regression and the characteristics of LTI systems, we provide a unified sparsity-promoting framework for estimating the parameters of the system with arbitrary block structures. To this goal, we have generalized our theoretical results in [29] and [30] to account for partially-sparse structures. Furthermore, we have relaxed certain assumptions on the structure of the true system that were initially required in [29], and provided comprehensive discussions and more relevant simulations on the performance of the proposed method.

Notations: For a matrix $M$, the symbols $\|M\|_{F}$, $\|M\|_{2}$, $\|M\|_{1}$, and $\|M\|_{\infty}$ denote its Frobenius, operator, $\ell_{1}/\ell_{1}$, and $\ell_{\infty}/\ell_{\infty}$ norms, respectively. Furthermore, $\kappa(M)$ refers to its 2-norm condition number, i.e., the ratio between its maximum and minimum singular values. Given integer sets $I$ and $J$, the notation $M_{IJ}$ refers to the submatrix of $M$ whose rows and columns are indexed by $I$ and $J$, respectively. Given the sequences $f_{1}(n)$ and $f_{2}(n)$, the notations $f_{1}(n)=O(f_{2}(n))$ and $f_{1}(n)=\Omega(f_{2}(n))$ imply that there exist $c_{1}<\infty$ and $c_{2}>0$ such that $f_{1}(n)\leq c_{1}f_{2}(n)$ and $f_{1}(n)\geq c_{2}f_{2}(n)$, respectively. Finally, $f_{1}(n)=o(f_{2}(n))$ is used to show that $f_{1}(n)/f_{2}(n)\rightarrow 0$ as $n\rightarrow\infty$. A zero-mean Gaussian distribution with covariance $\Sigma$ is shown as $N(0,\Sigma)$. Given a function $f(x)$, the expression $\arg\min f(x)$ refers to its minimizer. For a set $\mathcal{I}$, the symbol $|\mathcal{I}|$ denotes its cardinality.

II Problem Formulation

Consider the LTI system

[TABLE]

where $t$ is the time step, $A\in\mathbb{R}^{n\times n}$ is the state matrix, and $B\in\mathbb{R}^{n\times m}$ is the input matrix. Furthermore, $x[t]\in\mathbb{R}^{n}$, $u[t]\in\mathbb{R}^{m}$, and $w[t]\in\mathbb{R}^{n}$ are the state, input, and disturbance vectors at time $t$, respectively. The dimension of the system is defined as $m+n$. It is assumed that the input disturbance vectors are identically distributed and independent with distribution $N(0,\Sigma_{w})$ across different times. In this work, we assume that the matrices $A$ and $B$ are sparse and the goal is to estimate them 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. The $i^{\text{th}}$ sample trajectory $\{(x^{(i)}[\tau],u^{(i)}[\tau])\}_{\tau=0}^{T}$ is obtained by running the system from $t=0$ to $t=T$ and collecting the input and state vectors. The sparsity assumption on $A$ and $B$ is practical in many applications because of two reasons:

• •

In order to model many large-scale real-world systems accurately, one needs to consider an overwhelmingly large number of internal states. However, it is often the case that the interactions between different states and inputs obey a sparse structure, which translates into a sparse pattern in state and input matrices. An important example of these types of problems is multi-agent systems, where the agents (subsystems) interact with one another via a sparse communication network.

• •

More generally, it is well-known that there may not be a unique representation of the state-space model of the system. For example, one can entirely change the state and input matrices in (1) via linear/nonlinear transformations to arrive at a different, but equally accurate state-space equation. Furthermore, the dynamics of the system may have a sufficiently accurate sparse approximation. The question of interest is whether it is possible to design a data-driven method in order to estimate the sparsest state-space representation of the system. Answering this question is crucial, specially in the context of a distributed control problem where the goal is to design a decentralized controller whose structure respects the dynamics of the system [31, 32, 33].

Given the sample trajectories $\{(x^{(i)}[\tau],u^{(i)}[\tau])\}_{\tau=0}^{\top}$ for $i=1,2,...,d$, one can obtain an estimate of $(A,B)$ by solving the following least-squares optimization problem:

[TABLE]

In order to describe the behavior of the least-squares estimator, define

[TABLE]

for every sample trajectory $i=1,2,...,d$. Furthermore, let $Y$, $X$, and $W$ be defined as vertical concatenations of $Y^{(i)}$, $X^{(i)}$, and $W^{(i)}$ for $i=1,2,...,d$, respectively. Finally, denote $\Theta=[A\ \ B]^{\top}$ as the unknown system parameter and $\Theta^{*}$ as its true value. Based on these definitions, it follows from (1) that

[TABLE]

The system identification problem is then reduced to estimating $\Theta$ based on the observation matrix $Y$ and the design matrix $X$. Consider the following least-squares estimator:

[TABLE]

One can easily verify the equivalence of (2) and (5). The optimal solution of (5) can be written as

[TABLE]

Notice that ${\Theta}_{\text{ls}}$ is well-defined and unique if and only if $X^{\top}X$ is invertible, which necessitates $d\geq n+m$. The estimation error is then defined as

[TABLE]

Thus, one needs to study the behavior of $(X^{\top}X)^{-1}X^{\top}W$ in order to control the estimation error of the least-squares estimator. However, since the state of the system at time $t$ is affected by random input disturbances at times $0,1,...t-1$, the matrices $X$ and $W$ are correlated, which renders (7) hard to analyze. In order to circumvent this issue, [34] simplifies the estimator and considers only the state of the system at time $T$ in $Y^{(i)}$. By ignoring the first $T-1$ rows in $Y^{(i)}$, $X^{(i)}$, and $W^{(i)}$, one can ensure that the random matrix $(X^{\top}X)^{-1}X^{\top}$ is independent of $W$. Therefore, it is assumed in the sequel that

[TABLE]

With this simplification, [34] shows that, with input vectors $u^{(i)}[t]$ chosen randomly from $N(0,\Sigma_{u})$ for every $t=1,2,...,T-1$ and $i=1,2,...,d$, the least-squares estimator requires at least $d=\Omega(m+n+\log(1/\delta))$ sample trajectories to guarantee $\|E\|_{2}=\mathcal{O}\left(\sqrt{{(m+n)\log(1/\delta)}/{d}}\right)$ with probability of at least $1-\delta$. In what follows, a block-regularized estimator will be introduced that exploits the underlying sparsity structure of the system dynamics to significantly reduce the number of sample trajectories for an accurate estimation of the parameters.

III Main Results

Suppose that $A$ and $B$ can be partitioned as $A=[A^{(i,j)}]$ and $B=[B^{(k,l)}]$ where $(i,j)\in\{1,...,\bar{n}\}\times\{1,...,\bar{n}\}$ and $(k,l)\in\{1,...,\bar{n}\}\times\{1,...,\bar{m}\}$. $A^{(i,j)}$ is the $(i,j)^{\text{th}}$ block of $A$ with size $n_{i}\times n_{j}$. Similarly, $B^{(k,l)}$ is the $(k,l)^{\text{th}}$ block of $B$ with size $n_{k}\times m_{l}$. Note that $\sum_{i=1}^{\bar{n}}n_{i}=n$ and $\sum_{i=1}^{\bar{m}}m_{i}=m$. Suppose that it is known a priori that all elements in each block $A^{(i,j)}$ or $B^{(k,l)}$ are simultaneously zero or nonzero. This implies that, as long as one element in $A^{(i,j)}$ or $B^{(k,l)}$ is nonzero, there is no reason to promote sparsity in the remaining elements of the corresponding block. Clearly, this kind of block-sparsity constraint is not correctly reflected in (2). To streamline the presentation, we use the notation $\Theta=[A\ \ B]^{\top}$. Note that $\Theta^{(i,j)}=(A^{(j,i)})^{\top}$ for $i\in\{1,...,\bar{n}\}$ and $\Theta^{(i,j)}=(B^{(j,i-\bar{n})})^{\top}$ for $i\in\{\bar{n}+1,...,\bar{n}+\bar{m}\}$. In order to recover the true block-sparsity of $A$ and $B$, one can resort to an $\ell_{1}/\ell_{\infty}$ variant of the Lasso problem—known as the block-regularized least-squares (or simply block-regularized) problem:

[TABLE]

where $\|\Theta\|_{\mathrm{block}}$ is defined as the summation of $\|\Theta^{(i,j)}\|_{\infty}$ over $(i,j)\in\{1,...,\bar{n}+\bar{m}\}\times\{1,...,\bar{n}\}$. Suppose that $D$ denotes the maximum size of the blocks in $\Theta$, i.e., the multiplication of its number of rows and columns. Under the sparsity assumption on $(A,B)$, we will show that the non-asymptotic statistical properties of $\hat{\Theta}$ significantly outperform those of ${\Theta}_{\text{ls}}$. In particular, the primary objective is to prove that $\|\hat{\Theta}-\Theta^{*}\|_{\infty}$ decreases at the rate $\mathcal{O}(\sqrt{{D\log(n+m)+D^{2}\log(1/\delta)}/{{d}}})$ with probability of at least $1-\delta$ and with an appropriate scaling of the regularization coefficient, provided that $d=\Omega\left(k_{\max}^{2}\left(D\log(\bar{n}+\bar{m})+D^{2}\log(1/\delta)\right)\right)$. Here, $k_{\max}$ is the maximum number of nonzero elements in the columns of $[A\ \ B]^{\top}$. Comparing this number with the required lower bound $\Omega(n+m+\log(1/\delta))$ on the number of sample trajectories for the least-squares estimator, we conclude that the proposed method needs significantly less number of samples when $A$ and $B$ are sparse. The third objective is to prove that this method is able to find the correct block-sparsity structure of $A$ and $B$ with high probability. In contrast, it will be shown that the solution of the least-squares estimator is fully dense for any finite number of sample trajectories, and hence, it cannot correctly extract the sparsity structures of $A$ and $B$. We will showcase the superior performance of the block-regularized estimator both in sparsity identification and estimation accuracy in simulations.

To present the main results of this work, first note that

[TABLE]

where, without loss of generality, the initial state is assumed to be zero for every sample trajectory. The results can be readily extended to the case where the initial state is an unknown random vector with Gaussian distribution. Suppose that $u^{(i)}[t]$ and $w^{(i)}[t]$ are i.i.d samples of $N(0,\Sigma_{u})$ and $N(0,\Sigma_{w})$, respectively. Therefore, (III) and (8) imply that

[TABLE]

where $X_{i,:}$ is the $i^{\text{th}}$ row of $X$ and

[TABLE]

Define $\mathcal{A}_{j}(\Theta)=\{i:\Theta^{(i,j)}\not=0\}$. Unless stated otherwise, $\mathcal{A}_{j}$ is used to refer to $\mathcal{A}_{j}(\Theta^{*})$. Define $\mathcal{A}^{c}_{j}$ as the complement of $\mathcal{A}_{j}$. For $\mathcal{T}\subseteq\{1,...,\bar{n}+\bar{m}\}$, denote ${I}(\mathcal{T})$ as the index set of rows in $\Theta^{*}$ corresponding to the blocks $\{{\Theta^{*}}^{(i,:)}:i\in\mathcal{T}\}$. For an index set $\mathcal{U}$, define $X_{\mathcal{U}}$ as a $d\times|\mathcal{U}|$ submatrix of $X$ after removing the columns with indices not belonging to $\mathcal{U}$. With a slight abuse of notation, $X_{(i)}$, $X_{\mathcal{A}_{j}}$, and $X_{\mathcal{A}_{j}^{c}}$ are used to denote $X_{{I}(\{i\})}$, $X_{{I(\mathcal{A}_{j})}}$, and $X_{{I}(\mathcal{A}_{j}^{c})}$ when there is no ambiguity. Similarly, $\Sigma_{(i),\mathcal{A}_{j}}$ and $\Sigma_{\mathcal{A}_{j},\mathcal{A}_{j}}$ are used in lieu of $\Sigma_{{I}(\{i\}),I(\mathcal{A}_{j})}$ and $\Sigma_{I(\mathcal{A}_{j}),I(\mathcal{A}_{j})}$, respectively. Denote $k_{j}$ as the maximum number of nonzero elements in different columns of ${\Theta^{*}}^{(:,j)}$ which is the $j^{\text{th}}$ block column of $\Theta^{*}$. Finally, define

[TABLE]

The following set of assumptions plays a key role in deriving the main result of this paper:

Assumption 1**.**

By fixing the time horizon $T$, the following conditions hold for all finite system dimensions:

• A1.

(Mutual Incoherence Property): There exists a number $\gamma\in(0,1]$ such that

[TABLE]

• A2.

(Bounded eigenvalue): There exist numbers $0<\Lambda_{\min}<\infty$ and $0<\Lambda_{\max}<\infty$ such that

[TABLE]

• A3.

(Bounded minimum value): There exists a number $t_{\min}>0$ such that

[TABLE]

• A4.

(Block sizes): There exist numbers $\alpha_{n},\alpha_{m}<\infty$ such that

[TABLE]

The mutual incoherence property in Assumption A1 is a commonly known condition for the exact recovery of unknown parameters in compressive sensing and classical Lasso problems [27, 35, 36, 37]. This assumption entails that the effect of those submatrices of $\tilde{\Sigma}$ corresponding to zero (unimportant) elements of $\Theta$ on the remaining entries of $\tilde{\Sigma}$ should not be large. Roughly speaking, this condition guarantees that the unknown parameters are recoverable in the noiseless scenario, i.e. when $W=0$. If the recovery cannot be guaranteed in the noise-free setting, then there is little hope for the block-regularized estimator to recover the true structure of $A$ and $B$ when the system is subject to noise. This assumption is satisfied in all of our simulations.

The bounded eigenvalue condition in Assumption A2 entails that the condition number of $\tilde{\Sigma}$ is bounded away from [math] and $\infty$ for all finite system dimensions. However, as will be shown later, the values of $\Lambda_{\min}$ and $\Lambda_{\max}$ can change with respect to the time horizon $T$. In particular, it will be later shown that for highly unstable systems, $\tilde{\Sigma}$ becomes severely ill-conditioned as the time horizon increases, which in turn makes the system identification problem difficult to solve. Furthermore, this assumption implies that there exists a constant $\bar{\sigma}_{\max}^{2}<\infty$ such that $\max_{1\leq i\leq n+m}\Sigma_{ii}\leq\bar{\sigma}_{\max}^{2}$ for every finite system dimension.

Assumption A3 implies that, independent of the system dimensions, there always exists a strictly positive gap between the zero and nonzero elements of $A$ and $B$. This assumption holds in almost all practical settings and will facilitate the exact sparsity recovery of the parameters of the system. It is worthwhile to mention that this condition can be relaxed to the case where $t_{\min}$ is a decreasing function of the system dimension, provided that its decrease rate is slow enough. To streamline the presentation, we did not consider this case in the paper.

Finally, Assumption A4 requires that the maximum size of the blocks in $\Theta^{*}$ be polynomially bounded by the number of its block columns. For instance, $\bar{n}=O(1)$ and $\bar{m}=O(1)$ violate this assumption since it implies that $n_{\max}=\Omega((\bar{n}+\bar{m})^{\log n})$ and $m_{\max}=\Omega((\bar{n}+\bar{m})^{\log m})$. It is worthwhile to mention that Assumption A4 results in $k_{\max}=O((\bar{n}+\bar{m})^{\alpha_{k}})$ for some number $\alpha_{k}<\infty$; this will be used later in the paper.

Remark 1**.**

Note that, due to Assumption A2, $\kappa(\tilde{\Sigma})=O(1)$. However, this quantity will not be removed from the big-$O$ analysis of our subsequent theorems and corollaries to demonstrate its effect on the high-dimensional properties of the developed estimator.

The first theorem characterizes the non-asymptotic behavior of the Lasso-type estimator for the system identification problem where each block has the size $1\times 1$.

Theorem 1** (element-wise regularization).**

Suppose that $n_{\max}=m_{\max}=1$. Upon choosing

[TABLE]

the following statements hold with probability of at least $1-\delta$:

$\hat{\Theta}$* is unique and has the same nonzero elements as $\Theta^{*}$.*

• 2.

We have

[TABLE]

Theorem 1 states that, for sparse system dynamics, the support of $\hat{\Theta}$ will be recovered exactly, provided that the number of samples exceeds a threshold that is only logarithmic in the system dimension. This is specifically useful for large-scale systems where $k_{\max}$ is significantly smaller than $n+m$. Intuitively, this corresponds to those systems whose states are affected only by a few, and possibly local, states and inputs.

Define $D=p_{\max}n_{\max}$, which is the maximum size of the blocks in $\Theta$. The next theorem extends the result of Theorem 1 to arbitrary block sizes.

Theorem 2** (block-wise regularization).**

Upon choosing

[TABLE]

the following statements hold with probability of at least $1-\delta$:

$\hat{\Theta}$* is unique and has the same nonzero blocks as $\Theta^{*}$.*

• 2.

We have

[TABLE]

Theorem 2 shows that the minimum number of required sample trajectories is a quadratic function of the maximum block size. Therefore, only a small number of samples is enough to guarantee the uniqueness, exact block-sparsity recovery, and small estimation error for sparse systems, assuming that the sizes of the blocks are significantly smaller than the system dimensions.

Another important observation can be made from Theorems 1 and 2. Notice that the minimum number of required sample trajectories and the element-wise error of the estimated parameter depends on $\kappa(\tilde{\Sigma})$. Upon assuming that $\Sigma_{u}$ and $\Sigma_{w}$ are identity matrices, $\kappa(\tilde{\Sigma})$ reduces to $\kappa(F_{T}F_{T}^{\top}+G_{T}G_{T}^{\top})$. $F_{T}$ and $G_{T}$ are commonly known as finite-time controllability matrices for the input and disturbance noise, respectively. Roughly speaking, $\kappa(F_{T}F_{T}^{\top}+G_{T}G_{T}^{\top})$ quantifies the ratio between the eigenvalues of the easiest- and hardest-to-identify modes of the system. Therefore, Theorems 1 and 2 imply that only a small number of samples is required to accurately identify the dynamics of the system if all of its modes are easily excitable. The dependency of the estimation error on the modes of the system is also reflected in the non-asymptotic error bound of the least-squares estimator in [34]. This is completely in line with the conventional results on the identifiability of dynamical systems: independent of the method in use, it is significantly harder to identify the parameters of the system accurately if it possesses nearly-hidden modes. The connection between the identifiability of the system and the number of required sample trajectories to guarantee a small estimation error will be elaborated through different case studies in Section V.

In what follows, the effect of different scalings of block sizes on the minimum number of required sample trajectories and the element-wise error of the proposed estimator will be elaborated. In particular, it will be shown that, as the size of the blocks grows, a proportional increase in the number of sample trajectories is required for the successful recovery of the system parameters.

Corollary 1**.**

Assume that $n_{\max}=O(1)$ and $m_{\max}=O(1)$. Then, the block-regularized estimator behaves similarly to the element-wise regularized estimator, i.e., the respective choices of (18a) and (18b) for $\lambda_{d}$ and $d$ are enough to guarantee exact block-sparsity recovery of $\Theta^{*}$ and the element-wise error of (19) with probability of at least $1-\delta$.

Proof.

The proof is immediate by noting that $D=O(1)$ and $O(\log(\bar{n}+\bar{m}))=O(\log({n}+{m}))$. ∎

Corollary 1 shows that the block-regularized estimator behaves similarly to its element-wise counterpart if the sizes of the blocks do not change as the system dimension grows.

Corollary 2**.**

Assume that $n_{\max}=O(\log n)$ and $m_{\max}=O(\log m)$. Then,

[TABLE]

is enough to guarantee the exact block-sparsity recovery of $\Theta^{*}$ and

[TABLE]

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

Proof.

First, note that $D=O(\log^{2}(n+m))$. Therefore, (20a) is equivalent to (22a). Furthermore, (22b) implies that

[TABLE]

This reduces (21) to (23), thereby completing the proof. ∎

Corollary 2 studies the behavior of the proposed estimator when the sizes of the blocks increase logarithmically in the system dimensions. The effect of such growth in the size of the blocks is reflected in the minimum number of necessary sample trajectories. However, it can be seen that under the scenarios where $k_{\max}$ is small, i.e. when the system is sparse, this number can still be significantly smaller than the system dimensions.

Corollary 3**.**

Assume that $n_{\max}=O(n^{\beta_{n}})$ and $m_{\max}=O(m^{\beta_{m}})$ for some $\beta_{n}>0$ and $\beta_{m}>0$. Then,

[TABLE]

is enough to guarantee the exact sparsity recovery of $\Theta^{*}$ and

[TABLE]

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

Proof.

The proof is similar to that of Corollary 2. ∎

Corollary 3 is the analog of Corollaries 1 and 2 for the polynomial scaling of the block size. Similar to the previous cases, it can be seen that the size of the required sample trajectories heavily depends on the growth rate of the maximum block size of $\Theta$. Although the sampling rate is still sublinear when $\beta_{n}+\beta_{m}<1/2$, it may surpass the system dimension if $\beta_{n}+\beta_{m}>1/2$. A question arises as to whether one can resort to the ordinary least-squares estimator in lieu of the proposed block-regularized estimator for the cases where $\beta_{n}+\beta_{m}>1/2$ since the proposed estimator requires $d=\Omega((n+m)^{1+\epsilon}\log(1/\delta))$ for some $\epsilon>0$ whereas $d=\Theta(n+m+\log(1/\delta))$ in enough to guarantee the uniqueness of the least-squares estimator. However, in what follows, we will prove that the least-squares estimator does not extract the correct sparsity structure of $\Theta$ for any finite number of sample trajectories.

Theorem 3**.**

If $A$ and $B$ are not fully dense matrices, ${\Theta}_{ls}$ does not recover the support of $\Theta^{*}$ for any finite number of sample trajectories with probability 1.

Proof.

Define $R=((X^{\top}X)^{-1}X^{\top})^{\top}$, and note that $R$ and $W$ are independent random variables due to the construction of $X$. Now, suppose that $\Theta_{ij}^{*}=0$. We show that, with probability zero, $E_{ij}=|(\Theta_{\text{ls}})_{ij}-\Theta_{ij}^{*}|=0$ holds. Note that $E_{ij}=R_{:,i}^{\top}W_{:,j}$. If $R_{:,i}\not=0$, then $E_{ij}$ is a linear combination (with at least one nonzero coefficient) of identically distributed normal random variables with mean zero and variance $(\Sigma_{w})_{jj}$. Since $R_{:,i}$ and $W_{:,j}$ are independent, we have $E_{ij}=0$ with probability zero. Now, assume that $R_{:,i}=0$. This means that the $i^{\text{th}}$ row of $R^{\top}$ is a zero vector. This, in turn, implies that the $i^{\text{th}}$ row of $R^{\top}X$ is zero. However, $R^{\top}X=(X^{\top}X)^{-1}X^{\top}X=I$, which is a contradiction. This completes the proof. ∎

Define $h(n,m)=\sqrt{{(n+m)\log(1/\delta)}/{d}}$ and recall that $\|\Theta_{\mathrm{ls}}-\Theta^{*}\|_{2}=O(h(n,m))$. In the next corollary, we show that, under additional sparsity conditions, the operator norm of the estimation error for $\hat{\Theta}$ becomes arbitrarily smaller than $h(n,m)$ as the system dimension grows.

Corollary 4**.**

Assume that the number of nonzero elements at different rows and columns of $\Theta^{*}$ is upper bounded by $k_{\max}$. Furthermore, suppose that $\lambda_{d}$ satisfies (20a) and

[TABLE]

Then, we have

[TABLE]

with probability of at least $1-\delta$. Furthermore, we have

[TABLE]

provided that

[TABLE]

Proof.

One can use Holder’s inequality to write

[TABLE]

On the other hand, (24) can be verified under (27). Combined with (31) and Theorem 2, this certifies the validity of (28). It remains to prove the correctness of (29). Note that under (30), we have

[TABLE]

Combined with the definitions of $h(n,m)$ and $v(n,m)$, this completes the proof. ∎

Corollary 4 describes the settings under which our proposed method significantly outperforms the least-squares estimator in terms of the operator norm of the errors. This improvement is more evident for those systems where the states and inputs have sparse interactions and the block sizes in $A$ and $B$ are smaller than the system dimensions. A class of such systems is multi-agent networks where the agents interact only locally and their total number dominates the dimension of each individual agent.

IV Proofs

In this section, the proof of Theorem 2 will be presented. Furthermore, it will be shown that Theorem 1 is a special case of this theorem. A number of preliminary definitions and lemmas are required to present the proof of Theorem 2.

Definition 1** (sub-Gaussian random variable).**

A zero-mean random variable $x$ is sub-Gaussian with parameter $\sigma^{2}$ if there exists a constant number $c<\infty$ such that

[TABLE]

Lemma 1**.**

Given a set of zero-mean sub-Gaussian variables $x_{i}$ with parameters $\sigma_{i}$ for $i=1,2,...,m$, the inequality

[TABLE]

holds for some constant $c<\infty$.

Define $I_{d}$ as the $d\times d$ identity matrix. The next two lemmas are borrowed from [28] and [26], respectively.

Lemma 2**.**

Given a set of random vectors $X_{i}\sim N(0,\sigma_{i}^{2}I_{d})$ for $i=1,2,...,m$ and $d>2\log m$, the inequality

[TABLE]

holds, where $\sigma=\max_{i}\sigma_{i}$.

Lemma 3**.**

Consider a matrix $X\in\mathbb{R}^{m\times n}$ whose rows are drawn from $N(0,\Sigma)$. Assuming that $n\leq m$, we have

[TABLE]

for every $n\leq t\leq m$.

The basic inequalities given below will be used frequently in our subsequent arguments.

Lemma 4**.**

The following statements hold true:

• •

Given a number of (not necessarily independent) events $\mathcal{T}_{i}$ for $i=1,2,...,n$, the following inequality is satisfied:

[TABLE]

• •

Given events $\mathcal{B}$ and $\mathcal{C}$ together with the complement of $\mathcal{C}$, denoted as $\mathcal{C}^{c}$, the following inequality holds:

[TABLE]

The next lemma characterizes the first-order optimality conditions for (9).

Lemma 5** (KKT conditions).**

$\hat{\Theta}$* is an optimal solution for (9) if and only if it satisfies*

[TABLE]

for some $\hat{S}\in\mathbb{R}^{(n+m)\times n}\in\partial\|\hat{\Theta}\|_{\mathrm{block}}$, where $\partial\|\hat{\Theta}\|_{\mathrm{block}}$ denotes the sub-differential of $\|\cdot\|_{\mathrm{block}}$ at $\hat{\Theta}$.

Proof.

The proof is straightforward and omitted for brevity. ∎

$\hat{S}_{\mathcal{A}}$ and $\hat{S}_{\mathcal{A}^{c}}$ are obtained by removing those blocks of $\hat{S}$ with indices not belonging to $\mathcal{A}$ and $\mathcal{A}^{c}$, respectively. The equation (4) can be reformulated as the set of linear equations

[TABLE]

where $Y^{(:,j)}$, $\Theta^{(:,j)}$, and $W^{(:,j)}$ are the $j^{\text{th}}$ block column of $Y$, $\Theta$, and $W$, respectively. Based on this definition, consider the following set of block-regularized subproblems:

[TABLE]

Define $D_{j}=p_{\max}n_{j}$. The next two lemmas are at the core of our proof for Theorem 2.

Lemma 6** (No false positives).**

Given arbitrary constants $c_{1},c_{2}>1$, suppose that $\lambda_{d}$ and $d$ are chosen such that

[TABLE]

Then, with probability of at least

[TABLE]

$\hat{\Theta}^{(:,j)}$* is unique and its nonzero blocks exclude the zero blocks of ${\Theta^{*}}^{(:,j)}$. In other words, $\hat{\Theta}^{(:,j)}$ does not have any false positives.*

Recall that due to Assumption A4, one can write $n_{\max}=O\left((\bar{n}+\bar{m})^{\alpha_{n}}\right)$ and $k_{\max}=O\left((\bar{n}+\bar{m})^{\alpha_{k}}\right)$ for some $\alpha_{n}\geq 0$ and $\alpha_{k}\geq 0$.

Lemma 7** (Element-wise error).**

Given arbitrary constants $c_{3}>0$ and $c_{4}>1$, suppose that $\hat{\Theta}$ is unique and the set of its nonzero blocks excludes the zero blocks of $\Theta^{*}$. Then, with probability of at least

[TABLE]

we have

[TABLE]

Furthermore, the zero blocks of $\hat{\Theta}^{(:,j)}$ exclude the nonzero blocks of ${\Theta^{*}}^{(:,j)}$ if $\min_{i\in\mathcal{A}_{j}}\|\Theta^{(i,j)}\|_{\infty}>g_{j}$. In other words, $\hat{\Theta}^{(:,j)}$ does not have any false negatives if $\min_{i\in\mathcal{A}_{j}}\|\Theta^{(i,j)}\|_{\infty}>g_{j}$.

In what follows, we will present some preliminaries that are essential in proving Lemmas 6 and 7. Notice that $\hat{S}$ and $W$ have the same dimensions as $\hat{\Theta}$, and hence, can be similarly partitioned into different blocks. Since Lemmas 6 and 7 hold for any given column block index $j$, $\Theta^{(i,j)}$ and $\mathcal{A}_{j}$ will be referred to as $\Theta^{(i)}$ and $\mathcal{A}$ in order to streamline the presentation.

Lemma 8**.**

$Q\in\partial\|\tilde{\Theta}\|_{\mathrm{block}}$* if and only if the following conditions are satisfied for every $i\in\{1,2,...,\bar{n}+\bar{m}\}$:*

• •

If $\|\tilde{\Theta}^{(i)}\|_{\infty}\not=0$, define $M^{(i)}=\{(k,l):\tilde{\Theta}^{(i)}_{kl}=\|\tilde{\Theta}^{(i)}\|_{\infty}\}$. Then, $Q^{(i)}_{kl}=\eta_{kl}\cdot\mathrm{sign}(\tilde{\Theta}^{(i)}_{kl})$, where $\sum_{(k,l)\in M^{(i)}}\eta_{kl}=1$ and $\eta_{kl}=0$ if $(k,l)\not\in M^{(i)}$.

• •

If $\|\tilde{\Theta}^{(i)}\|_{\infty}=0$, then $\|Q^{(i)}\|_{1}\leq 1$.

The proofs of Lemmas 6 and 7 are based on the well-known primal-dual witness approach introduced in [26, 28], which is defined as follows:

Primal-dual witness approach ([26, 28]):

• •

Step 1: Define the restricted regularized problem as

[TABLE]

whose solution is unique if $X_{\mathcal{A}}^{\top}X_{\mathcal{A}}$ is invertible.

• •

Step 2: With a slight abuse of notation, $\tilde{\Theta}$ can be written as $(\tilde{\Theta}_{\mathcal{A}},0)$. Choose $\tilde{S}_{\mathcal{A}}$ as an element of the sub-differential $\partial\|\tilde{\Theta}_{\mathcal{A}}\|_{\mathrm{block}}$.

• •

Step 3: Find $\tilde{S}_{\mathcal{A}}^{c}$ by solving the KKT equations (39), given $\tilde{\Theta}$ and $\tilde{S}_{\mathcal{A}}$. Then, verify

[TABLE]

If (47) can be verified in the last step, it is said that the primal-dual witness (PDW) approach succeeds. The next lemma unveils a close relationship between the block-regularized estimator, the PDW approach, and the true regression parameter $\Theta^{*}$.

Lemma 9**.**

The following statements hold:

• •

If the PDW approach succeeds, then $\tilde{\Theta}$ is the unique optimal solution of (9), i.e. $\hat{\Theta}=\tilde{\Theta}$.

• •

Conversely, suppose that $\hat{\Theta}$ is the optimal solution of (9) such that $\hat{\Theta}^{(i)}=0$ for every $i\in\mathcal{A}^{c}$. Then, the PDW approach succeeds.

Proof.

The proof is a simple generalization of Lemma 2 in [28]. The details are omitted for brevity. ∎

Lemma 9 is the building block of our proofs for Lemmas 6 and 7. In particular, Lemma 9 indicates that in order to show that the solution of (40) is unique and excludes false positive errors, it is enough to verify that the PDW approach succeeds with high probability. Then, conditioned on the success of the PDW approach, our focus can be devoted to the optimal solution of the restricted problem (46) and bounding its difference from the true parameters.

Lemma 10**.**

Define $\tilde{\Theta}-\Theta=E$. The following equalities hold:

[TABLE]

Proof.

To verify (48b) and (48c), note that the KKT condition in Lemma 5 reduces to

[TABLE]

Solving (49a) with respect to $E_{\mathcal{A}}$ and substituting the solution in (49b) completes the proof. ∎

IV-A Proof of Lemma 6:

As shown in Lemma 9, it is enough to prove that the PDW succeeds with high probability. To this goal, we show that $\max_{i\in\mathcal{A}^{c}}\|\tilde{S}^{(i)}\|_{1}<1$ with high probability, which results in the success of the PDW approach. Lemma 10 yields that

[TABLE]

Similar to [28], we will show that $\max_{i\in\mathcal{A}^{c}}Z_{1}^{(i)}<\gamma/2$ and $\max_{i\in\mathcal{A}^{c}}Z_{2}^{(i)}<1-\gamma/2$ with high probability. First, consider $\max_{i\in\mathcal{A}^{c}}Z_{1}^{(i)}$. We have

[TABLE]

Given $X$, note that $R_{kl}^{(i)}$ is Gaussian with variance

[TABLE]

Moreover, $X_{\mathcal{A}}(X_{\mathcal{A}}^{\top}X_{\mathcal{A}})^{-1}X_{\mathcal{A}}^{\top}$ is an orthogonal projection onto the range of $X_{\mathcal{A}}$. Therefore,

[TABLE]

Define $p_{i}=n_{i}$ if $1\leq i\leq\bar{n}$ and $p_{i}=m_{i}$ if $\bar{n}+1\leq i\leq\bar{n}+\bar{m}$. Due to Lemma 2, the last inequality is upper bounded by ${4\sigma_{w}^{2}\sigma_{\max}^{2}}/{d\lambda_{d}^{2}}$ for every $k\in\{1,...,p_{i}\}$ with probability of at least $1-\exp(-d/2+\log p_{i})$ for $d>2\log p_{i}$. Conditioned on this event, one can write:

[TABLE]

which means that $\sum_{(k,l)\in\Theta^{(i)}}\epsilon_{kl}R_{kl}^{(i)}$ is sub-Gaussian with the parameter ${4D_{j}\sigma_{w}^{2}\sigma_{\max}^{2}}/{d\lambda_{d}^{2}}$. This implies that

[TABLE]

where we have used Lemma 1, the second statement of Lemma 4 and the facts that $p_{i}\leq p_{\max}$ and $|\mathcal{A}^{c}|\leq\bar{n}+\bar{m}$ in the last inequality. Now, setting $\zeta=\gamma/2$ and

[TABLE]

for some arbitrary constant $c_{1}>1$ yields that

[TABLE]

where the last inequality is due to the lower bound (42b) on $d$. Next, an upper bound on $\max_{i\in\mathcal{A}^{c}}Z_{2}^{(i)}$ will be derived. Since each row of $X$ is drawn from $N(0,\tilde{\Sigma})$, one can write the distribution of $X_{\mathcal{A}^{c}}^{\top}$, conditioned on $X_{\mathcal{A}}$ as

[TABLE]

Based on (58), one can verify that $\frac{1}{d}X_{\mathcal{A}^{c}}^{\top}X_{\mathcal{A}}(\frac{1}{d}X_{\mathcal{A}}^{\top}X_{\mathcal{A}})^{-1}\tilde{S}_{\mathcal{A}}$ has the same distribution as

[TABLE]

where $V$ is a random matrix with zero mean, covariance $\tilde{\Sigma}_{\mathcal{A}^{c}|\mathcal{A}}$, and independent of $X$. In light of the definition of $\tilde{\Sigma}_{\mathcal{A}^{c}|\mathcal{A}}$, it can be easily seen that the elements of $V$ are sub-Gaussian with parameters of at most $\sigma_{\max}^{2}$. This implies that

[TABLE]

where we have used the mutual incoherence property and the fact that $\|\tilde{S}^{(i)}\|_{1}=1$ for every $i\in\mathcal{A}$. Now, it remains to show that $\max_{i\in\mathcal{A}^{c}}Z_{3}^{(i)}<\gamma/2$ with high probability. Similar to $Z_{1}^{(i)}$, one can write:

[TABLE]

Given $X$, note that $T^{(i)}_{kl}$ is Gaussian with variance

[TABLE]

Also, $\|(\tilde{S}_{\mathcal{A}})_{:,l}\|_{2}^{2}\leq k_{j}$. Therefore, Lemma 3 can be used to bound (62) as follows:

[TABLE]

with probability of at least $1-2\exp(-d/{2})$. Similar to the arguments made for bounding $\max_{i\in\mathcal{A}^{c}}Z_{1}^{(i)}$, one can verify that

[TABLE]

Now, choosing

[TABLE]

for some arbitrary constant $c_{2}>1$ results in

[TABLE]

Therefore, $\max_{i\in\mathcal{A}^{c}}\|\tilde{S}^{(i)}\|_{1}<1$ and, hence, PDW succeeds with a probability that is lower bounded by (43).∎

IV-B Proof of Lemma 7:

In order to bound the estimation error, an upper bound on $\|E\|_{\infty}$ will be derived, conditioning on the success of the PDW approach. Note that $E_{\mathcal{A}^{c}}=0$ according to Lemma 10 and, hence, it suffices to bound $\|E_{\mathcal{A}}\|_{\infty}$. Again, due to Lemma 10, one can write:

[TABLE]

for $k=1,2,...,n_{j}$. For bounding $Z_{5}^{k}$, it can be argued similarly to (IV-A) that

[TABLE]

for some $c_{3}>0$ with probability of at least $1-2\exp(-({k_{j}n_{j}+c_{3}\log(\bar{n}+\bar{m})})/{2})$, where we have used the matrix norm properties and Lemma 3 with $t={k_{j}n_{j}+c_{3}\log(\bar{n}+\bar{m})}$ (note that $|I(\mathcal{A})|\leq k_{j}n_{j}$). Now, it remains to bound $\max_{k=1,...,n_{j}}Z_{4}^{k}$. This can be carried out similar to the previous arguments, i.e., by making use of (IV-A) and obtaining a sub-Gaussian parameter for $(\frac{1}{d}X_{\mathcal{A}}^{\top}X_{\mathcal{A}})^{-1}\frac{1}{d}X_{\mathcal{A}}^{\top}W_{:,k}$. For brevity, only the final key relation is stated below:

[TABLE]

where the last inequality is due to the assumption that $n_{j}\leq n_{\max}=O\left((\bar{n}+\bar{m})^{\alpha_{n}}\right)$ and $k_{j}\leq k_{\max}=O\left((\bar{n}+\bar{m})^{\alpha_{k}}\right)$. Now, setting

[TABLE]

for an arbitrary constant $c_{4}>1$, together with the inequality $\log r_{1}\leq\log(k_{j}D_{j})$, leads to

[TABLE]

with probability of at least

[TABLE]

Combining this inequality with (IV-B) results in the elementwise error bound (45) with probability of at least (44). This concludes the proof.∎

IV-C Proof of Theorem 2:

First, we present the sketch of the proof in a few steps:

We decompose the block-regularized problem (9) into $\bar{n}$ disjoint block-regularized subproblems defined in (41).

• 2.

For each of these subproblems, we consider the event that Lemmas 6 and 7 hold.

• 3.

We consider the intersection of these $\bar{n}$ events and show that, together with (20a) and (20b), they lead to the element-wise error (21) with probability of at least $1-\delta$.

Step 1: (9) can be rewritten as follows:

[TABLE]

The above optimization problem can be naturally decomposed into $\bar{n}$ disjoint block-regularized subproblems in the form of (41).

Step 2: Assume that (42b) and (42a) hold for every $1\leq j\leq\bar{n}$. Upon defining $\mathcal{T}_{j}$ as the event that Lemmas 6 and 7 hold, one can write:

[TABLE]

For every $1\leq j\leq\bar{n}$.

Step 3: Assume that $c_{1},c_{2},c_{4}>2$ and $c_{3}>1$. Consider the event $\mathcal{T}=\mathcal{T}_{1}\cap\mathcal{T}_{2}\cap\cdots\cap\mathcal{T}_{n}$. Based on (74) and Lemma 4, one can write:

[TABLE]

for some constants $K_{1},K_{2},K_{3},K_{4}$. One can easily verify that the following equalities are enough to guarantee that the right hand side of (75) is equal to $1-\delta$:

[TABLE]

Substituting (IV-C) in Lemmas 6 and 7 leads to two observations:

If $\lambda_{d}$ and $d$ satisfy (20a) and (20b), then they also satisfy (42a) and (42b).

• -

The parameter $g$ defined in (21) is greater than or equal to $g_{j}$ for every $j=1,...,\bar{n}$.

Therefore, (20a) and (20b) guarantee that: 1) $\hat{\Theta}$ is unique and does not have any false positive in its blocks, and 2) its element-wise error is upper bounded by (21). Now, it only remains to show that $\hat{\Theta}$ excludes false negatives (the blocks that are mistakenly estimated to have nonzero values). To this goal, it suffices to show that (20b) guarantees $g<t_{\min}$. Suppose that

[TABLE]

In what follows, we will show that $C_{\Theta}=O(1)$ is enough to have $g<t_{\min}$. The lower bound on $d$ in (20b) yields that

[TABLE]

for some constant $K$. Therefore,

[TABLE]

is enough to ensure $g<t_{\min}$. This completes the proof. ∎

IV-D Proof of Theorem 1:

The proof immediately follows from Theorem 2 by noting that $D=n_{\max}=1$, $\bar{n}=n$, and $\bar{m}=m$.∎

V Numerical Results

In this section, we illustrate the performance of the block-regularized estimator and compare it with its least-squares counterpart. We consider three case studies on synthetically generated systems, physical mass-spring networks, and multi-agent systems. The simulations are run on a laptop computer with an Intel Core i7 quad-core 2.50 GHz CPU and 16GB RAM. The reported results are for a serial implementation in MATLAB R2017b, and the function lasso is used to solve (9). Define the (block) mismatch error as the total number of false positives and false negatives in the (block) sparsity pattern of the estimator. Moreover, define relative number of sample trajectories (RST) as the number of sample trajectories normalized by the dimension of the system, and relative (block) mismatch error (RME) as the mismatch error normalized by total number of elements (blocks) in $\Theta$. To verify the developed theoretical results, $\lambda_{d}$ is set to

[TABLE]

in all of the experiments. Note that this choice of $\lambda_{d}$ does not require any additional fine-tuning.

V-A Case Study 1: Synthetically Generated Systems

Given the numbers $n$ and $w$, and for each instance of the problem, the state and input matrices are constructed as follows: The size of each block in $A$ and $B$ is set to 1. The diagonal elements of $A\in\mathbb{R}^{n\times n}$ and $B\in\mathbb{R}^{n\times n}$ are set to 1 (the dimensions of the inputs and states are chosen to be equal). The elements of the first $w$ upper and lower diagonals of $A$ and $B$ are set to $0.3$ or $-0.3$ with equal probability. Furthermore, at each row of $A$, another $w$ elements are randomly chosen from the elements not belonging to the first $w$ upper and lower diagonals and set to $0.3$ or $-0.3$ with equal probability. The structure of $[A\ \ B]$ is visualized in Figure 1 for $w=2,3,4$. Based on this procedure, the number of nonzero elements at each row of $[A\ \ B]$ varies between $3w+2$ and $5w+2$. We set $\Sigma_{u}=I$ and $\Sigma_{w}=0.5I$. The mutual incoherence property is satisfied for all constructed instances.

In the first set of experiments, we consider the mismatch error of $\hat{\Theta}$ with respect to the number of sample trajectories and for different system dimensions. The length of the time horizon $T$ is set to $3$. The results are illustrated in Figure 2a for $n+m$ equal to $200$, $600$, $1200$, and $2000$. In all of these test cases, $w$ is chosen in such a way that the number of nonzero elements in each column of $\Theta$ is between $(n+m)^{0.3}$ and $(n+m)^{0.4}$. It can be observed that as the dimension of the system increases, a higher number of sample trajectories is required to have a small mismatch error in the block-regularized estimator. Conversely, the required value of RST to achieve a small RME reduces as the dimension of the system grows. More precisely, RST should be at least $1.80$, $1.13$, $0.37$, and $0.20$ to guarantee $\text{RME}\leq 0.1\%$, when $m+n$ is equal to $200$, $600$, $1200$, and $2000$, respectively.

In the next set of experiments, we consider the mismatch error for different time horizons $T=3,4,...,7$, by fixing $m+n=600$ and $w=2$. As mentioned before, large values of $T$ tend to inflate the easily identifiable modes of the system and suppress the nearly hidden ones, thereby making it hard to obtain an accurate estimation of the parameters. It is pointed out that $\kappa(F_{T}F_{T}^{\top}+G_{T}G_{T}^{\top})$ is a good indicator of the gap between these modes. This relationship is clearly reflected in Figures 2b and 2c. As can be observed in Figure 2b, $330$ sample trajectories are enough to guarantee $\text{RME}\leq 0.1\%$ for $T=3$. However, for $T=7$, RME cannot be reduced below $0.42\%$ even with 1000 sample trajectories. To further elaborate on this dependency, Figure 2c is used to illustrate the value of $\kappa(F_{T}F_{T}^{\top}+G_{T}G_{T}^{\top})$ with respect to $T$ in a $\log$-$\log$ scale. One can easily verify that $\kappa(F_{T}F_{T}^{\top}+G_{T}G_{T}^{\top})$ associated with $T=7$ is $485$ times greater than this parameter for $T=3$.

Finally, we study the block-regularized estimator for different per-column numbers of nonzero elements in $\Theta$ and compare its accuracy to the least-squares estimator. Fixing $T=3$ and $m+n=600$, Figure 3a depicts the mismatch error of the block-regularized estimator when the maximum number of nonzero elements at each column of $\Theta$ ranges from $7$ (corresponding to $w=1$) to $27$ (corresponding to $w=5$). Not surprisingly, the required number of samples to achieve a small mismatch error increases as the number of nonzero elements in each column of $\Theta$ grows. On the other hand, the least-squares estimator is fully dense in all of these experiments, regardless of the number of sample trajectories. To have a better comparison between the two estimators, we consider the $2$-norm of the estimation errors normalized by the $2$-norm of $\Theta^{*}$, for different numbers of nonzero elements in each column of $\Theta^{*}$. As it is evident in Figure 3b, the block-regularized estimator significantly outperforms the least-squares one for any number of sample trajectories. Furthermore, the least-squares estimator is not defined for $d<600$.

V-B Case Study 2: Mass-Spring Systems

In this case study, we conduct simulations on mass-spring networks with different sizes to elucidate the performance of the block-regularized estimator on physical systems. Consider $N$ identical masses connected via a path of springs. Figure 4 exemplifies this system for $N=2$. The state-space equation of this system in continuous domain can be written as

[TABLE]

where $A_{c}\in\mathbb{R}^{2N\times 2N}$, $B_{c}\in\mathbb{R}^{2N\times N}$, and $x_{c}(t)$ consists of two parts: the first $N$ elements correspond to the locations of the masses while the second $N$ elements capture their velocities. With the unit masses and spring constants, we have

[TABLE]

where $S\in\mathbb{R}^{n\times n}$ is a tridiagonal matrix whose diagonal elements are set to $-2$ and its first upper and lower diagonal elements are set to $1$ [38]. This continuous system is discretized using the forward Euler method with sampling time of $0.2$ seconds. Similar to the previous case, $\Sigma_{u}$ and $\Sigma_{w}$ are set to $I$ and $0.5I$, respectively. Furthermore, $T$ is set to $3$ and $\lambda_{d}$ is chosen as (42a). The mutual incoherence property holds in all of these simulations. Notice that RST is equal to $3d/N$ in this case study, since $n=2N$ and $m=N$. Figure 5a depicts the required RST to achieve $\text{RME}\leq 0.1\%$ for different numbers of masses. If $N=30$, the minimum number of sample trajectories to guarantee $\text{RME}\leq 0.1\%$ is $5$ times higher than the dimension of the system, whereas this number is dropped to $0.02$ for $N=500$. Furthermore, Figure 5b shows the normalized estimation errors of the block-regularized and least-squares estimators for $N=300$ and $T=3$. Similar to the previous case study, the proposed estimator results in significantly smaller estimation errors in medium- and large-sampling regimes, while it remains as the only viable estimator when the number of available sample trajectories is smaller than $900$.

V-C Case Study 3: Switching Networks

In this case study, we study a network of multi-agent systems that are interconnected through a switching information exchange topology. Recently, a special attention has been devoted to multi-agent systems with a time-varying network topology; in many communication networks, each sensor has access only to the information of its neighbors. Therefore, when the location of these sensors changes over time, so does the topology of the interconnecting links [39]. The dwell time is defined as the time interval in which the network topology is unchanged. The goal is to identify the structure of the network within the dwell time. The state-space equation of agent $i$ admits the following general form:

[TABLE]

where, as before, $A^{(i,j)}\in\mathbb{R}^{n_{i}\times n_{i}}$ and $B^{(i,j)}\in\mathbb{R}^{n_{i}\times m_{i}}$ are the $(i,j)^{\text{th}}$ blocks of $A$ and $B$. Furthermore, $\mathcal{N}_{x}(i)$ and $\mathcal{N}_{u}(i)$ are the sets of neighbors of agent $i$ whose respective state and input actions affect the state of agent $i$.

We consider 200 agents connected through a randomly generated sparse network. In particular, we assume that each agent is connected to 5 other agents. If $j\in\mathcal{N}_{x}(i)$ or $j\in\mathcal{N}_{u}(i)$, then each element of $A^{(i,j)}$ or $B^{(i,j)}$ is randomly selected from $[-0.4\ -0.3]\cup[0.3\ 0.4]$. The behavior of the proposed block-regularized estimator will be examined for different dimensions of the agents. In particular, $(n_{i},m_{i})$ will be chosen from $\{(5,5),(8,8),(11,11)\}$. This entails that $D\in\{25,64,121\}$ and $(n,m)\in\{(1000,1000),(1600,1600),(2200,2200)\}$. Similar to the previous case study, $T$ is set to $3$ and the system is discretized using the forward Euler method with the sampling time of $0.2$ seconds. This implies that each sample trajectory is collected within $0.6$ seconds. The respective number of block mismatch and normalized estimation errors is depicted in Figures 6a and 6b with respect to the dwell time. It can be seen that as the dwell time becomes longer, a larger number of sample trajectories can be obtained, which in turn results in smaller block mismatch and normalized estimation errors. Furthermore, as the size of the blocks grows, an accurate identification of the system parameters requires more sample trajectories. In particular $195$, $278$, and $434$ sample trajectories are needed to achieve $\text{RME}\leq 0.1\%$ when $D$ is equal to $25$, $64$, and $121$, respectively. In contrast, the least-squares estimator requires at least $2000$, $3200$, and $4400$ sample trajectories to be uniquely defined.

VI Conclusion

We consider the problem of identifying the parameters of linear time-invariant (LTI) systems. In many real-world problems, the state-space equation describing the evolution of the system admits a block-sparse representation due to localized or internally limited interactions of its states and inputs. In this work, we leverage this property and introduce a block-regularized estimator to identify the sparse representation of the system. Using modern high-dimensional statistics, we derive sharp non-asymptotic bounds on the minimum number of input-state data samples to guarantee a small element-wise estimation error. In particular, we show that the number of available sample trajectories can be significantly smaller than the system dimension and yet, the proposed block-regularized estimator can correctly recover the block-sparsity of the state and input matrices and result in a small element-wise error. Through different case studies on synthetically generated systems, mass-spring networks, and multi-agent systems, we demonstrate substantial improvements in the accuracy of the proposed estimator, compared to its well-known least-squares counterpart.