Missing Data in Sparse Transition Matrix Estimation for Sub-Gaussian Vector Autoregressive Processes

TL;DR

This paper develops new theoretical tools for estimating sparse transition matrices in high-dimensional vector autoregressive models with missing data, addressing challenges posed by dependent covariates and partial observations.

Contribution

It introduces a novel concentration inequality and a new interaction measure to handle missing data in dependent, high-dimensional time series models.

Findings

Established a concentration result for innovation noise with convex concentration property

Derived bounds for transition matrix estimation under missing data

Provided theoretical guarantees for sparse VAR models with partial observations

Abstract

High-dimensional time series data exist in numerous areas such as finance, genomics, healthcare, and neuroscience. An unavoidable aspect of all such datasets is missing data, and dealing with this issue has been an important focus in statistics, control, and machine learning. In this work, we consider a high-dimensional estimation problem where a dynamical system, governed by a stable vector autoregressive model, is randomly and only partially observed at each time point. Our task amounts to estimating the transition matrix, which is assumed to be sparse. In such a scenario, where covariates are highly interdependent and partially missing, new theoretical challenges arise. While transition matrix estimation in vector autoregressive models has been studied previously, the missing data scenario requires separate efforts. Moreover, while transition matrix estimation can be studied from a high-dimensional sparse linear regression perspective, the covariates are highly dependent and existing results on regularized estimation with missing data from i.i.d.~covariates are not applicable. At the heart of our analysis lies 1) a novel concentration result when the innovation noise satisfies the convex concentration property, as well as 2) a new quantity for characterizing the interactions of the time-varying observation process with the underlying dynamical system.

Equations

$$w_{t+1}=B_0{w}_t+{ϵ}_t,t=0,\dots ,n-1,$$

$$\underset{\mathcal{Y}}{\underbrace{\left[\begin{array}{c}{w}_1\cdots w_n\end{array}\right]}}=B_0\underset{\mathcal{X}}{\underbrace{\left[\begin{array}{c}{w}_0\cdots w_{n-1}\end{array}\right]}}+\underset{\mathcal{E}}{\underbrace{\left[\begin{array}{c}{ϵ}_0\cdots {ϵ}_{n-1}\end{array}\right]}}$$

$$(\overline{w}{}_t{)}_i=\{\begin{array}{ll}(w_t{)}_i& \text{with probability}1-\delta \\ 0& \text{with probability}\delta ,\end{array}$$

$${\mathrm{argmin}}_{\Vert B{\Vert }_1\le b_0\sqrt{k}}\frac{1}{n}\Vert B\overline{\mathcal{X}}-\overline{\mathcal{Y}}{\Vert }_F^2-\frac{\delta }{n}\Vert B\overline{\mathcal{D}}{\Vert }_F^2+(1-\delta {)}^2{\lambda }_n\Vert B{\Vert }_1$$

$${\mathrm{argmin}}_{\Vert B{\Vert }_1\le \Vert B_0{\Vert }_1}\frac{1}{n}\Vert B\overline{\mathcal{X}}-\overline{\mathcal{Y}}{\Vert }_F^2-\frac{\delta }{n}\Vert B\overline{\mathcal{D}}{\Vert }_F^2.$$

$${\mathrm{argmin}}_{\beta }\frac{1}{n}\Vert X\beta -y{\Vert }_2^2+\lambda \Vert \beta {\Vert }_1$$

$$\rho (B)=\max \{|\lambda |:\lambda \in \sigma (B)\}$$

$$|||B||{|}_{\iota \to \varrho }:=\underset{x\ne 0}{\sup }\frac{\Vert Bx{\Vert }_{\varrho }}{\Vert x{\Vert }_{\iota }}.$$

$$\det (I-B_{0}z)\ne 0\text{for all}|z|\le 1,$$

$${\vartheta }_0(B)$$

$${\vartheta }_1(B)$$

$${\vartheta }_2(B)$$

$${\vartheta }_2(B)\le {\vartheta }_1(B)\le \sqrt{2k}{\vartheta }_2(B).$$

$${\theta }_0:=\frac{{\vartheta }_2(B_0{)}^2}{{\vartheta }_1(B_0{)}^2}\in [\frac{1}{2k},1].$$

$${\kappa }_{ϵ}$$

$${\kappa }_0$$

$$h:=\frac{b_0}{7(|||B_0||{|}_{2\to \infty }^2+1)\sqrt{k}},\zeta :=\frac{1+\delta {\theta }_{0}k}{hk}.$$

$$\frac{\Vert B_0{\Vert }_F}{(|||B_0||{|}_{2\to \infty }^2+1)\sqrt{k}}$$

$$\sqrt{\frac{n}{\log p}}\ge \frac{{\kappa }_{ϵ}{\kappa }_0\zeta }{(1-\delta {)}^2}\frac{1}{(\frac{1}{27}\zeta -\delta {\theta }_0{)}^2}$$

$$\Phi :=\frac{c_0{b}_0\sqrt{k}}{7}\frac{{\kappa }_{ϵ}{\kappa }_0\zeta }{(1-\delta {)}^2}\sqrt{\frac{\log p}{n}}.$$

$${\phi }_0$$

$$\Vert \hat{B}-B_0{\Vert }_F\le 2\sqrt{k}{\phi }_0{\lambda }_n,\Vert \hat{B}-B_0{\Vert }_1\le 16k{\phi }_0{\lambda }_n$$

$$|\mathrm{supp}(\tilde{B})\setminus \mathrm{supp}(B_0)|\le 112k{\phi }_0.$$

$${\vartheta }_0(B)\le 1+|||B||{|}_2\le 1+\frac{|||B||{|}_1+|||B||{|}_{\infty }}{2}.$$

$${\vartheta }_1(B)\le \frac{1}{1-\rho (B)}|||R||{|}_2|||R^{-1}||{|}_2$$

$${\vartheta }_2(B)\ge \underset{|z|=1}{\max }|||(I-Bz)||{|}_{1\to 2}^{-1}\ge (1+|||B||{|}_{1\to 2}{)}^{-1}.$$

$$\sqrt{{\theta }_0}=\frac{{\vartheta }_2(B_0)}{{\vartheta }_1(B_0)}\ge \frac{1-\rho (B)}{(1+|||B||{|}_{1\to 2})|||R||{|}_2|||R^{-1}||{|}_2}$$

$$\sqrt{{\kappa }_0}={\vartheta }_0(B){\vartheta }_1(B)\le \frac{|||R||{|}_2|||R^{-1}||{|}_2}{1-\rho (B)}(1+\frac{|||B||{|}_1+|||B||{|}_{\infty }}{2}).$$

$$B_0=\left[\begin{array}{cc}a& \frac{1}{a}\\ 0& a\end{array}\right]$$

$$w_t=B_1{w}_{t-1}+B_2{w}_{t-2}+\dots +B_d{w}_{t-d}+{ϵ}_t,$$

$$\left[\begin{array}{c}{w}_t\\ w_{t-1}\\ ⋮\\ w_{t-d+1}\end{array}\right]=\underset{B}{\underbrace{\left[\begin{array}{ccccc}{B}_1& B_2& \cdots & B_{d-1}& B_d\\ I_p& 0& \cdots & 0& 0\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& \cdots & I_p& 0\end{array}\right]}}\left[\begin{array}{c}{w}_{t-1}\\ w_{t-2}\\ ⋮\\ w_{t-d}\end{array}\right]+\left[\begin{array}{c}{ϵ}_t\\ 0\\ ⋮\\ 0\end{array}\right].$$

Full text

Introduction

Consider a $p$-dimensional covariance-stationary vector autoregressive model of lag one, namely a $\textup{VAR}(1)$, as

[TABLE]

where $B_{0}\in\mathbb{R}^{p\times p}$ is the corresponding transition matrix, and each $\epsilon_{t}$ is a $p$-dimensional vector of innovations, with zero mean and covariance $\Sigma_{\epsilon}$, that is temporally uncorrelated with other noise vectors. The goal is to estimate $B_{0}$ from partial observations of entries of $w_{0},\ldots,w_{n}$, given prior knowledge on $B_{0}$ being sparse. Concatenating all the vectors in (1) as columns of matrices yields

[TABLE]

where each of the brackets represent a $(p\times n)$-dimensional matrix. The available information from (1) is in the form of entries in $W=\begin{bmatrix}w_{0}\;\cdots\;w_{n}\end{bmatrix}$ that are missing according to i.i.d. Bernoulli random variables with probability $0\leq\delta<1$. Consider a new process $\{\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{t}\}$ where, for any $i=1,\ldots,p$,

[TABLE]

and observation is independent for different $i=1,\ldots,p$ and different $t=1,\ldots,n$. We use the bar notation for other objects constructed from $\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{0},\ldots,\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{n}$. For example, $\mkern 1.0mu\overline{\mkern-1.0muW\mkern-1.0mu}\mkern 1.0mu=\begin{bmatrix}\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{0}\;\cdots\;\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{n}\end{bmatrix}$, $\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu=\begin{bmatrix}\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{0}\;\cdots\;\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{n-1}\end{bmatrix}$, and so forth. For simplicity, we consider a centered process, hence $w_{0}=0$.

We handle the missing data by modifying the LASSO [Tib96] using population information on the observation pattern, as described in the Appendix. More specifically, similar to [LW12], we solve either of the two following constrained versions of this program: either

[TABLE]

where $\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{D}\mkern-1.0mu}\mkern 1.0mu=({\operatorname{diag}}(\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu^{\prime}))^{1/2}\in{\mathbb{R}^{p\times p}}$ is a diagonal matrix of sample autocovariances for each of the $p$ covariates, $k$ is the number of nonzero entries of $B_{0}$, i.e., $k=\|B_{0}\|_{0}$, $b_{0}$ is any value at least equal to $\|B_{0}\|_{F}$, and $\lambda_{n}$ is the regularization parameter that will be chosen according to the parameters of the problem, or

[TABLE]

Note that, with a possibly non-convex quadratic optimization program, we need a suitably constrained feasible set to avoid an unbounded optimization problem and hope for recovering the target model, hence the constraints in (3) and (4).

Through employing a well-known machinery for the analysis of LASSO (summarized as Theorem 6 in the Appendix), and by developing new concentration results for the random processes of interest in this work (Proposition 5), we provide guarantees on the $\ell_{1}$- and $\ell_{2}$-norm estimation errors for (3) and (4) in Theorem 1. Before stating the main result, we discuss the prior art and how our setup leads to new challenges. We then discuss certain characteristics of the processes in (1) and (2) that are used in our guarantees, in Section 1.3. More specifically, we introduce a new quantity, namely $\vartheta_{2}(B_{0})$ in (8), which is used in characterizing the interplay between the dependence among covariates and the difficulty of recovery from partial information. We elaborate on these characteristics in Section 3. Section 4 provides a sketch of the proof for providing error bounds on LASSO and its variants, and is similar to many other works on LASSO in the literature. Section 5 contains our main contribution on establishing the required concentration inequalities for providing error bounds for LASSO through concentration of sub-Gaussian quadratic forms.

Prior Art and New Challenges

Estimators (3) and (4) can be seen as modifications of the LASSO [Tib96], with constraints that help remedy the possible non-convexity of the estimator. Similar estimators have been considered in the literature for several sparse regression tasks and a similar framework has been used to analyze such estimators; e.g., in [LW12, BRT09, RWY10, BM15], and many more. However, the distinguishing aspect of different works in this area is the difference in the data generation processes and the required concentration analysis. In a simple data generation scheme as $y=X\beta_{0}+\epsilon$ where $\epsilon\sim\mathcal{N}(0,I)$ and $X$ has i.i.d. random entries drawn from $\mathcal{N}(0,1)$ independently from $\epsilon$, the analysis of LASSO,

[TABLE]

boils down to understanding the spectrum of the random matrix $X$. A more complicated case of correlated Gaussian designs is considered in [RWY10]. More involved data generation scenarios require more involved probabilistic arguments to establish the conditions that guarantee (near-)optimality of $\beta_{0}$ for (5). For example, [BM15] extends the above to transition matrix estimation in Gaussian vector autoregressive models. Authors in [LW12] extend (5) to the case of non-convex quadratic optimization programs when $X$ has one of several interesting dependency patterns.

We note that our focus is different from [LW12, Corollary 4] which considers sparse regression with missing data when the design matrix is generated by an autoregressive process with known transition matrix $A$ satisfying $|\kern-1.07639pt|\kern-1.07639pt|A|\kern-1.07639pt|\kern-1.07639pt|_{2}<1$. In our case, the interactions among covariates depend on the unknown transition matrix, making the problem even more challenging.

Assuming Gaussian innovations, $\epsilon_{0},\ldots,\epsilon_{n-1}$, makes $\{w_{t}\}$ of (1) a Gaussian process. However, the partially observed process $\{\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{t}\}$ will no longer be a Gaussian process. Nonetheless, $\{\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{t}\}$ belongs to the family of sub-Gaussian processes for which many properties are known. In this work, we consider a subset of sub-Gaussian processes for the innovations: those with the convex concentration property in Definition 3.

Matrix Quantities

We now review some important quantities associated to the transition matrix of interest, $B_{0}$. Consider $B_{0}$ as the adjacency matrix of a weighted directed graph on $p$ nodes as in the left panel of Figure 1. It is intuitive that not only the number of edges in such graph but also the configuration of edges, e.g., the degree distribution, plays an important role in any inverse problem for identifying such graph. In the following, we present important quantities associated with $B_{0}$ that help in reflecting nuances of the corresponding graph. From a dynamical systems point of view, we review the notions of the spectral radius and the spectral norm, as well as new quantities presented in Section 1.3, which allow us to capture the stability of the autoregressive process as well as the information content of our observations from this system.

Denote by $\mathbb{C}^{p\times p}$ the space of all $p$ by $p$ complex matrices. The spectral radius of $B\in\mathbb{C}^{p\times p}$ is the non-negative real number

[TABLE]

where $|\cdot|$ denotes the magnitude and $\sigma(B)$ is the set of all eigenvalues of $B$. A norm on $\mathbb{C}^{p\times p}$ is called a matrix norm if it satisfies the submultiplicative property as $|\kern-1.07639pt|\kern-1.07639pt|AB|\kern-1.07639pt|\kern-1.07639pt|\leq|\kern-1.07639pt|\kern-1.07639pt|A|\kern-1.07639pt|\kern-1.07639pt||\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|$ for all $A,B\in\mathbb{C}^{p\times p}$. If $|\kern-1.07639pt|\kern-1.07639pt|\cdot|\kern-1.07639pt|\kern-1.07639pt|$ is any matrix norm, then $\rho(B)\leq|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|$. For any $\varrho,\iota\geq 1$, consider the corresponding $\ell_{\varrho}$ and $\ell_{\iota}$ vector norms. The corresponding induced operator norm is then defined as

[TABLE]

For example, $|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|_{1\to 2}$ is the maximum $\ell_{2}$ norm among the columns of $B$, and $|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|_{2\to\infty}$ is the maximum $\ell_{2}$ norm among the rows of $B$. When $\iota=\varrho$, we simply denote $|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|_{\iota}\coloneqq|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|_{\iota\to\iota}$. For example, $|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|_{1}=\max_{1\leq j\leq p}\sum_{i=1}^{p}|B_{ij}|$ and $|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|_{\infty}=\max_{1\leq i\leq p}\sum_{j=1}^{p}|B_{ij}|$. When $\iota=2$, $|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|_{2}=\sqrt{\Lambda_{\max}(B^{\prime}B)}$ is the spectral norm which is also simply referred to as the operator norm. Notice that $|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|_{2}\neq\rho(B)$ in general. The extension of vector $\ell_{\iota}$ norms to matrices is denoted by $\|B\|_{\iota}\coloneqq\|{\operatorname{vec}}(B)\|_{\iota}$. Finally, we abuse notation to denote by $\|B\|_{0}$ the number of nonzero elements in the vector or matrix input. The interested reader is referred to Section 5.6 of [HJ90] for a more comprehensive account of these matrix quantities.

Going back to the directed graph interpretation of $B$, we can view $|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|_{1\to 2}$ as the maximum energy that any one node can exert on other nodes which is related to the maximum out-degree of nodes, while $|\kern-1.07639pt|\kern-1.07639pt|B|\kern-1.07639pt|\kern-1.07639pt|_{2\to\infty}$ can be viewed as the maximum energy that is being exerted upon each node and is related to the maximum in-degree among nodes.

Salient Characteristics of the Model

The autoregressive process in (1) is called stable if and only if $\rho(B_{0})<1$: all eigenvalues of $B_{0}$ have modulus less than one. This is equivalent to

[TABLE]

where $z$ is a complex-valued scalar variable. In such case, we define three main quantities

[TABLE]

The first two quantities are related to the least and the largest singular values of the transfer function on the unit circle. In other words, for all $z$ with $|z|=1$, they quantify the least and the largest values of $\|(I-Bz)^{-1}u\|_{2}$ when $\|u\|_{2}=1$. The third quantity, on the other hand, characterizes the largest value of $\|(I-Bz)^{-1}u\|_{2}$ when $\|u\|_{1}=1$, for all $z$ with $|z|=1$. Moreover, for $k=\|B\|_{0}$, we have (see Section 3.2)

[TABLE]

While (6) and (7) have been considered in other works on autoregressive models (e.g., see $\mu_{\max}$ and $\mu_{\min}$ in Equation (2.6) of [BM15]), the definition of (8) in the context of autoregressive model estimation is, to the best of our knowledge, new and motivated by the missing data setup.

Main Results

In this section, we state our main result followed by a discussion on the main quantities, a sketch of the proof, and a list of ingredients for this proof that we establish in the subsequent sections. In essence, we would show that the error scales with

[TABLE]

We define a few more quantities to simplify the presentation of our main result. First, consider an innovation condition number defined as

[TABLE]

where $c_{a}$ is a global constant and $c_{\epsilon}$ is a function of the innovation process $\epsilon$ and will be defined later. We also consider a a transition condition number defined as

[TABLE]

Finally, the nonzero pattern of $B_{0}$ and the quality of our choice for $b_{0}$ can be measured through the followings:

[TABLE]

Note that $\|B_{0}\|_{F}\leq\sqrt{k}|\kern-1.07639pt|\kern-1.07639pt|B_{0}|\kern-1.07639pt|\kern-1.07639pt|_{2\to\infty}\leq\sqrt{k}|\kern-1.07639pt|\kern-1.07639pt|B_{0}|\kern-1.07639pt|\kern-1.07639pt|_{2\to\infty}^{2}/\rho(B_{0})$, which helps in understanding $h$ through

[TABLE]

as $b_{0}$ is chosen to be at least $\|B_{0}\|_{F}$.

Theorem 1** (main result).**

Consider the $p$-dimensional autoregressive process in (1) satisfying $\rho(B_{0})<1$. Suppose $\|B_{0}\|_{0}=k$ and the innovations are temporally uncorrelated with zero mean and a positive definite covariance matrix $\Sigma_{\epsilon}$, and satisfy the convex concentration property (Definition 3) with constant $c_{\epsilon}$. Suppose we have partially observed the process, with missing probability $\delta$, for time length $n$ satisfying

[TABLE]

while $\zeta>27\delta\theta_{0}$. Define

[TABLE]

Consider any $b_{0}\geq\|B_{0}\|_{F}$ and any $\lambda_{n}$ satisfying $\lambda_{n}\geq\frac{2\Phi}{\varphi_{0}}$ where

[TABLE]

and $c_{0}$ and $c_{1}$ are universal constants. Then, with probability at least $1-10p^{-1}$, for any optimal $\widehat{B}$ in (3) we have

[TABLE]

and, for \tilde{B}\coloneqq\big{[}\widehat{B}_{ij}\mathbf{1}_{|\widehat{B}_{ij}|>\lambda_{n}}\big{]}_{i,j=1,\ldots,p} we have

[TABLE]

The same bounds, where $\lambda_{n}$ is replaced by $\frac{2\Phi}{\varphi_{0}}$, apply to (4).

As mentioned before, there is a well-developed theory for providing guarantees on the performance of the LASSO and its variants which we summarize as Theorem 6 in the Appendix. This framework requires establishing certain concentration properties for the underlying data generation process. We provide the required concentration results in Section 5 and combine them with the framework of Theorem 6 to prove Theorem 1. In the following, we first provide further intuition into the results of Theorem 1 in Section 3. We then elaborate on the required conditions for proving Theorem 1 and motivate our concentration result in Section 4.

Remarks on the Main Quantities

The quantities appearing in Theorem 1 worth further discussion. In this section, we provide further details on different quantities we defined in relation to the transition matrix of interest, $B_{0}$.

The Support

First, as we will see in the concentration result, $\vartheta_{2}(B_{0})$ appears because of the missing data setup; specifically, due to the term $-\|B\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{D}\mkern-1.0mu}\mkern 1.0mu\|_{F}^{2}$ in (3) and (4). Intuitively, we expect that the support of $B_{0}$, and how each covariate affects the value of other covariates in the next time step (see Figure 1), should play an important role in our ability in recovery from missing data. For example, if $B_{0}$ is diagonal, then covariates are temporally uncorrelated (do not directly affect each other over time) and the entries of $B_{0}$ have to be estimated independently. On the other hand, for more distributed supports of $B_{0}$, we experience two competing phenomena:

• •

when each covariate is influenced by many covariates from the previous time point, residuals between $B_{0}\mathcal{X}$ and $B_{0}\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu$ are generally smaller because the value of missing covariates play less of a role, so recovery is more robust to missing data. This is captured by $\theta_{0}$ (defined in (10)) in our results.

• •

higher dependence among covariates makes observations more highly correlated and the resulting inverse problem becomes more ill-posed even when we have complete data. This is captured by $\kappa_{0}$ (defined in (11)) in our results.

In short, not all $k$-sparse $B_{0}$ are equally easy or difficult to infer from incomplete data. For example, if only one column of $B_{0}$ is nonzero (in-star graph in Figure 2), then one element of $w_{t}$ is influenced by the previous realizations of the process, while the other covariates are not. If only one row of $B_{0}$ is nonzero (out-star graph in Figure 2), then all covariates are being influenced by the same single covariate and there are no other influences. Finally, if $B_{0}$ is nonzero on a single off-diagonal, then the $i$-th covariate is only influencing covariate $i+1$, for $i=1,\ldots,p-1$, corresponding to a chain graph representation of influence structure.

Dimension-independence

Denote the set of nonzero rows of $B$ by $J_{r}\subseteq\{1,\ldots,p\}$ and the set of its nonzero columns by $J_{c}\subseteq\{1,\ldots,p\}$. For $k=\|B\|_{0}$, it is easy to see that $|J_{r}|\leq k$ and $|J_{c}|\leq k$, with $|\cdot|$ denoting the size of the set. Moreover, for $J\coloneqq J_{r}\cup J_{c}$, all of the nonzero entries of $B$ are in a principal submatrix indexed by $J$. Therefore, for any integer value $t\geq 1$, all of the nonzero entries of $B^{t}$ are in the same principal sub-matrix indexed by $J$. Considering the Neumann series $(I-A)^{-1}=\sum_{t=0}^{\infty}A^{t}$ when $\rho(A)<1$, the above implies that $\vartheta_{0}(B)$, $\vartheta_{1}(B)$, and $\vartheta_{2}(B)$, are only concerned with the smallest principal submatrix of $B$ containing all of its nonzero entries and are independent of the dimension of $B$: embedding $B$ into a larger zero matrix does not change these values, as desired.

It is worth mentioning that the same conclusion, of independence from the ambient dimension, cannot be made about the quantity $\mathcal{M}(f_{w},s)$ used in [BM15] as the innovations could make the time series fully supported over all entries.

Bounds

While $B$ and the quantities in (6), (7), and (8), do not directly scale with each other, they have a close relationship that can be used in better understanding the main theorem.

Lemma 2** (Proposition 2.2 in [BM15]).**

Suppose $\det(I-Bz)\neq 0$ for all $|z|\leq 1$. Then,

[TABLE]

Moreover, if $B$ is diagonalizable, then

[TABLE]

where the columns of $R$ are the eigenvectors of $B$.

Moreover, submultiplicativity of induced operator norms provides

[TABLE]

Therefore, when $B$ is diagonalizable,

[TABLE]

and

[TABLE]

A Restrictive Assumption We Avoid

In this paper, we only assume stability, i.e., $\rho(B_{0})<1$. This assumption is milder than the more stringent condition $|\kern-1.07639pt|\kern-1.07639pt|B_{0}|\kern-1.07639pt|\kern-1.07639pt|_{2}<1$ prevalent in the literature. Only requiring the milder assumption used in this paper has important practical consequences. While $|\kern-1.07639pt|\kern-1.07639pt|B_{0}|\kern-1.07639pt|\kern-1.07639pt|_{2}<1$ implies $\rho(B_{0})<1$ (hence the stability of the corresponding autoregressive process), $|\kern-1.07639pt|\kern-1.07639pt|B_{0}|\kern-1.07639pt|\kern-1.07639pt|_{2}<1$ is necessary only when $B_{0}$ is symmetric. In other words, an asymmetric matrix $B_{0}$ with $|\kern-1.07639pt|\kern-1.07639pt|B_{0}|\kern-1.07639pt|\kern-1.07639pt|_{2}\geq 1$ can correspond to a stable autoregressive process; e.g., see Lemma E.1 in [BM15]. For example, for some $0<a<1$, the matrix

[TABLE]

has eigenvalues equal to $a$, hence a spectral radius of $a<1$, but an operator norm that is slightly larger than $\frac{1}{a}>1$. Not assuming a spectral norm bound on the transition matrix becomes important in the study of vector autoregressive processes with a lag larger than one, defined as

[TABLE]

where $d\geq 1$ is the lag. It is easy to see that the above can be reformulated as a vector autoregressive process with lag one, as

[TABLE]

Lemma E.2 in [BM15] establishes the fact that $d>1$ implies $|\kern-1.07639pt|\kern-1.07639pt|{B}|\kern-1.07639pt|\kern-1.07639pt|_{2}\geq 1$, even when $\rho({B})<1$, illustrating the restrictiveness of operator norm bound assumptions.

Estimation Error for Non-convex LASSO

Both (3) or (4) can be viewed as constrained quadratic optimization programs,

[TABLE]

for corresponding choices of the constraint set $\mathcal{B}\subset\mathbb{R}^{p\times p}$ and regularization parameter $\lambda$, where

[TABLE]

In this work, we are not concerned with the possible non-convexity of these estimators from a computational point of view and focus on the statistical performance. Nonetheless, simple algorithms such as variants of projected gradient descent can be used for convergence to a small neighborhood of the set of all global minimizers, similar to [LW12]. We postpone such convergence guarantees to future work.

LASSO [Tib96] and its variants have been used and studied extensively in the literature. We specifically use a popular approach for providing guarantees on the estimation performance of LASSO and its variants presented in [BRT09]. We present a version of this result, tailored to norm-constrained $\ell_{1}$-regularized non-convex quadratic optimization, as Theorem 6 in the Appendix. Theorem 6 is essentially the same in any work on guarantees for LASSO and its variants, but lumps all the mechanical, and now well-known, parts of the process in one theorem and is discussed for clarity of our exposition. For example, the theorem can be seen as an extension of [WTL16, Theorem 5] and [BM15, Proposition 4.1] for non-convex LASSO and an extension of [LW12, Theorem 1] for transition matrix estimation in multivariate time series.

To provide $\ell_{1}$ and $\ell_{2}$ norm error bounds for such estimates, following the framework developed in [BRT09], we need to establish the so-called lower restricted eigenvalue condition, stated equivalently [LW12] as

[TABLE]

as well as a deviation bound,

[TABLE]

where $c$ depends on the parameters of the problem instance.

Since $Q$ and $L$ come from samples generated by the partial observation of a vector autoregressive model, they are random objects and reasonable values of $\alpha_{\textup{low}}$, $\tau_{\textup{low}}$, and $c$, in (14) and (15), may be used only with high probability. Therefore, we use relevant concentration results to establish these bounds with high probability.

In the following, we expand the conditions in (14) and (15) and represent them as simple functions of the autoregressive process, which will then be bounded in Section 5 using results on concentration of sub-Gaussian quadratic forms. Let us fix some notation first. For a $p$-dimensional discrete-time, centered, covariance-stationary (wide-sense stationary) process $\{w_{t}\}$, denote the autocovariance function by

[TABLE]

For a matrix $A$, the transpose is denoted by $A^{\prime}$. Denote by $\odot$ and $\oslash$ the Hadamard (element-wise) product and division respectively, and by $\otimes$ the Kronecker product. The covariance matrix for the Bernoulli mask characterized in (2) is given by

[TABLE]

so that $Q=\frac{1}{n}\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu^{\prime}\oslash P$.

Restricted Eigenvalue Condition

Observe that $\mathbb{E}Q=\Gamma_{w}(0)=\Gamma_{\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu}(0)\oslash P$, which gives

[TABLE]

Then, bounding $|v^{\prime}(Q-\mathbb{E}Q)v|$, for all $v\in\mathbb{R}^{p}$, allows for establishing (14) through the application of the triangle inequality. Suppose we established the following condition for a fixed value of $s$ which will be determined later:

1. (C1)

For any fixed $v\in\mathbb{R}^{p}$ with $\|v\|_{0}\leq 2s$ and $\|v\|_{2}=1$, there exists $\eta(s)$ such that $|v^{\prime}(Q-\mathbb{E}Q)v|\leq\eta(s)$ with probability at least $1-p_{1}(s)$.

Then, such concentration can be stated over the set of $2s$-sparse vectors using a discretization argument, as in Lemma F.2 of [BM15], followed by a simple argument that relates the set of sparse vectors to those with a bounded $\ell_{1}$ norm, as in Lemma 12 of [LW12]. As the above calculations depend on the free parameter $s$, it should be chosen in a way that makes $\eta(s)$ as small as possible while maintaining the probability for (14), which depends on $p_{1}(s)$ and the union bound in the discretization step, at a desired level. We specify our choice of $s$ for the proof of Theorem 1 right after the statement of Theorem 6 in the Appendix.

Deviation Bound

The matrix of interest in (15) is given by

[TABLE]

where we used the fact that $B_{0}\Gamma_{w}(0)=\Gamma_{w}(1)^{\prime}$ and $\Gamma_{w}(1)=\frac{1}{(1-\delta)^{2}}\Gamma_{\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu}(1)$. The first assertion considers full information and is related to the interaction of $\{w_{t}\}$ and $\{\epsilon_{t}\}$ processes. In fact, using the original process in (1) we get

[TABLE]

which is zero in our setup. For clarity, we state (15) as another condition:

1. (C2)

There exists $\varphi>0$ such that $\|B_{0}Q-L\|_{\infty}\leq\varphi$ with probability at least $1-p_{2}$.

Therefore, to derive the desired bounds in (14) and (15) and provide $\ell_{1}$ and $\ell_{2}$ norm error bounds for (3) and (4), we can establish (C1) and (C2); a complete description of this procedure is stated as Theorem 6 in the Appendix. This amounts to computing concentration bounds on the four terms in (16) and (17). To that end, we rewrite our process in matrix form and leverage classical linear time invariant dynamical systems to establish several quantities that characterize the process in (1) and will appear in those concentration bounds. These relationships are summarized in Lemma 4 and lead to the main concentration result given in Proposition 5.

Concentration of Sub-Gaussian Quadratic Forms

As mentioned before, establishing either of the conditions (16) and (17) relies on certain concentration properties for the underlying data generation process that defines $Q$ and $L$. In the following, we make this relationship concrete and provide the main concentration result in Proposition 5.

To establish (C1) for (16) (to get (14)), we are interested in the concentration of $v^{\prime}\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu^{\prime}v$ and $v^{\prime}(\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu^{\prime}\odot I)v$ around their mean, for any fixed $v$. In the following, we manipulate these quantities into convex quadratic forms in terms of the noise vector

[TABLE]

Define $I_{\Omega}\in\{0,1\}^{pn\times pn}$ as the diagonal matrix whose $(pt+j)$-th diagonal entry is one if $(w_{t})_{j}$ is observed and zero otherwise, for $t=0,\ldots,n-1$ and $j=1,\ldots,p$. Moreover, define

[TABLE]

which is a block-Toeplitz matrix. Then, $v^{\prime}\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu^{\prime}v=\mathfrak{e}_{n}^{\prime}\Psi_{(1)}^{\prime}\Psi_{(1)}\mathfrak{e}_{n}$ and $v^{\prime}(\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu^{\prime}\odot I)v=\mathfrak{e}_{n}^{\prime}\Psi_{(2)}^{\prime}\Psi_{(2)}\mathfrak{e}_{n}$ where

[TABLE]

The latter is because of the following,

[TABLE]

The concentration of the above two quadratic forms, in $\Psi_{(1)}$ and $\Psi_{(2)}$, can be studied when we assume the so-called convex concentration property on noise vectors $\epsilon_{t}$, for $t=0,\ldots,n-1$, or equivalently on the noise vector $\mathfrak{e}_{n}$.

Definition 3** (Convex concentration property, [Ada15]).**

Let $x$ be a random vector in $\mathbb{R}^{n}$. We will say that $x$ has the convex concentration property with constant $c_{x}$ if for every 1-Lipschitz convex function $g:\mathbb{R}^{n}\to\mathbb{R}$, we have $\mathbb{E}|g(x)|<\infty$ and for every $t>0$,

[TABLE]

If the above tail bound holds for all functions $g(x)=\langle x,u\rangle$ where $u\in\mathbb{R}^{n}$ is any vector with $\|u\|_{2}=1$, then $x$ is called a sub-Gaussian random vector [Ver12]. However, the convex concentration property requires such tail bound to hold for every 1-Lipschitz convex function, and characterizes a subclass for the sub-Gaussian random vectors. See [Ada15, VW15] for examples of such random vectors. As pointed out by [Ada15], $2c_{x}^{2}\geq|\kern-1.07639pt|\kern-1.07639pt|\Sigma_{x}|\kern-1.07639pt|\kern-1.07639pt|_{2}$ always holds.

Improving upon a bound in [VW15], Theorem 2.5 in [Ada15] allows for bounding the deviations of our quadratic forms from their mean, as a function of $|\kern-1.07639pt|\kern-1.07639pt|\Psi_{(i)}^{\prime}\Psi_{(i)}|\kern-1.07639pt|\kern-1.07639pt|_{2}=|\kern-1.07639pt|\kern-1.07639pt|\Psi_{(i)}|\kern-1.07639pt|\kern-1.07639pt|_{2}^{2}$ and $\|\Psi_{(i)}^{\prime}\Psi_{(i)}\|_{F}^{2}$, which is at most $n\|v\|_{0}|\kern-1.07639pt|\kern-1.07639pt|\Psi_{(i)}|\kern-1.07639pt|\kern-1.07639pt|_{2}^{4}$, for $i=1,2$.

These operator norms can be related to certain norms of the block-Toeplitz matrix $\Psi_{n}(B)$. This matrix can in turn be related to a a transfer function that is used in the definitions of $\vartheta_{1}(B)$ and $\vartheta_{2}(B)$. The result is summarized in the next lemma whose proof is given in the Appendix.

Lemma 4**.**

With the above notation, the followings hold

[TABLE]

All in all, we get the following concentration result.

Proposition 5**.**

Consider the autoregressive time series in (1) where all $\epsilon_{t}$, for $t=0,1,\ldots,n-1$, are temporally uncorrelated, have zero mean and variance $\Sigma_{\epsilon}$, and satisfy the convex concentration property with constant $c_{\epsilon}$. Moreover, consider $\{\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{t}\}$ as the partially observed time series corresponding to $\{w_{t}\}$, as characterized in (2). Then, for any fixed vector $v\in\mathbb{R}^{p}$ with $\|v\|_{0}\geq 2\|B_{0}\|_{0}$, any $t>0$, and any $r>0$,

[TABLE]

and

[TABLE]

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

Sketch of Proof of Proposition 5.

The proof is by plugging the bounds of Lemma 4 in Theorem 2.5 of [Ada15] followed by some algebraic manipulations. For the first bound, we bound the operator norm by a scaled Frobenius norm via $\|\Psi_{(1)}^{\prime}\Psi_{(1)}\|_{F}^{2}\leq n|\kern-1.07639pt|\kern-1.07639pt|\Psi_{(1)}|\kern-1.07639pt|\kern-1.07639pt|_{2}^{4}$. For the second bound, we use

[TABLE]

As it is evident from (18) and (19), they can be directly used to bound the quadratic forms in both terms in (16) and in the first term in (17). For the last term in (17), we can derive another concentration result from (18). Observe that

[TABLE]

Remember $\mkern 1.0mu\overline{\mkern-1.0muW\mkern-1.0mu}\mkern 1.0mu=\begin{bmatrix}\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{0}\;\cdots\;\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{n}\end{bmatrix}$ and observe that $\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu$ and $\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu$ are simply subsets of this matrix. Hence, $u^{\prime}\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu$ and $v^{\prime}\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu$ can be expressed similarly through $\Psi_{n+1}(B)$, and choosing certain rows (corresponding to $\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu$ and $\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{Y}\mkern-1.0mu}\mkern 1.0mu$ being subsets of $\mkern 1.0mu\overline{\mkern-1.0muW\mkern-1.0mu}\mkern 1.0mu$) does not increase the operator norm.

Conclusions

This paper presented a new methodology and associated performance guarantees for estimating the parameters of linear vector autoregressive processes by leveraging 1) ideas from sparse regression and the LASSO, 2) estimators designed for robustness to missing data, and 3) concentration results from empirical process theory. Note that optimization problems in (3) and (4) are possibly non-convex because of the $-\|B\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{D}\mkern-1.0mu}\mkern 1.0mu\|_{F}^{2}$ term. Without this term we would have a convex formulation, but would not have a consistent estimator. Our approach generalizes to other measurement schemes beyond multiplication by i.i.d. Bernoulli masks as in (2). In fact, we can adapt our analysis to any covariance-stationary observation process independent of the underlying process whose autocovariance matrices of lag [math] and $1$ have no zero entries.

Appendix A Derivation of the Estimators in (3) and (4)

In this section, we motivate the design of the proposed estimators in (3) and (4). First, we review the relevant notation. For a $p$-dimensional discrete-time, centered, covariance-stationary (wide-sense stationary) process $\{w_{t}\}$, denote the autocovariance function by $\Gamma_{w}(h)=\textup{cov}(w_{t},w_{t+h})$. For a matrix $A$, the transpose is denoted by $A^{\prime}$ and the conjugate transpose is denoted by $A^{\dagger}$. Denote by $\odot$ and $\oslash$ the Hadamard (element-wise) product and division respectively, and by $\otimes$ the Kronecker product. Assuming the process is stationary and ignoring for the moment the fact that $\Sigma_{\epsilon}$ might not be the identity matrix, for any $t$, the best linear estimator for $B_{0}$ is given by

[TABLE]

and replacing the autocovariance with its sample approximation yields

[TABLE]

Given the prior information that $B_{0}$ is sparse, and provided that we have complete information on $\mathcal{X}$ and $\mathcal{Y}$, we can solve either of the following convex optimization problems to estimate $B_{0}$:

[TABLE]

Guarantees on support recovery as well as different error measures for $\widehat{B}_{\textup{full}}$ with respect to $B_{0}$ can be derived through establishing the now well-known lower restricted eigenvalue condition and deviation bound for sample statistics $\mathcal{X}\mathcal{X}^{\prime}$ and $\mathcal{X}\mathcal{Y}^{\prime}$ [BRT09, vdGB09].

Multiplicative Corruption.

When $\mathcal{X}$ and $\mathcal{Y}$ are not fully observed, the above estimators cannot be used anymore. However, going back to (21), we can design a new estimator from scratch if we can estimate the autocovariance matrices $\Gamma_{w}(0)$ and $\Gamma_{w}(1)$ from the given partial data. Suppose that the underlying process $\{w_{t}\}$ is observed through the lens of another covariance-stationary process $\{m_{t}\}$, independent of $\{w_{t}\}$:

[TABLE]

In this case, for any integer value $h$, we have:

[TABLE]

where we used the fact that $\mathbb{E}w_{t}=0$ and the independence of $w_{t}$ and $m_{t}$ implies $\mathbb{E}\mkern 1.0mu\overline{\mkern-1.0muw\mkern-1.0mu}\mkern 1.0mu_{t}=0$ regardless of $m_{t}$ being centered or not. Suppose that for the observation process $\{m_{t}\}$, the autocovariance matrices $\Gamma_{m}(0)$ and $\Gamma_{m}(1)$ have no zero entries. In this case, we have

[TABLE]

which can be plugged in (21) to yield

[TABLE]

whose approximation via the sample autocovariance matrices gives

[TABLE]

Observe that the Hadamard division by $\Gamma_{m}(0)$ can make the quadratic term non-convex.

Missing Data.

A simple scenario for partial observations is when each $m_{t}$ in (22) has entries drawn i.i.d. from a Bernoulli distribution of parameter $1-\delta$, for some $\delta\in[0,1)$. In this case,

[TABLE]

have no zero entries and (23) can be simply expressed as

[TABLE]

where $\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{D}\mkern-1.0mu}\mkern 1.0mu=(\frac{1}{n}{\operatorname{diag}}(\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu\mkern 1.0mu\overline{\mkern-1.0mu\mathcal{X}\mkern-1.0mu}\mkern 1.0mu^{\prime}))^{1/2}\in{\mathbb{R}^{p\times p}}$ is a diagonal matrix of sample autocovariances for each of the $p$ covariates. In this case, with a possibly non-convex quadratic optimization program, we need a constrained optimization program to hope for recovering the target model. Again, given the prior information that $B_{0}$ is sparse, we can use regularization or an $\ell_{1}$-norm constraint. With this consideration, we arrive at the problems in (3) and (4).

Appendix B Estimation Error for Non-convex LASSO

Theorem 6**.**

Consider two random matrices, a symmetric matrix $Q\in\mathbb{R}^{p\times p}$ and a matrix $L\in\mathbb{R}^{p\times p}$, as well as a reference matrix $B_{0}\in\mathbb{R}^{p\times p}$ with $\|B_{0}\|_{0}=k$, and an integer $s\geq 1$. Suppose the following conditions hold:

1. (C1)

For any $v\in\mathbb{R}^{p}$ with $\|v\|_{0}\leq 2s$ and $\|v\|_{2}=1$, there exists $\eta(s)$ such that $|v^{\prime}(Q-\mathbb{E}Q)v|\leq\eta(s)$ with probability at least $1-p_{1}(s)$, 2. (C2)

There exists $\varphi>0$ such that $\|B_{0}Q-L\|_{\infty}\leq\varphi$ with probability at least $1-p_{2}$.

Consider either of the following estimators,

[TABLE]

Consider the largest value of $s$ that satisfies

[TABLE]

while

[TABLE]

and define $\alpha_{\textup{low}}\coloneqq\Lambda_{\min}(\mathbb{E}Q)-27\eta(s)$. Then, for any $B_{0}$ with $\|B_{0}\|_{0}\leq k$, there is a universal positive constant $c_{0}$ such that any global optimum $\widehat{B}$ of (24) with any $b_{0}\geq\|B_{0}\|_{F}$ and $\lambda\geq 2\varphi$ satisfies the bounds

[TABLE]

with probability at least $1-p_{3}(s)-p_{2}$. The same bounds, where $\lambda$ is replaced by $\varphi$, apply to (25). Further, a threshold variant of (24), defined as $\tilde{B}=\{\widehat{B}_{ij}\mathbf{1}_{|\widehat{B}_{ij}|>\lambda}\}_{i,j=1,\ldots,p}\,$, satisfies

[TABLE]

We omit the proof of Theorem 6 for brevity.

Choosing an appropriate $s$ in establishing (C1) might require a lot of algebraic manipulations. Hence, we mention our choice of $s$ in the proof of Theorem 1 using Theorem 6:

[TABLE]

where all notations have been defined in Section 2.

Appendix C Proof of Lemma 4

Using the submultiplicativity of operator norms, we have

[TABLE]

as well as

[TABLE]

in which we used the commutativity of diagonal matrices. Notice that we do not bound the last term with $|\kern-1.07639pt|\kern-1.07639pt|\Psi_{n}(B)|\kern-1.07639pt|\kern-1.07639pt|_{2}$ as a possibly tighter bound is possible. It is easy to see that Lemma 7 implies

[TABLE]

Therefore, it remains to upper bound $|\kern-1.07639pt|\kern-1.07639pt|\Psi_{n}(B)|\kern-1.07639pt|\kern-1.07639pt|_{2}$ and $|\kern-1.07639pt|\kern-1.07639pt|\Psi_{n}(B)|\kern-1.07639pt|\kern-1.07639pt|_{1\to 2}$. However, a closer look reveals that these two quantities are input-output gains, in specific norms, of the following discrete-time linear time-invariant system,

[TABLE]

Since we have assumed $\rho(B)<1$, this system is stable. Moreover, the transfer matrix from $u$ to $x$ is given by

[TABLE]

where $z$ is a complex number. Therefore, we get the right-most set of inequalities in Lemma 4.

Lemma 7**.**

Given a matrix $A$, for any $v$ we have

[TABLE]

where $|\kern-1.07639pt|\kern-1.07639pt|\cdot|\kern-1.07639pt|\kern-1.07639pt|_{1\to 2}$ denotes the largest $\ell_{2}$ norm of columns.

Proof of Lemma 7.

For $v\in\mathbb{R}^{p}$ with $\|v\|_{2}=1$, observe that

[TABLE]

Then, $\|u\odot v\|_{1}=\sum_{i=1}^{p}|u_{i}v_{i}|=|u|^{\prime}|v|\leq\||u|\|_{2}\cdot\||v|\|_{2}=1$ establishes the claim. ∎