An analysis of the SPARSEVA estimate for the finite sample data case

TL;DR

This paper develops a theoretical upper bound for the SPARSEVA estimation error in finite sample scenarios, applicable to strongly convex cost functions and decomposable regularizers, with numerical validation.

Contribution

It introduces a general upper bound for SPARSEVA error applicable to broad conditions and demonstrates its application to sparse regression problems.

Findings

The derived bound effectively estimates the SPARSEVA error.

Numerical results confirm the bound's accuracy and usefulness.

The approach extends SPARSEVA analysis to finite sample data cases.

Abstract

In this paper, we develop an upper bound for the SPARSEVA (SPARSe Estimation based on a VAlidation criterion) estimation error in a general scheme, i.e., when the cost function is strongly convex and the regularized norm is decomposable for a pair of subspaces. We show how this general bound can be applied to a sparse regression problem to obtain an upper bound for the traditional SPARSEVA problem. Numerical results are used to illustrate the effectiveness of the suggested bound.

Equations

$${\hat{\theta }}_{{\lambda }_N}\in \underset{\theta \in {\mathbb{R}}^n}{\text{arg min}}\{\mathcal{L}(\theta ;Z_1^N)+{\lambda }_N\mathcal{R}(\theta )\},$$

$${\hat{\theta }}_{{ϵ}_N}\in \underset{\theta \in {\mathbb{R}}^n}{\text{arg min}}$$

$$\mathcal{L}(\theta ;Z_1^N)\le \mathcal{L}({\hat{\theta }}_{NR};Z_1^N)(1+{ϵ}_N),$$

$${\hat{\theta }}_{NR}\in \underset{\theta \in {\mathbb{R}}^n}{\text{arg min}}\mathcal{L}(\theta ;Z_1^N).$$

$${\overline{\mathcal{M}}}^{\perp }=\{v\in {\mathbb{R}}^n|⟨u,v⟩=0\text{for all}u\in \overline{\mathcal{M}}\},$$

$$\mathcal{R}(\theta +\gamma )=\mathcal{R}(\theta )+\mathcal{R}(\gamma )$$

$$\mathcal{M}(S)=\{\theta \in {\mathbb{R}}^n|{\theta }_j=0\text{for all}j\notin S\}.$$

$${\overline{\mathcal{M}}}^{\perp }(S)=\{\gamma \in {\mathbb{R}}^n|{\gamma }_j=0\text{for all}j\in S\}.$$

$${\mathcal{R}}^{\ast }(v)=\underset{u\in {\mathbb{R}}^n\setminus \{0\}}{\sup }\frac{⟨u,v⟩}{\mathcal{R}(u)}=\underset{\mathcal{R}(u)\le 1}{\sup }⟨u,v⟩,$$

$${▽}^2\mathcal{L}(\theta )⪰mI,$$

$$\mathcal{L}(\theta +\Delta )\ge \mathcal{L}(\theta )+▽\mathcal{L}(\theta {)}^T\Delta +\frac{m}{2}\Vert \Delta {\Vert }_2^2.$$

$$\Psi (\mathcal{M})=\underset{u\in \mathcal{M}\setminus \{0\}}{\sup }\frac{\mathcal{R}(u)}{\Vert u\Vert }.$$

$${\Pi }_{\mathcal{M}}(u)=\underset{v\in \mathcal{M}}{\mathrm{argmin}}\Vert u-v{\Vert }_2.$$

$${\lambda }_{{ϵ}_N}\le 1/{\mathcal{R}}^{\ast }(\nabla \mathcal{L}({\theta }^{\ast })),$$

$$\Vert {\hat{\theta }}_{{ϵ}_N}-{\theta }^{\ast }{\Vert }_2^2$$

$${\lambda }_{{ϵ}_N}>1/{\mathcal{R}}^{\ast }(\nabla \mathcal{L}({\theta }^{\ast })),$$

$$\Vert {\hat{\theta }}_{{ϵ}_N}-{\theta }^{\ast }{\Vert }_2^2$$

$$+\frac{8}{{\kappa }_{\mathcal{L}}^2}\{{\mathcal{R}}^{\ast }(\nabla \mathcal{L}({\theta }^{\ast })){\}}^2{\Psi }^2({\overline{\mathcal{M}}}^{\perp })$$

$$+\frac{4}{{\kappa }_{\mathcal{L}}{\lambda }_{{ϵ}_N}}\mathcal{R}({\theta }_{{\mathcal{M}}^{\perp }}^{\ast }).$$

$$Y_N={\Phi }_N^T{\theta }^{\ast }+e,$$

$${\mathbb{B}}_q(R_q):=\left\{\theta \in {\mathbb{R}}^p|\sum _{i=1}^p|{\theta }_i{|}^q\le R_q\right\},$$

$$\mathcal{L}(\theta )=\frac{1}{2N}\Vert Y_N-{\Phi }_N^T\theta {\Vert }_2^2,$$

$${\hat{\theta }}_{{ϵ}_N}\in \underset{\theta \in {\mathbb{R}}^n}{\text{arg min}}$$

$$\mathcal{L}(\theta )\le \mathcal{L}({\hat{\theta }}_{NR})(1+{ϵ}_N),$$

$${▽}^2\mathcal{L}(\theta )=\frac{1}{N}{\Phi }_N{\Phi }_N^T.$$

$$\exists {\kappa }_{\mathcal{L}}>0\text{s.t.}\frac{1}{N}{\Phi }_N{\Phi }_N^T⪰2{\kappa }_{\mathcal{L}}I.$$

$$\Vert {\hat{\theta }}_{{ϵ}_N}-{\theta }^{\ast }{\Vert }_2^2\le \max (a_1,a_2),$$

$$a_1$$

$$+\frac{\sqrt{32{\sigma }_e^2{s}_{\text{max}}{\chi }_{\beta }^2(N-n)(1+{ϵ}_N)\ln (2/\beta )}}{{\kappa }_{\alpha }N}\Vert {\theta }_{[n_{\eta }+1:n]}^{\ast }{\Vert }_1,$$

$$a_2$$

$$+\frac{\sqrt{32{\sigma }_e^2{s}_{\text{max}}{\chi }_{\beta }^2(N-n)(1+{ϵ}_N)\ln (2/\beta )}}{{\kappa }_{\alpha }N}\Vert {\theta }_{[n_{\eta }+1:n]}^{\ast }{\Vert }_1.$$

$$S_{\eta }^c=\{1,2,...,n\}\setminus S_{\eta };$$

$$\mathcal{M}(S_{\eta })$$

$${\mathcal{M}}^{\perp }(S_{\eta })$$

Full text

††thanks: The material in this paper was not presented at any conference.

, , ,

1 Introduction

Regularization is a well known technique for estimating model parameters from measured input-output data. Its applications are in any fields that are related to constructing mathematical models from observed data, such as system identification, machine learning and econometrics. The idea of the regularization technique is to solve a convex optimization problem constructed from a cost function and a weighted regularizer (regularized M-estimators). There are various types of regularizers that have been suggested so far, such as the $l_{1}$ [19], $l_{2}$ [20] and nuclear norms [5] [6].

During the last few decades, in the system identification community, regularization has been utilised extensively [15], to impose properties of smoothness and sparsity in the estimated models (see, e.g., [13, 22]). Most of this work has focused on analysing the asymptotic properties of an estimator, i.e., when the length of the data goes to infinity. The purpose of this type of analysis is to evaluate the performance of the estimation method to determine if the estimate is acceptable. However, in practice, the data sample size for any estimation problem is always finite, hence, it is difficult to judge the performance of the estimated parameters based on asymptotic properties, especially when the data length is short.

Recently, a number of authors have published research ([1], [4], [12]) aimed at analysing estimation error properties of the regularized M-estimators when the sample size of the data is finite. Specifically, they develop upper bounds on the estimation error for high dimensional problems, i.e., when the number of parameters is comparable to or larger than the sample size of the data. Most of these activities are from the statistics and machine learning communities. Among these works, the paper [12] provides a very elegant and interesting framework for establishing consistency and convergence rates of estimates obtained from a regularized procedure under high dimensional scaling. It determines a general upper bound for regularized M-estimators and then shows how it can be used to derive bounds for some specific scenarios.

Here in this paper we utilize the framework suggested in [12] to develop an upper bound for the estimation error of the M-estimators used in a system identification problem. Here, the M-estimator problems are implemented using the SPARSEVA (SPARSe Estimation based on a VAlidation criterion) framework [16], [17]. The approach in [12] has been developed for penalized estimators, so it has to be suitably modified for SPARSEVA, which is not a penalized estimator, but the solution of a constrained optimization problem. Our aim is to derive an upper bound for the estimation error of the general SPARSEVA estimate. We then apply this bound to a sparse linear regression problem to obtain an upper bound for the traditional SPARSEVA problem with some assumptions on the regression matrix. These assumptions can be considered as the price in order to derive the upper bound. In addition, we also provide numerical simulation results to illustrate the suggested bound of the SPARSEVA estimation error.

The paper is organized as follows. Section 2 formulates the problem. Section 3 provides definitions and properties required for the later analysis. The general bound for the SPARSEVA estimation error is then developed in Section 4. In Section 5, we apply the general bound to the special case when the model is cast in a linear regression framework. Section 6 illustrates the developed bound by numerical simulation. Finally, Section 7 provides conclusions.

1.1 Notation

In this paper, we will use the following notation:

• •

$f(x|0,\sigma^{2})=(2\pi\sigma^{2})^{-1/2}\exp(-x^{2}/2\sigma^{2})$ denotes the probability density function (pdf) of the Normal distribution $\mathcal{N}(0,\sigma^{2})$.

• •

$\chi^{2}_{\beta}(N)$ denotes the value that $P(X<\chi^{2}_{\beta}(N))=1-\beta$, where $X$ is Chi square distributed with $N$ degrees of freedom.

2 Problem Formulation

Let $Z_{1}^{N}=\{Z_{1},...,Z_{N}\}{\in\mathcal{Z}^{N}}$ denote $N$ identically distributed observations with marginal distribution $\mathbb{P}$ in $\mathcal{Z}\subseteq\mathbb{R}^{k}$. $\mathcal{L}:\mathbb{R}^{n}\times\mathcal{Z}^{N}\rightarrow\mathbb{R}$ denotes a convex and differentiable cost function. Let $\theta^{*}\in\text{argmin}_{\theta\in\mathbb{R}^{n}}\overline{\mathcal{L}}(\theta)$ be a minimizer of the population risk $\overline{\mathcal{L}}(\theta)=\mathbb{E}_{Z_{1}^{N}}[\mathcal{L}(\theta;Z_{1}^{N})]$.

The task here is to estimate the unknown parameter $\theta^{*}$ from the data $Z_{1}^{N}$. A well known approach to this problem is to use a regularization technique, i.e., to solve the following convex optimization problem,

[TABLE]

where $\lambda_{N}>0$ is a user-defined regularization parameter and $\mathcal{R}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{+}$ is a norm.

A difficulty when estimating the parameter $\theta^{*}$ using the above regularization technique is that one needs to find the regularization parameter $\lambda_{N}$. The traditional method to choose $\lambda_{N}$ is to use cross validation, i.e., to estimate the parameter $\theta^{*}$ with different values of $\lambda_{N}$, then select the value of $\lambda_{N}$ that provides the best fit to the validation data. This cross validation method is quite time consuming and very dependent on the data. Here we are specifically interested in the SPARSEVA (SPARSe Estimation based on a VAlidation criterion) framework, suggested in [16] and [17], which provides automatic tuning of the regularization parameters. Utilizing the SPARSEVA framework, an estimate of $\theta^{*}$ can be computed using the following convex optimization problem:

[TABLE]

where $\epsilon_{N}>0$ is the regularization parameter and $\hat{\theta}_{NR}$ is the “non-regularized” estimate obtained from minimizing the cost function $\mathcal{L}(\theta;Z_{1}^{N})$, i.e.

[TABLE]

It can be shown [17] that (1) and (2) are equivalent in the sense that there exists a bijection between $\lambda_{N}$ and $\epsilon_{N}$ such that both estimators coincide. However, as discussed in [17, Section V.D], that bijection is data-dependent and it does not seem possible to derive an explicit expression for it. The advantage of the SPARSEVA framework, with respect to (1), is that there are some natural choices of the regularization parameter $\epsilon_{N}$ based the chosen validation criterion. For example, as suggested in [16] [17], $\epsilon_{N}$ can be chosen as $2n/N$ (Akaike Information Criterion (AIC)), $n\text{log}(N)/N$ (Bayesian Information Criterion (BIC)); or as suggested in [8], $n/N$ (Prediction Error Criterion).

For the traditional regularization method described in (1), [12] recently developed an upper bound on the estimation error between the estimate $\hat{\theta}_{\lambda_{N}}$ and the unknown parameter vector $\theta^{*}$. This bound is a function of some constants related to the nature of the data, the regularization parameter $\lambda_{N}$, the cost function $\mathcal{L}$ and the data length $N$. The beauty of this bound is that it quantifies the relationship between the estimation error and the finite data length $N$. Through this relationship, it is easy to confirm most of the properties of the estimate $\hat{\theta}_{\lambda_{N}}$ in the asymptotic scenario, i.e. $N\rightarrow\infty$, which were developed in the literature some time ago ([9], [11]).

Inspired by [12], our goal is to derive a similar bound for the SPARSEVA estimate $\hat{\theta}_{\epsilon_{N}}$, i.e, we want to know how much the SPARSEVA estimate $\hat{\theta}_{\epsilon_{N}}$ differs from the true parameter $\theta^{*}$ when the data sample size $N$ is finite. Note that the notation and techniques used in this paper are similar to [12]; however, in [12], the convex optimization problem is posed in the traditional regularization framework (1), while in this paper, the optimization problem is based on the SPARSEVA regularization (2).

3 Definitions and Properties of the Norm $\mathcal{R}(\theta)$ and the Cost Function $\mathcal{L}(\theta)$

In this section, we provide descriptions of some definitions and properties of the norm $\mathcal{R}(\theta)$ and the cost function $\mathcal{L}(\theta;Z_{1}^{N})$, needed to establish an upper bound on the estimation error. Note that we only provide a brief summary such that the research described in this paper can be understood. Readers can find a more detailed discussion in [12].

3.1 Decomposability of a Norm

Let us consider a pair of arbitrary linear subspaces of $\mathbb{R}^{n}$, $(\mathcal{M},\overline{\mathcal{M}})$, such that $\mathcal{M}\subseteq\overline{\mathcal{M}}$. The orthogonal complement of the space $\overline{\mathcal{M}}$ is then defined as,

[TABLE]

where $\langle\cdot,\cdot\rangle$ is the inner product that maps $\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}$.

The norm $\mathcal{R}$ is said to be decomposable with respect to $(\mathcal{M},\overline{\mathcal{M}}^{\perp})$ if

[TABLE]

for all $\theta\in\mathcal{M}$ and $\gamma\in\overline{\mathcal{M}}^{\perp}$.

There are many combinations of norms and vector spaces that satisfy this property (cf. [12]). An example is the $l_{1}$ norm and the sparse vector space defined (5). For any subset $S\subseteq\{1,2,\dots,n\}$ with cardinality $s$, define the model subspace $\mathcal{M}$ as,

[TABLE]

Now if we define $\overline{\mathcal{M}}(S)=\mathcal{M}(S)$, then the orthogonal complement $\overline{\mathcal{M}}(S)$, with respect to the Euclidean inner product, can be computed as follows,

[TABLE]

As shown in [12], the $l_{1}$-norm is decomposable with respect to the pair $(\mathcal{M}(S),\mathcal{M}^{\perp}(S))$.

3.2 Dual Norm

For a given inner product $\langle\cdot,\cdot\rangle$, the dual of the norm $\mathcal{R}$ is defined by,

[TABLE]

where sup is the supremum operator.

Based on the above definition, one can easily see that the dual of the $l_{1}$ norm, with respect to the Euclidean inner product, is the $l_{\infty}$ norm [12].

3.3 Strong Convexity

A twice differentiable function $\mathcal{L}(\theta):\mathbb{R}^{n}\rightarrow\mathbb{R}$ is strongly convex on $\mathbb{R}^{n}$ when there exists an $m>0$ such that its Hessian $\bigtriangledown^{2}\mathcal{L}(\theta)$ satisfies,

[TABLE]

for all $\theta\in\mathbb{R}^{n}$ [3]. This is equivalent to the statement that the minimum eigenvalue of $\bigtriangledown^{2}\mathcal{L}(\theta)$ is not smaller than $m$ for all $\theta\in\mathbb{R}^{n}$.

An interesting consequence of the strong convexity property in (7) is that for all $\theta,\Delta\in\mathbb{R}^{n}$, we have,

[TABLE]

The inequality in (8) has a geometric interpretation in that the graph of the function $\mathcal{L}(\theta)$ has a positive curvature at any $\theta\in\mathbb{R}^{n}$. The term $m/2$ for the largest $m$ satisfying (7) is typically known as the curvature of $\mathcal{L}(\theta)$.

3.4 Subspace Compatibility Constant

For a given norm $\mathcal{R}$ and an error norm $\|\cdot\|$, the subspace compatibility constant of a subspace $\mathcal{M}\subseteq\mathbb{R}^{n}$ with respect to the pair $(\mathcal{R},\|\cdot\|)$ is defined as,

[TABLE]

This quantity measures how well the norm $\mathcal{R}$ is compatible with the error norm $\|.\|$ over the subspace $\mathcal{M}$. As shown in [12], when $\mathcal{M}$ is $\mathbb{R}^{s}$, the regularized norm $\mathcal{R}$ is the $l_{1}$ norm, and the error norm is the $l_{2}$ norm, then the subspace compatibility constant is $\Psi(\mathcal{M})=\sqrt{s}$. Notice also that $\Psi(\mathcal{M})$ is finite, due to the equivalence of finite dimensional norms.

3.5 Projection Operator

The projection of a vector $u$ onto a space $\mathcal{M}$, with respect to the Euclidean norm, is defined by the following,

[TABLE]

In the sequel, to simplify the notation, we will write $u_{\mathcal{M}}$ to denote $\Pi_{\mathcal{M}}(u)$.

4 Analysis of the Regularization Technique using the SPARSEVA

In this section, we apply the properties described in Section 3 to derive an upper bound on the error between the SPARSEVA estimate $\hat{\theta}_{\epsilon_{N}}$ and the unknown parameter $\theta^{*}$. This upper bound is described in the following theorem.

Theorem 4.1

Assume $\mathcal{R}$ is a norm and is decomposable with respect to the subspace pair ($\mathcal{M},\overline{\mathcal{M}}^{\perp}$) and the cost function $\mathcal{L}(\theta)$ is differentiable and strongly convex with curvature $\kappa_{\mathcal{L}}$. Consider the SPARSEVA problem in (2), then the following properties hold:

1. i.

*When $\epsilon_{N}>0$, there exists a Lagrange multiplier, $\lambda_{N}=\lambda_{\epsilon_{N}}$, such that (1) and (2) have the same solution. * 2. ii.

*Any optimal solution $\hat{\theta}_{\epsilon_{N}}\neq 0$ of the SPARSEVA problem (2) satisfies the following inequalities: *

• •

If $\epsilon_{N}$ is chosen such that

[TABLE]

then

[TABLE]

• •

If $\epsilon_{N}$ is chosen such that

[TABLE]

then

[TABLE]

Proof. See the Appendix (Section A.2). $\square$

Remark 4.1

Note that Theorem 4.1 is intended to provide an upper bound on the estimation error for the general SPARSEVA problem (2). At this stage, it is hard to evaluate, or quantify, the value on the right hand side of the inequalities (11) and (12) as they still contain the term $\lambda_{\epsilon_{N}}$ and other abstract terms. However, in the later sections of this paper, from this general upper bound, we will provide bounds on the estimation errors for some specific scenarios.

Remark 4.2

The bound in Theorem 4.1 is actually a family of bounds. For each choice of the pair of subspaces ($\mathcal{M},\overline{\mathcal{M}}^{\perp}$), there is one bound for the estimation error. Hence, in the usual sense, to apply Theorem 4.1 for any specific scenario, the goal is to choose $\mathcal{M}$ and $\overline{\mathcal{M}}^{\perp}$ to obtain an optimal rate of the bound.

5 An Upper Bound for Sparse Regression

In this section, we illustrate how to apply Theorem 4.1 to derive an upper bound of the error between the SPARSEVA estimate $\hat{\theta}_{\epsilon_{N}}$ and the true parameter $\theta^{*}$ for the following linear regression model,

[TABLE]

where $\theta^{*}\in\mathbb{R}^{n}$ is the unknown parameter that is required to be estimated; $e\in\mathbb{R}^{N}$ is the disturbance noise; $\Phi_{N}\in\mathbb{R}^{n\times N}$ is the regression matrix and $Y_{N}\in\mathbb{R}^{N}$ is the output vector. Here, we make the following assumption on the true parameter $\theta^{*}$,

Assumption 5.1

The true parameter $\theta^{*}$ is “weakly” sparse, i.e. $\theta^{*}\in\mathbb{B}_{q}(R_{q})$, where,

[TABLE]

with $q\in[0,1]$ being a constant.

Using the SPARSEVA framework in (2) with $\mathcal{R}$ chosen as the $l_{1}$ norm and the cost function $\mathcal{L}(\theta)$ chosen as,

[TABLE]

then an estimate of $\theta^{*}$ in (13) can be found by solving the following problem,

[TABLE]

with $\hat{\theta}_{NR}=(\Phi_{N}\Phi_{N}^{T})^{-1}\Phi_{N}Y_{N}$ and $\epsilon_{N}>0$ being the user-defined regularization parameter. Now $\epsilon_{N}$ can be chosen as either $2n/N$ or $\text{log}(N)n/N$ as suggested in [16]; or $n/N$ as suggested in [8].

Remark 5.1

Note that the sparse regression problem is very common in system identification and is often used to obtain a low order linear model by regularization.

Remark 5.2

For Assumption 5.1, note that when $q=0$, under the convention that $0^{0}=0$, the set in (14) corresponds to an exact sparsity set, where all the elements belonging to the set have at most $R_{0}$ non-zero entries. Generally, for $q\in(0,1]$, the set $\mathbb{B}_{q}(R_{q})$ forces the ordered absolute values of $\theta^{*}$ to decay with a certain rate.

5.1 An Analysis on the Strong Convexity Property and the Curvature of the $l_{2}$ norm Cost Function

Consider the convex optimization problem in (16), the Hessian matrix of the cost function $\mathcal{L}(\theta)$ is computed as,

[TABLE]

To prove that $\mathcal{L}(\theta)$ is strongly convex, we need to prove,

[TABLE]

We see that the requirement in (17) coincides with the requirement of persistent excitation of the input signal in a system identification problem. If an experiment is well-designed, then the input signal $u(t)$ needs to be persistently exciting of order $n$, i.e., the matrix $\Phi_{N}\Phi_{N}^{T}$ is a positive definite matrix. This means that the condition in (17) is always satisfied for any linear regression problem derived from a well posed system identification problem. This means that for any choice of the regression matrix $\Phi_{N}$ that satisfies the persistent excitation condition, there exists a positive curvature $\kappa_{\mathcal{L}}$ of the cost function $\mathcal{L}(\theta)$.

Consider $\Phi_{N}\in\mathbb{R}^{n\times N}$ to be a matrix where each row $\Phi_{N,j}$ is sampled from a Normal distribution of zero mean and covariance matrix $\Sigma\in\mathbb{R}^{N\times N}$, i.e., $\Phi_{N,j}\sim\mathcal{N}(0,\Sigma),\ \ \ \forall j=1,..,n.$ We then denote the distribution of the smallest eigenvalue of $N^{-1}\Phi_{N}\Phi_{N}^{T}$ to be $P(x|\Sigma,N,n)$, means that given a probability $1-\alpha\ ,\ 0\leq\alpha\leq 1$, there exists a value $w_{\text{min}}$ such that $N^{-1}\Phi_{N}\Phi_{N}^{T}\succeq w_{\text{min}}I$, for any matrix $\Phi_{N}$ constructed following the above assumption. Then the global curvature $\kappa$, i.e. the curvature that satisfies (17) for any regression matrix $\Phi_{N}$, can be expressed as $(1/2)w_{\text{min}}$. For the rest of the paper, we will denote by $\kappa_{\alpha}$ lower bound on the global curvature $\kappa$ with probability $1-\alpha\ ,\ 0\leq\alpha\leq 1$.

5.2 Assumptions

For the linear regression in (13), the following assumptions are made:

Assumption 5.2

The rows $\Phi_{N,j},j=1,...n$ of the regressor matrix $\Phi_{N}$ are distributed as $\Phi_{N,j}\sim\mathcal{N}(0,\Sigma)$, where $\Sigma\in\mathbb{R}^{N\times N}$ is a constant, symmetric, positive definite matrix.

Note that an obvious practical case where Assumption 5.2 is satisfied is when the model is FIR and the input signal being white noise or coloured noise.

Assumption 5.3

The noise vector $\text{e}\in\mathbb{R}^{N}$ is Gaussian with i.i.d. $\mathcal{N}(0,\sigma_{e}^{2})$ entries.111The assumption of Gaussian noise is fairly standard in system identification. However, this assumption can be relaxed to ‘sub-Gaussian’ noise (i.e., when the tails of the noise distribution decay like $e^{-\alpha x^{2}}$) at the expense of longer derivations.

5.3 Developing the Upper Bound

The following theorem provides an upper bound on the estimation error $\|\hat{\theta}_{\epsilon_{N}}-\theta^{*}\|_{2}$ for the optimization problem in (16) in the case of weakly sparse estimates.

Theorem 5.1

*Suppose Assumptions 5.2, 5.3 and 5.1 hold, when N is large, then with probability $(1-\alpha)(1-4n\beta)$ ($0\leq\alpha\leq 1,\ 0\leq\beta\leq 1$), if $\hat{\theta}_{\epsilon_{N}}\neq 0$ we have the following inequality

[TABLE]

where

[TABLE]

where $\kappa_{\alpha}$ is a lower bound on the curvature of the regression matrix (i.e., half the smallest eigenvalue of $N^{-1}\Phi_{N}^{T}\Phi_{N}$) with probability $1-\alpha$, $n_{\eta}$ is any integer between $1$ and $n$, $\theta_{[n_{\eta}+1:n]}^{*}$ is the vector formed from the $n-n_{\eta}$ smallest (in magnitude) entries of $\theta^{*}$, and $s_{\text{max}}$ is the maximum singular value of the matrix $\Sigma$.

Proof. This proof relies on three preliminary results introduced in Appendix A.3. For an integer $n_{\eta}\in\{1,\dots,n\}$, define $S_{\eta}$ as the set of the indices of the $n_{\eta}$ largest (in magnitude) entries of $\theta^{*}$, and its complementary set $S^{c}_{\eta}$ as

[TABLE]

with the corresponding subspaces $\mathcal{M}(S_{\eta})$ and $\mathcal{M}^{\perp}(S_{\eta})$ as,

[TABLE]

Using the definition of the subspace compatibility constant described in Section 3, we have,

[TABLE]

where $|S|$ denotes the cardinality of $S$.

Now, for Theorem 4.1 to generate an upper bound for the problem (16), we need to establish an upper bound on $\|\theta^{*}_{\mathcal{M}^{\perp}(S_{\eta})}\|_{1}$. Based on the definition of the subspace $\mathcal{M}^{\perp}(S_{\eta})$, we have,

[TABLE]

where $\theta_{[n_{\eta}+1:n]}^{*}$ denotes the vector formed from the $n-n_{\eta}$ smallest (in magnitude) entries of $\theta^{*}$. Define $\kappa_{\alpha}$ as a lower bound on the global curvature of the regression matrix, i.e. half the smallest eigenvalue of $\Phi_{N}^{T}\Phi_{N}$, with probability $1-\alpha,\ 0\leq\alpha\leq 1$. Substituting the results of Propositions A.1-A.3 from Appendix A.3, (21) and (22) into the bound in Theorem 4.1, then with $n_{\eta}$ being any integer between $1$ and $n$, we have the following bounds:

• •

If $\epsilon_{N}$ is chosen such that

[TABLE]

then, with probability at least $(1-\alpha)(1-2n\beta)$,

[TABLE]

• •

If $\epsilon_{N}$ is chosen such that

[TABLE]

then, with probability at least $(1-\alpha)(1-4n\beta)$,

[TABLE]

Therefore, for $n_{\eta}$ being any integer between $1$ and $n$, with probability at least ($1-\alpha$)($1-4n\beta)$, we have

[TABLE]

where

[TABLE]

[TABLE]

$\square$

Remark 5.3

The bound in Theorem 5.1 is also a family of bounds, one for each value of $n_{\eta}$.

Remark 5.4

When $\Sigma=\sigma_{u}I$, i.e. the model is FIR and the input $u(t)$ is white noise, then $s_{\text{max}}=\sigma_{u}$ and the generalized Chi square distribution $\chi^{2}(\Sigma,I)$ becomes the Chi square distribution $\sigma_{u}\chi^{2}(N)$.

Remark 5.5

Note that the developed bound in Theorem 5.1 depends on the true parameter $\theta^{*}$, which is unknown but constant. Using a similar proof as in Proposition 2.3 of [7], we can derive under Assumption 5.1 an upper bound for the term $\|\theta_{[n_{\eta}+1:n]}^{*}\|_{1}$. Specifically, we have,

[TABLE]

Since $S_{\eta}$ is the set of the indices of the $n_{\eta}$ largest (in magnitude) entries of $\theta^{*}$, i.e. $|\theta^{*}_{[i]}|\leq|\theta^{*}_{n_{\eta}}|,\ \forall i=n_{\eta}+1,...,n$, hence,

[TABLE]

Using the same argument, we have,

[TABLE]

Therefore,

[TABLE]

This means we can always place an upper bound on the term $\|\theta_{[n_{\eta}+1:n]}^{*}\|_{1}$ by a known constant which depends on the nature of the true parameter $\theta^{*}$. Therefore, from Theorem 5.1, we can see that the estimation error $\|\hat{\theta}_{\epsilon_{N}}-\theta^{*}\|^{2}_{2}=O_{p}(N^{-1/2})$ [17]. This confirms the result in [17], that in the asymptotic case, when $\epsilon_{N}>0$, the SPARSEVA estimate $\hat{\theta}_{\epsilon_{N}}$ converges to the true parameter $\theta^{*}$.

6 Numerical Evaluation

In this section, numerical examples are presented to illustrate the bound $\|\hat{\theta}_{\epsilon_{N}}-\theta^{*}\|^{2}_{2}$ as stated in Theorem 5.1. In Section 6.1, we consider the case when the input is Gaussian white noise whilst in Section 6.2, the input is a correlated signal with zero mean.

6.1 Gaussian White Noise Input

In this section, a random discrete time system with a random model order between 1 and 10 is generated using the command drss from Matlab. The system has poles with magnitude less than 0.9. Gaussian white noise is added to the system output to give different levels of SNR, e.g. 30dB, 20dB and 10dB. For each noise level, 50 different input excitation signals (Gaussian white noise with variance 1) and output noise realizations are generated. For each set of input and output data, the system parameters are estimated using a different sample size, i.e., $N=[450,1000,5000,10000,50000,100000]$.

The FIR model structure is used here in order to construct the SPARSEVA problem (16). The number of parameters $n$ of the FIR model is set to be 35. The regularization parameter $\epsilon_{N}$ is chosen as $n/N$ [8].

We then compute the upper bound of $\|\hat{\theta}_{\epsilon_{N}}-\theta^{*}\|_{2}$ using (18) with different values of $n_{\eta}$, i.e. $n_{\eta}=[10,15,25]$. The probability parameters $\alpha$ and $\beta$ are chosen to be $0.02$ and $0.001$ respectively. Related to the computation of the universal constant $\kappa$ corresponding to the distribution $\mathcal{N}(0,\Sigma)$, note that, in reality, it is very difficult to compute its exact distribution $P(x|\Sigma,N,n)$, hence, here we use an empirical method to compute the distribution $P(x|\Sigma,N,n)$. The idea is to generate a large number of random matrices $\Phi_{N}$, compute the smallest eigenvalue of $N^{-1}\Phi_{N}\Phi_{N}^{T}$, and then build a histogram of these values, which is an approximation of $P(x|\Sigma,N,n)$. Then we compute the value of $w_{\text{min}}$ to ensure the inequality $N^{-1}\Phi_{N}\Phi_{N}^{T}\succeq w_{\text{min}}I$ occurs with probability $1-\alpha$. Finally, $\kappa_{\alpha}$ is computed using the formula $\kappa_{\alpha}=(1/2)w_{\text{min}}$.

With the setting described above, the probability of the upper bound being correct is $(1-\alpha)(1-4n\beta)=0.84$. This upper bound will be compared with $\|\hat{\theta}_{\epsilon_{N}}-\theta^{*}\|_{2}$. Note that we plot both the upper bound and the true estimation errors on a logarithmic scale.

Plots of the estimation error versus the data length $N$ with different noise levels are displayed in Figures 1 to 3. In Figures 1 to 3, the red lines are the true estimation errors from 50 estimates using the SPARSEVA framework. The magenta, blue and cyan lines are the upper bounds developed in Theorem 5.1, which correspond to $n_{\eta}=[10,15,25]$, respectively. We can see that the plots confirm the bound developed in Theorem 5.1 for all noise levels. When $N$ becomes large, the estimation error and the corresponding upper bound become smaller. When $N$ goes to infinity, the estimation error will tend to 0. Note that the bounds are slightly different for the chosen values of $n_{\eta}$, however, not significantly. As can be seen, the bounds are relatively insensitive to the choice of $n_{\eta}$.

In addition, we plot another graph, shown in Fig. 4, to compare the proposed upper bound and the true estimation errors corresponding to different value of $\epsilon_{N}$, i.e. $\dfrac{n}{N}$ (PEC), $\dfrac{2n}{N}$ (AIC) and $\dfrac{log(N)n}{N}$ (BIC). The blue lines are the upper bounds developed in Theorem 5.1, which correspond to $n_{\eta}=25$, with the three different values of $\epsilon_{N}$. The magenta (BIC), green (AIC) and red (PEC) lines are the true estimation errors from 50 estimates (for each value of $\epsilon_{N}$) using the SPARSEVA framework. We can see that the plot again confirms the validity of the proposed upper bound for all choices of $\epsilon_{N}$. Note that the upper bound is not extremely tight, it is quite conservative, however, it is the price to usually pay for finite sample bounds with a general SPARSEVA setting, i.e. the regularized parameter $\epsilon_{N}$ can be any positive value. When $\epsilon_{N}$ is larger, the upper bound will be closer to the true estimate error.

6.2 Coloured Noise Input

In this section, a random discrete time system with a random model order between 1 and 10 is generated using the command drss from Matlab. The system has poles with magnitude less than 0.9. White noise is added to the system output with different levels of SNR, e.g. 30dB, 20dB and 10dB. For each noise level, 50 different input excitation signals and output noise realizations are generated. For each set of input and output data, the system parameters are estimated using different sample sizes, i.e. $N=[450,1000,5000,10000,50000]$.

Here, the input signal is generated by filtering a zero mean Gaussian white noise with unit variance through the filter,

[TABLE]

Due to this filtering, the covariance matrix of the regression matrix distribution will not be of a diagonal form. Note that this is a completely different scenario to that in Section 6.1.

The FIR model structure is used here in order to construct the linear regression for the SPARSEVA problem (16). The number of parameters $n$ of the FIR model is set to be 35. The regularization parameter, $\epsilon_{N}$, is chosen as $n/N$ [8].

We then compute the upper bound of $\|\hat{\theta}_{\epsilon_{N}}-\theta^{*}\|_{2}$ using (18) with different values of $n_{\eta}$, i.e. $n_{\eta}=[10,15,25]$. The probability parameters $\alpha$ and $\beta$ are chosen to be $0.02$ and $0.001$ respectively. With this setting, the probability of the upper bound being correct is $(1-\alpha)(1-4n\beta)=0.84$. This upper bound will be compared with $\|\hat{\theta}_{\epsilon_{N}}-\theta^{*}\|_{2}$.

Plots of the upper bound as stated in Theorem 5.1 and the true estimation error $\|\hat{\theta}_{\epsilon_{N}}-\theta^{*}\|_{2}$ are displayed in Figures 5 to 7. In Figures 5 to 7, the red lines are the true estimation errors from 50 estimates using the SPARSEVA framework. The magenta, blue and cyan lines are the upper bounds developed in Theorem 5.1, which correspond to $n_{\eta}=[10,15,25]$ respectively. We can see that the plots confirmed the bound developed in Theorem 5.1 for all noise levels. When $N$ becomes large, the estimation error and the corresponding upper bound become smaller. When $N$ goes to infinity, the estimation error will tend to 0. Note that the bounds are slightly different for the chosen values of $n_{\eta}$, however, not significantly. As can be seen, the bounds are relatively insensitive to the choice of $n_{\eta}$.

7 Conclusion

The paper provides an upper bound on the SPARSEVA estimation error in the general case, for any choice of strongly convex cost function and decomposable norm. We also evaluate the bound for a specific scenario, i.e., a sparse regression estimate problem. Numerical results confirm the validity of the developed bound for different input signals with different output noise levels for different choices of the regularization parameters.

Appendix A Appendix

A.1 Background knowledge

First, we cite a lemma directly from [12], to enable the proof of Theorem 4.1 to be constructed.

Lemma A.1

For any norm $\mathcal{R}$ that is decomposable with respect to $(\mathcal{M,\overline{\mathcal{M}}^{\perp}})$; and any vectors $\theta$, $\Delta$, we have

[TABLE]

Recall that $\Delta_{\overline{\mathcal{M}}^{\perp}}$ is the Euclidean projection of $\Delta$ onto $\overline{\mathcal{M}}^{\perp}$(see Section 3.5), and similarly for the other terms in (24).

Proof. See the supplementary material of [12]. $\square$

We now quote the following lemma from [2] with modification to fit with the notation used in the SPARSEVA problem (2). This lemma helps us to find some important properties related to the SPARSEVA estimate $\hat{\theta}_{\lambda_{N}}$. Based on these properties, in the next section we can derive the upper bound on the estimation error. Note that for notational simplicity, we will denote $\mathcal{L}(\theta;Z_{1}^{N})$ as $\mathcal{L}(\theta)$.

Lemma A.2

Consider the convex optimization problem in (2). Then the pair ($\hat{\theta}_{\epsilon_{N}},\lambda_{\epsilon_{N}}$), with $\hat{\theta}_{\epsilon_{N}}\neq 0$, has the property that $\hat{\theta}_{\epsilon_{N}}$ is the solution of the problem (2) and $\lambda_{\epsilon_{N}}$ is the Lagrange multiplier if and only if all of the following hold:

$\lambda_{\epsilon_{N}}\in R^{+}$; 2. 2.

the function $\mathcal{R}(\theta)+\lambda\{\mathcal{L}(\theta)-\mathcal{L}(\hat{\theta}_{NR})(1+\epsilon_{N})\}$ attains its minimum over $\mathbb{R}^{n}$ at $\hat{\theta}_{\epsilon_{N}}$; and 3. 3.

$\mathcal{L}(\hat{\theta}_{\epsilon_{N}})-\mathcal{L}(\hat{\theta}_{NR})(1+\epsilon_{N})=0$.

Proof. See the proof of Theorem 2.2 in [2]. The third condition in the cited theorem is a complementary slackness condition, which reduces to condition (3) here if $\hat{\theta}_{\epsilon_{N}}\neq 0$ (cf. [17, Lemma II.2]). $\square$

A.2 Proof of Theorem 4.1

First we need to prove that there exists a Lagrange multiplier for the SPARSEVA problem (2). We can assume without loss of generality that $\mathcal{L}(\hat{\theta}_{NR})\neq 0$, since otherwise we can take $\lambda_{\epsilon_{N}}=0$. According to [2], the Lagrange multiplier for a convex optimization problem with constraint exists when the Slater condition is satisfied. Specifically, for the SPARSEVA problem (2), the Lagrange multiplier $\lambda_{\epsilon_{N}}$ exists when there exists a $\theta_{1}$ such that $\mathcal{L}(\theta_{1})<\mathcal{L}(\hat{\theta}_{NR})(1+\epsilon_{N})$. If $\epsilon_{N}>0$, and $\mathcal{L}(\hat{\theta}_{NR})\neq 0$, there always exists a parameter vector $\theta_{1}$ such that $\mathcal{L}(\theta_{1})<\mathcal{L}(\hat{\theta}_{NR})(1+\epsilon_{N})$ (just take $\theta_{1}=\hat{\theta}_{NR}$). Therefore, there exists a Lagrange multiplier $\lambda_{\epsilon_{N}}$ for the SPARSEVA problem.

Now that we have confirmed the existence of the Lagrange multiplier $\lambda_{\epsilon_{N}}$, consider the function $\mathcal{F}(\Delta)$ defined as follows,

[TABLE]

Using the strong convexity condition of $\mathcal{L}(\theta)$,

[TABLE]

From (6), we have that

[TABLE]

Next, combining the inequality (27) and the triangle inequality, i.e., $\mathcal{R}(\Delta)\leq\mathcal{R}(\Delta_{\overline{\mathcal{M}}})+\mathcal{R}(\Delta_{\overline{\mathcal{M}}^{\perp}})$, we have,

[TABLE]

therefore,

[TABLE]

Now, combining (25), (26), (28), and Lemma A.1,

[TABLE]

Notice that, when $\hat{\theta}_{\epsilon_{N}}$ is the estimate of the SPARSEVA problem, then property 2 in Lemma A.2 states that the function $\mathcal{R}(\theta)+\lambda\{\mathcal{L}(\theta)-\mathcal{L}(\hat{\theta}_{NR})(1+\epsilon_{N})\}$ attains its minimum over $\mathbb{R}^{n}$ at $\hat{\theta}_{\epsilon_{N}}$, which means,

[TABLE]

Hence,

[TABLE]

or, taking $\theta=\theta^{*}$ and defining $\hat{\Delta}_{\epsilon_{N}}:=\hat{\theta}_{\epsilon_{N}}-\theta^{*}$,

[TABLE]

Combining (A.2) with (30), we then have,

[TABLE]

Now we consider two cases.

Case 1**:** $\ \lambda_{\epsilon_{N}}\leq 1/\mathcal{R}^{*}(\nabla\mathcal{L}(\theta^{*}))$

From (31), we have,

[TABLE]

By the definition of subspace compatibility,

[TABLE]

Now we also have that

[TABLE]

Therefore,

[TABLE]

Substituting this into (32) gives,

[TABLE]

Note that for a quadratic polynomial $f(x)=ax^{2}+bx+c$, with $a>0$, if there exists $x\in\mathbb{R}^{+}$ that makes $f(x)\leq 0$, then such $x$ must satisfy

[TABLE]

Since $(A+B)^{2}\leq 2A^{2}+2B^{2}$ for all $A,B\in\mathbb{R}$,

[TABLE]

Applying this inequality to (34), we have,

[TABLE]

Case 2**:** $\lambda_{\epsilon_{N}}>1/\mathcal{R}^{*}(\nabla\mathcal{L}(\theta^{*}))$

Using a similar analysis as in Case 1,

[TABLE]

Substituting (37) and (33) into (31), we obtain,

[TABLE]

Now using the inequality (35), yields,

[TABLE]

Applying the inequality $(A+B)^{2}\leq 2A^{2}+2B^{2}$ to the first term in (39) gives,

[TABLE]

Note that $0<1/\lambda_{\epsilon_{N}}<\mathcal{R}^{*}(\nabla\mathcal{L}(\theta^{*}))$, therefore,

[TABLE]

We also have,

[TABLE]

Therefore, combining (40) and (41),

[TABLE]

$\square$

A.3 Preliminary propositions for Theorem 5.1

In this Appendix we present three propositions that assist in the development of the proof of Theorem 5.1:

Proposition A.1

Consider the optimization problem in (16), and denote by $\lambda_{\epsilon_{N}}$ the corresponding Lagrange multiplier of its constraint. Then, if $\hat{\theta}_{\epsilon_{N}}\neq 0$, $\lambda_{\epsilon_{N}}$ can be computed as

[TABLE]

Proof. See Appendix A.4. $\square$

Proposition A.2

Suppose Assumptions 5.2 and 5.3 hold. Then, with probability $1-n\beta$ ($0\leq\beta\leq{1/n}$), we have

[TABLE]

where $s_{\text{max}}$ is the maximum element on the diagonal of the matrix $\Sigma$. In particular, choosing a specific value for $t$,

[TABLE]

with probability at least $1-2n\beta$ ($0\leq\beta\leq 1/2n$).

Proof. See Appendix A.5. $\square$

Proposition A.3

Suppose Assumptions 5.2 and 5.3 hold, then with probability at least $1-n\beta$ ($0\leq\beta\leq{1/n}$), we have

[TABLE]

where $s_{\text{max}}$ is the maximum element on the diagonal of the matrix $\Sigma$. In particular, choosing a specific value for $t$,

[TABLE]

with probability at least $1-2n\beta$ ($0\leq\beta\leq 1/2n$).

Proof. See Appendix A.6. $\square$

A.4 Proof of Proposition A.1

Let us rewrite the SPARSEVA problem (16),

[TABLE]

in the Lagrangian form (1) using Lemma A.2. The Lagrangian of the optimization problem (45) is,

[TABLE]

The subdifferential of $g(\theta,\lambda)$ can be computed as

[TABLE]

where $\text{v}=(v_{1},\dots,v_{m})^{T}$ is of the form

[TABLE]

Using property 2 of Lemma A.2, when $\hat{\theta}_{\epsilon_{N}}$ is a solution of the SPARSEVA problem (16) and $\lambda_{\epsilon_{N}}$ is a Lagrange multiplier, we have,

[TABLE]

for some v of the form in (48). Note that when ${\hat{\theta}_{\epsilon_{N}}}\neq\textbf{0}$, $\|\text{v}\|_{\infty}=1$, which means that

[TABLE]

Since $\nabla\mathcal{L}(\hat{\theta}_{\epsilon_{N}})=\dfrac{1}{N}\Phi_{N}(Y_{N}-\Phi_{N}^{T}\hat{\theta}_{\epsilon_{N}})$, we can also write $\lambda_{\epsilon_{N}}$ as

[TABLE]

$\square$

Note that this proof is similar to the one in [14], where an expression was derived for the Lagrange multiplier in the traditional $l_{1}$ norm regularization problem (the LASSO). Here we have derived the Lagrange multiplier for the SPARSEVA problem as given in (16).

A.5 Proof of Proposition A.2

For the linear regression (13) and the choice of $\mathcal{L}(\theta)$ in (15),

[TABLE]

Denote $R_{j}$ as the $j^{th}$ row of the matrix $\Phi_{N}$, then $\nabla\mathcal{L}(\theta^{*})$ can be computed as,

[TABLE]

consider the variable $Z=N^{-1}R_{j}e$, using Assumption 5.3 on the disturbance noise $e$, $e\sim\mathcal{N}(0,\sigma_{e}^{2})$, we have,

[TABLE]

Now in order to derive a bound for $\nabla\mathcal{L}(\theta^{*})$, we first derive an upper bound for the variance $N^{-2}\sigma_{e}^{2}R_{j}R_{j}^{T}$ of the distribution in (50). Since $R_{j}\sim\mathcal{N}(0,\Sigma)$, we have,

[TABLE]

where $\chi^{2}(\Sigma,I)$ is the generalized Chi squared with parameters $\Sigma$ and $I$. Hence, with probability $1-\beta,\ 0\leq\beta\leq 1$, we have,

[TABLE]

Hence, the variance of the distribution of the variable $N^{-1}R_{j}e$, is,

[TABLE]

with probability $1-\beta$.

Note that from (50), for any $t>0$, we have,

[TABLE]

where $f(x|0,N^{-2}\sigma_{e}^{2}R_{j}R_{j}^{T})$ denotes the pdf of the Normal distribution $\mathcal{N}(0,N^{-2}\sigma_{e}^{2}R_{j}R_{j}^{T})$. This gives,

[TABLE]

This expression can be bounded from below using the standard result that $P(|\mathcal{N}(0,\sigma^{2})|>t)\leq 2\exp(-t^{2}/2\sigma^{2})$ [21, Eq. (5.5)], to obtain

[TABLE]

The expression in parentheses on the right hand side is monotonically decreasing in $R_{j}R_{j}^{T}$, so using (52) gives

[TABLE]

which holds with probability222This bound follows because the events $A_{j}$ that (52) holds are not necessarily independent, but their joint probability can be bounded like $P(A_{1}\cap\cdots\cap A_{n})=1-P(A_{1}^{C}\cup\cdots\cup A_{n}^{C})\geq 1-P(A_{1}^{C})-\cdots-P(A_{n}^{C})=1-n\beta$. at least $1-n\beta$.

In particular, taking $t=\sqrt{2\sigma_{e}^{2}\chi^{2}_{\beta}(\Sigma,\ I)\ln(2/\beta)}/N$ gives

[TABLE]

with probability at least $1-n\beta$, or equivalently,

[TABLE]

with probability at least $1-2n\beta$. $\square$

A.6 Proof of Proposition A.3

When $\hat{\theta}_{\epsilon_{N}}$ is the solution of the problem in (16), we have,

[TABLE]

Denote $e_{\epsilon_{N}}=Y_{N}-\Phi_{N}^{T}\hat{\theta}_{\epsilon_{N}}$, and $R_{j}$ as the $j^{th}$ row of the matrix $\Phi_{N}$, then (55) becomes,

[TABLE]

From Assumption 5.2, and using the same argument as in Proposition A.2, we have that each element of $R_{j}$ is distributed as $\mathcal{N}(0,\Sigma(j,j))$.

Consider the variable $Z=N^{-1}e_{\epsilon_{N}}^{T}R_{j}^{T}$, Since $R_{j}\sim\mathcal{N}(0,\ \Sigma)$,

[TABLE]

Since $\Sigma$ is symmetric and positive definite matrix, hence using singular value decomposition, we can find a diagonal matrix $D$ that satisfies,

[TABLE]

where $Q$ is the unitary matrix, i.e. $QQ^{T}=I$. Therefore, we have,

[TABLE]

where $s_{\text{max}}$ is the maximum element on the diagonal of matrix $D$, i.e. maximum singular value of matrix $\Sigma$. Note that,

[TABLE]

From Section 4.4 in [18], we have,

[TABLE]

which gives,

[TABLE]

with probability $1-\beta,\ 0\leq\beta\leq 1$. Combining this inequality with (57) and (58) gives, with probability $1-\beta$,

[TABLE]

Hence,

[TABLE]

with probability $1-\beta$. This means,

[TABLE]

with probability at least $1-n\beta$, following the same reasoning as in the proof of Proposition A.2. Therefore,

[TABLE]

with probability at least $1-n\beta$.

Taking $t=\sqrt{2\sigma_{e}^{2}s_{\text{max}}\chi^{2}_{\beta}(N-n)(1+\epsilon_{N})\ln(2/\beta)}/N$ gives

[TABLE]

with probability at least $1-n\beta$, or, equivalently,

[TABLE]

with probability at least $1-2n\beta$. $\square$