High-dimensional regression with noisy and missing data: Provable guarantees with nonconvexity

Po-Ling Loh, Martin J. Wainwright

Introduction

In standard formulations of prediction problems, it is assumed that the covariates are fully-observed and sampled independently from some underlying distribution. However, these assumptions are not realistic for many applications, in which covariates may be observed only partially, observed subject to corruption or exhibit some type of dependency. Consider the problem of modeling the voting behavior of politicians: in this setting, votes may be missing due to abstentions, and temporally dependent due to collusion or “tit-for-tat” behavior. Similarly, surveys often suffer from the missing data problem, since users fail to respond to all questions. Sensor network data also tends to be both noisy due to measurement error, and partially missing due to failures or drop-outs of sensors.

There are a variety of methods for dealing with noisy and/or missing data, including various heuristic methods, as well as likelihood-based methods involving the expectation–maximization (EM) algorithm (e.g., see the book LitRub87 and references therein). A challenge in this context is the possible nonconvexity of associated optimization problems. For instance, in applications of EM, problems in which the negative likelihood is a convex function often become nonconvex with missing or noisy data. Consequently, although the EM algorithm will converge to a local minimum, it is difficult to guarantee that the local optimum is close to a global minimum.

In this paper, we study these issues in the context of high-dimensional sparse linear regression—in particular, in the case when the predictors or covariates are noisy, missing, and/or dependent. Our main contribution is to develop and study simple methods for handling these issues, and to prove theoretical results about both the associated statistical error and the optimization error. Like EM-based approaches, our estimators are based on solving optimization problems that may be nonconvex; however, despite this nonconvexity, we are still able to prove that a simple form of projected gradient descent will produce an output that is “sufficiently close”—as small as the statistical error—to any global optimum. As a second result, we bound the statistical error, showing that it has the same scaling as the minimax rates for the classical cases of perfectly observed and independently sampled covariates. In this way, we obtain estimators for noisy, missing, and/or dependent data that have the same scaling behavior as the usual fully-observed and independent case. The resulting estimators allow us to solve the problem of high-dimensional Gaussian graphical model selection with missing data.

There is a large body of work on the problem of corrupted covariates or error-in-variables for regression problems (e.g., see the papers and books Hwa86 , CarEtal95 , ItuEtal99 , XuYou07 , as well as references therein). Much of the earlier theoretical work is classical in nature, meaning that it requires that the sample size n{n} diverges with the dimension p{p} fixed. Most relevant to this paper is more recent work that has examined issues of corrupted and/or missing data in the context of high-dimensional sparse linear models, allowing for n≪p{n}\ll{p}. Städler and Bühlmann StaBuh10 developed an EM-based method for sparse inverse covariance matrix estimation in the missing data regime, and used this result to derive an algorithm for sparse linear regression with missing data. As mentioned above, however, it is difficult to guarantee that EM will converge to a point close to a global optimum of the likelihood, in contrast to the methods studied here. Rosenbaum and Tsybakov RosTsy10 studied the sparse linear model when the covariates are corrupted by noise, and proposed a modified form of the Dantzig selector (see the discussion following our main results for a detailed comparison to this past work, and also to concurrent work RosTsy11 by the same authors). For the particular case of multiplicative noise, the type of estimator that we consider here has been studied in past work XuYou07 ; however, this theoretical analysis is of the classical type, holding only for n≫p{n}\gg{p}, in contrast to the high-dimensional models that are of interest here.

The remainder of this paper is organized as follows. We begin in Section 2 with background and a precise description of the problem. We then introduce the class of estimators we will consider and the form of the projected gradient descent algorithm. Section 3 is devoted to a description of our main results, including a pair of general theorems on the statistical and optimization error, and then a series of corollaries applying our results to the cases of noisy, missing, and dependent data. In Section 4, we demonstrate simulations to confirm that our methods work in practice, and verify the theoretically-predicted scaling laws. Section 5 contains proofs of some of the main results, with the remaining proofs contained in the supplementary Appendix LohWai11 .

Background and problem setup

In this section, we provide background and a precise description of the problem, and then motivate the class of estimators analyzed in this paper. We then discuss a simple class of projected gradient descent algorithms that can be used to obtain an estimator.

This setup applies to various disturbances to the covariates, including: {longlist}[(a)]

We work within a high-dimensional framework that allows the number of predictors p{p} to grow and possibly exceed the sample size n{n}. Of course, consistent estimation when n≪p{n}\ll{p} is impossible unless the model is endowed with additional structure—for instance, sparsity in the parameter vector β∗\beta^{*}. Consequently, we study the class of models where β∗\beta^{*} has at most k{k} nonzero parameters, where k{k} is also allowed to increase to infinity with pp and nn.

2 M𝑀M-estimators for noisy and missing covariates

As long as the constraint radius RR is at least ∥β∗∥1\|\beta^{*}\|_{1}, the unique solution to this convex program is β^=β∗{\widehat{\beta}}=\beta^{*}. Of course, this program is an idealization, since in practice we may not know the covariance matrix Σx{\Sigma_{x}}, and we certainly do not know Σxβ∗{\Sigma_{x}}\beta^{*}—after all, β∗\beta^{*} is the quantity we are trying to estimate!

or alternatively, the regularized version

where λn>0{\lambda_{n}}>0 is a user-defined regularization parameter. Note that the two problems are equivalent by Lagrangian duality when the objectives are convex, but not in the case of a nonconvex objective. The Lasso Tib96 , CheEtal98 is a special case of these programs, obtained by setting

In the presence of nonconvexity, it is generally impossible to provide a polynomial-time algorithm that converges to a (near) global optimum, due to the presence of local minima. Remarkably, we are able to prove that this issue is not significant in our setting, and a simple projected gradient descent algorithm applied to the programs (4) or (7) converges with high probability to a vector extremely close to any global optimum.

Let us illustrate these ideas with some examples. Recall that (Γ^,γ^)({\widehat{\Gamma}},\widehat{\gamma}) serve as unbiased estimators for (Σx,Σxβ∗)(\Sigma_{x},\Sigma_{x}\beta^{*}).

Suppose we observe Z=X+WZ=X+W, where WW is a random matrix independent of XX, with rows wiw_{i} drawn i.i.d. from a zero-mean distribution with known covariance Σw{\Sigma_{w}}. We consider the pair

where :⊖{}:\ominus{} denotes elementwise division. A small calculation shows that these are unbiased estimators of Σx\Sigma_{x} and Σxβ∗\Sigma_{x}\beta^{*}, respectively. The estimators (10) have been studied in past work XuYou07 , but only under classical scaling (n≫p{n}\gg{p}).

As a special case of the estimators (10), suppose the entries uiju_{ij} of UU are independent Bernoulli⁡(1−ρj)\operatorname{Bernoulli}(1-\rho_{j}) random variables. Then the observed matrix Z=X⊙UZ=X\odot U corresponds to a missing-data matrix, where each element of the jjth column has probability ρj\rho_{j} of being missing. In this case, the estimators (10) become

\boldsρ{\bolds\rho} is the parameter vector containing the ρj\rho_{j}’s, and 1\mathbf{1} is the vector of all 1’s. In this way, we obtain a generalization of the estimator discussed in Example 2.

3 Restricted eigenvalue conditions

The matrix Γ^{\widehat{\Gamma}} satisfies a lower restricted eigenvalue condition with curvature α1>0{\alpha_{1}}>0 and tolerance τ(n,p)>0{\tau}({n},{p})>0 if

Finally, although such upper bounds are not necessary for statistical consistency, our algorithmic results make use of the analogous upper restricted eigenvalue condition, formalized in the following:

The matrix Γ^{\widehat{\Gamma}} satisfies an upper restricted eigenvalue condition with smoothness α2>0{\alpha_{2}}>0 and tolerance τ(n,p)>0{\tau}({n},{p})>0 if

4 Gradient descent algorithms

In addition to proving results about the global minima of the (possibly nonconvex) programs (4) and (5), we are also interested in polynomial-time procedures for approximating such optima. In this paper, we analyze some simple algorithms for solving either the constrained program (4) or the Lagrangian version (7). Note that the gradient of the quadratic loss function takes the form ∇L(β)=Γ^β−γ^\nabla{\mathcal{L}}(\beta)={\widehat{\Gamma}}\beta-\widehat{\gamma}. In application to the constrained version, the method of projected gradient descent generates a sequence of iterates {βt,t=0,1,2,…}\{\beta^{t},t=0,1,2,\ldots\} by the recursion

Main results and consequences

We now state our main results and discuss their consequences for noisy, missing, and dependent data.

To aid intuition, note that inequality (16) holds whenever the following two deviation conditions are satisfied:

Suppose the surrogates (Γ^,γ^)({\widehat{\Gamma}},\widehat{\gamma}) satisfy the deviation bound (16), and the matrix Γ^{\widehat{\Gamma}} satisfies the lower-RE condition (12) with parameters (α1,τ)({\alpha_{1}},{\tau}) such that

Then for any vector β∗\beta^{*} with sparsity at most k{k}, there is a universal positive constant c0{c}_{0} such that any global optimum β^{\widehat{\beta}} of the Lagrangian program (7) with any b0≥∥β∗∥2b_{0}\geq\|\beta^{*}\|_{2} satisfies the bounds

The same bounds (without λn{\lambda_{n}}) also apply to the constrained program (4) with radius choice R=∥β∗∥1R=\|\beta^{*}\|_{1}.

Note that in the presence of nonconvexity, it is possible in principle for the optimization problems (4) and (7) to have many global optima that are separated by large distances. Interestingly, Theorem 1 guarantees that this unpleasant feature does not arise under the stated conditions: given any two global optima β^{\widehat{\beta}} and β~{\widetilde{\beta}} of the program (4), Theorem 1 combined with the triangle inequality guarantees that

Finally, as noted by a reviewer, the constraint R=∥β∗∥1R=\|\beta^{*}\|_{1} in the program (4) is rather restrictive, since β∗\beta^{*} is unknown. Theorem 1 merely establishes a heuristic for the scaling expected for this optimal radius. In this regard, the Lagrangian estimator (7) is more appealing, since it only requires choosing b0b_{0} to be larger than ∥β∗∥2\|\beta^{*}\|_{2}, and the conditions on the regularizer λn{\lambda_{n}} are the standard ones from past work on the Lasso.

1.2 Optimization error

Although Theorem 1 provides guarantees that hold uniformly for any global minimizer, it does not provide guidance on how to approximate such a global minimizer using a polynomial-time algorithm. Indeed, for nonconvex programs in general, gradient-type methods may become trapped in local minima, and it is impossible to guarantee that all such local minima are close to a global optimum. Nonetheless, we are able to show that for the family of programs (4), under reasonable conditions on Γ^{\widehat{\Gamma}} satisfied in various settings, simple gradient methods will converge geometrically fast to a very good approximation of any global optimum. The following theorem supposes that we apply the projected gradient updates (14) to the constrained program (4), or the composite updates (15) to the Lagrangian program (7), with stepsize η=2α2\eta=2{\alpha_{2}}. In both cases, we assume that n≿klog⁡p{n}\succsim{k}\log{p}, as is required for statistical consistency in Theorem 1.

Under the conditions of Theorem 1: {longlist}[(a)]

For any global optimum β^{\widehat{\beta}} of the constrained program (4), there are universal positive constants (c1,c2)({c}_{1},{c}_{2}) and a contraction coefficient γ∈(0,1){\gamma}\in(0,1), independent of (n,p,k)(n,p,k), such that the gradient descent iterates (14) satisfy the bounds

Letting ϕ\phi denote the objective function of Lagrangian program (7) with global optimum β^{\widehat{\beta}}, and applying composite gradient updates (15), there are universal positive constants (c1,c2)({c}_{1},{c}_{2}) and a contraction coefficient γ∈(0,1){\gamma}\in(0,1), independent of (n,p,k)({n},{p},{k}), such that

where T:=c2log⁡(ϕ(β0)−ϕ(β^))δ2/log⁡(1/γ)T:=c_{2}\log\frac{(\phi(\beta^{0})-\phi({\widehat{\beta}}))}{\delta^{2}}/\log(1/{\gamma}).

Remarks. As with Theorem 1, these claims are deterministic in nature. Probabilistic conditions will enter into the corollaries, which involve proving that the surrogate matrices Γ^{\widehat{\Gamma}} used for noisy, missing and/or dependent data satisfy the lower- and upper-RE conditions with high probability. The proof of Theorem 2 itself is based on an extension of a result due to Agarwal et al. AgaEtal11 on the convergence of projected gradient descent and composite gradient descent in high dimensions. Their result, as originally stated, imposed convexity of the loss function, but the proof can be modified so as to apply to the nonconvex loss functions of interest here. As noted following Theorem 1, all global minimizers of the nonconvex program (4) lie within a small ball. In addition, Theorem 2 guarantees that the local minimizers also lie within a ball of the same magnitude. Note that in order to show that Theorem 2 can be applied to the specific statistical models of interest in this paper, a considerable amount of technical analysis remains in order to establish that its conditions hold with high probability.

Experimentally, we have found that the predictions of Theorem 2 are borne out in simulations. Figure 2 shows the results of applying the projected gradient descent method to solve the optimization problem (4) in the case of additive noise [panel (a)], and missing data [panel (b)]. In each case, we generated a random problem instance, and then applied the projected gradient descent method to compute an estimate β^{\widehat{\beta}}. We then reapplied the projected gradient method to the same problem instance 1010 times, each time with a random starting point, and measured the error ∥βt−β^∥2\|\beta^{t}-\widehat{\beta}\|_{2} between the iterates and the first estimate (optimization error), and the error ∥βt−β∗∥2\|\beta^{t}-\beta^{*}\|_{2} between the iterates and the truth (statistical error). Within each panel, the blue traces show the optimization error over 1010 trials, and the red traces show the statistical error. On the logarithmic scale given, a geometric rate of convergence corresponds to a straight line. As predicted by Theorem 2, regardless of the starting point, the iterates {βt}\{\beta^{t}\} exhibit geometric convergence to the same fixed point.To be precise, Theorem 2 states that the iterates will converge geometrically to a small neighborhood of all the global optima. The statistical error contracts geometrically up to a certain point, then flattens out.

2 Some consequences

We begin with the case of i.i.d. samples with additive noise, as described in Example 1.

(b) We may also compare the results in (a) with bounds from past work on high-dimensional sparse regression with noisy covariates RosTsy11 . In this work, Rosenbaum and Tsybakov derive similar concentration bounds on sub-Gaussian matrices. The tolerance parameters are all O(log⁡pn){\mathcal{O}}(\sqrt{\frac{\log p}{n}}), with prefactors depending on the sub-Gaussian parameters of the matrices. In particular, in their notation,

leading to the bound (cf. Theorem 2 of Rosenbaum and Tsybakov RosTsy11 )

Based on the estimator Σ^w{\widehat{\Sigma}}_{w}, we form the pair (Γ~,γ~)({\widetilde{\Gamma}},\widetilde{\gamma}) such that γ~=1nZTy\widetilde{\gamma}=\frac{1}{n}Z^{T}y and Γ~=ZTZn−Σ^w{\widetilde{\Gamma}}=\frac{Z^{T}Z}{n}-{\widehat{\Sigma}}_{w}. In the proofs of Section 5, we will analyze the case where Σ^w=1nW0TW0{\widehat{\Sigma}}_{w}=\frac{1}{n}W_{0}^{T}W_{0} and show that the result of Corollary 1 still holds when Σw\Sigma_{w} must be estimated from the data. Note that the estimator in equation (23) will also yield the same result, but the analysis is more complicated.

2.2 Bounds for missing data: i.i.d. case

Next, we turn to the case of i.i.d. samples with missing data, as discussed in Example 3. For a missing data parameter vector \boldsρ{\bolds\rho}, we define ρmax⁡:=max⁡jρj\rho_{\max}:=\max_{j}\rho_{j}, and assume ρmax⁡<1\rho_{\max}<1.

Remarks. Suppose XX is a Gaussian random matrix and ρj=ρ\rho_{j}=\rho for all jj. In this case, the ratio σx2λmin⁡(Σx)=λmax⁡(Σx)λmin⁡(Σx)=κ(Σx)\frac{\sigma_{x}^{2}}{\lambda_{\min}(\Sigma_{x})}=\frac{\lambda_{\max}(\Sigma_{x})}{\lambda_{\min}(\Sigma_{x})}=\kappa(\Sigma_{x}) is the condition number of Σx\Sigma_{x}. Then

a quantity that depends on both the conditioning of Σx\Sigma_{x}, and the fraction ρ∈[0,1)\rho\in[0,1) of missing data. We will consider the results of Corollary 2 applied to this example in the simulations of Section 4.

We will show in Section 5 that Corollary 2 holds when \boldsρ{\bolds\rho} is estimated by \boldsρ^{\widehat{\bolds\rho}}.

2.3 Bounds for dependent data

Turning to the case of dependent data, we consider the setting where the rows of XX are drawn from a stationary vector autoregressive (VAR) process according to

Note that we may extend the cases of dependent data to situations when Σw\Sigma_{w} and \boldsρ{\bolds\rho} are unknown and must be estimated from the data. The proofs of these extensions are identical to the i.i.d case, so we will omit them.

3 Application to graphical model inverse covariance estimation

where εj\varepsilon^{j} is a vector of i.i.d. Gaussians and εj⊥ ⁣ ⁣ ⁣ ⁣⊥X−j\varepsilon^{j}\perp\!\!\!\!\perp X^{-j} for each jj. If we define aj:=−(Σjj−Σj,−jθj)−1a_{j}:=-(\Sigma_{jj}-\Sigma_{j,-j}\theta^{j})^{-1}, we can verify that Θj,−j=ajθj\Theta_{j,-j}=a_{j}\theta^{j}. Our algorithm, described below, forms estimates θ^j\widehat{\theta}{}^{j} and a^j{\widehat{a}}_{j} for each jj, then combines the estimates to obtain an estimate Θ^j,−j=a^jθ^j{\widehat{\Theta}}_{j,-j}={\widehat{a}}_{j}\widehat{\theta}{}^{j}.

In the additive noise case, we observe the matrix Z=X+WZ=X+W. From the equations (26), we obtain Zj=X−jθj+(εj+Wj)Z^{j}=X^{-j}\theta^{j}+(\varepsilon^{j}+W^{j}). Note that δj=εj+Wj\delta^{j}=\varepsilon^{j}+W^{j} is a vector of i.i.d. Gaussians, and since X⊥ ⁣ ⁣ ⁣ ⁣⊥WX\perp\!\!\!\!\perp W, we have δj⊥ ⁣ ⁣ ⁣ ⁣⊥X−j\delta^{j}\perp\!\!\!\!\perp X^{-j}. Hence, our results on covariates with additive noise allow us to recover θj\theta^{j} from ZZ. We can verify that this reduces to solving the program (4) or (7) with the pair (Γ^(j),γ^(j))=(Σ^−j,−j,1nZ−jTZj)({\widehat{\Gamma}}^{(j)},\widehat{\gamma}{}^{(j)})=({\widehat{\Sigma}}_{-j,-j},\frac{1}{n}Z^{-jT}Z^{j}), where Σ^=1nZTZ−Σw{\widehat{\Sigma}}=\frac{1}{n}Z^{T}Z-\Sigma_{w}.

When ZZ is a missing-data version of XX, we similarly estimate the vectors θj\theta^{j} via equation (26), using our results on the Lasso with missing covariates. Here, both covariates and responses are subject to missing data, but this makes no difference in our theoretical results. For each jj, we use the pair

where {\widehat{\Sigma}}=\frac{1}{n}Z^{T}Z\mbox{\,{}:\ominus{}\,}M, and MM is defined as in Example 3.

To obtain the estimate Θ^{\widehat{\Theta}}, we therefore propose the following procedure, based on the estimators {(Γ^(j),γ^(j))}j=1p\{({\widehat{\Gamma}}^{(j)},\widehat{\gamma}{}^{(j)})\}_{j=1}^{p} and Σ^{\widehat{\Sigma}}.

(1) Perform p{p} linear regressions of the variables ZjZ^{j} upon the remaining variables Z−jZ^{-j}, using the program (4) or (7) with the estimators (Γ^(j),γ^(j))({\widehat{\Gamma}}^{(j)},\widehat{\gamma}{}^{(j)}), to obtain estimates θ^j\widehat{\theta}{}^{j} of θj\theta^{j}.

(2) Estimate the scalars aja_{j} using the quantity a^j:=−(Σ^jj−Σ^j,−jθ^j)−1{\widehat{a}}_{j}:=-({\widehat{\Sigma}}_{jj}-{\widehat{\Sigma}}_{j,-j}\widehat{\theta}{}^{j})^{-1}, based on the estimator Σ^{\widehat{\Sigma}}. Form Θ~{\widetilde{\Theta}} with Θ~j,−j=a^jθ^j{\widetilde{\Theta}}_{j,-j}={\widehat{a}}_{j}\widehat{\theta}{}^{j} and Θ~jj=−a^j{\widetilde{\Theta}}_{jj}=-{\widehat{a}}_{j}.

(3) Set Θ^=arg⁡min⁡Θ∈Sp∣ ⁣∣ ⁣∣Θ−Θ~∣ ⁣∣ ⁣∣1{\widehat{\Theta}}=\arg\min_{\Theta\in S^{p}}|\!|\!|\Theta-{\widetilde{\Theta}}|\!|\!|_{{1}}, where SpS^{p} is the set of symmetric matrices.

Note that the minimization in step (3) is a linear program, so is easily solved with standard methods. We have the following corollary about Θ^{\widehat{\Theta}}:

Suppose the columns of the matrix Θ\Theta are kk-sparse, and suppose the condition number κ(Θ)\kappa(\Theta) is nonzero and finite. Suppose we have

and suppose we have the following additional deviation condition on Σ^{\widehat{\Sigma}}:

Finally, suppose the lower-RE condition holds uniformly over the matrices Γ^(j){\widehat{\Gamma}}^{(j)} with the scaling (18). Then under the estimation procedure of Algorithm 3.1, there exists a universal constant c0c_{0} such that

Note that Corollary 5 is again a deterministic result, with parallel structure to Theorem 1. Furthermore, the deviation bounds (27) and (28) hold for all scenarios considered in Section 3.2 above, using Corollaries 1–4 for the first two inequalities, and a similar bounding technique for ∥Σ^−Σ∥max⁡\|{\widehat{\Sigma}}-\Sigma\|_{\max}; and the lower-RE condition holds over all matrices Γ^(j){\widehat{\Gamma}}^{(j)} by the same technique used to establish the lower-RE condition for Γ^{\widehat{\Gamma}}. The uniformity of the lower-RE bound over all sub-matrices holds because

Hence, the error bound in Corollary 5 holds with probability at least 1−c1exp⁡(−c2log⁡p)1-c_{1}\exp(-c_{2}\log p) when n≿klog⁡pn\succsim{k}\log{p}, for the appropriate values of φ\varphi and α1{\alpha_{1}}.

Simulations

In order to verify this theoretical prediction, we plotted σw\sigma_{w} versus the rescaled error 1+σw2σw+0.5∥β^−β∗∥2\frac{\sqrt{1+\sigma_{w}^{2}}}{\sigma_{w}+0.5}\|{\widehat{\beta}}-\beta^{*}\|_{2}. As shown by Figure 4(a), the curve is roughly constant, as predicted by the theory.

The plot of ρ{\rho} versus the rescaled error ∥β^−β∗∥21+0.5(1−ρ)\frac{\|{\widehat{\beta}}-\beta^{*}\|_{2}}{1+0.5(1-{\rho})} is shown in Figure 4(b). The curve is again roughly constant, agreeing with theoretical results.

Finally, we studied the behavior of the inverse covariance matrix estimation algorithm on three types of Gaussian graphical models: {longlist}[(a)]

Proofs

In this section, we prove our two main theorems. For the more technical proofs of the corollaries, see the supplementary Appendix LohWai11 .

Let L(β)=12βTΓ^β−⟨γ^,β⟩+λn∥β∥1{\mathcal{L}}(\beta)=\frac{1}{2}\beta^{T}{\widehat{\Gamma}}\beta-\langle\widehat{\gamma},\beta\rangle+{\lambda_{n}}\|\beta\|_{1} denote the loss function to be minimized. This definition captures both the estimator (4) with λn=0{\lambda_{n}}=0 and the estimator (7) with the choice of λn{\lambda_{n}} given in the theorem statement. For either estimator, we are guaranteed that β∗\beta^{*} is feasible and β^{\widehat{\beta}} is optimal for the program, so L(β^)≤L(β∗){\mathcal{L}}({\widehat{\beta}})\leq{\mathcal{L}}(\beta^{*}). Indeed, in the regularized case, the k{k}-sparsity of β∗\beta^{*} implies that ∥β∗∥1≤k∥β∗∥2≤b0k\|\beta^{*}\|_{1}\leq\sqrt{{k}}\|\beta^{*}\|_{2}\leq{b_{0}}\sqrt{{k}}. Defining the error vector ν^:=β^−β∗{\widehat{\nu}}:={\widehat{\beta}}-\beta^{*} and performing some algebra leads to the equivalent inequality

In the remainder of the proof, we first derive an upper bound for the right-hand side of this inequality. We then use this upper bound and the lower-RE condition to show that the error vector ν^{\widehat{\nu}} must satisfy the inequality

Finally, we combine inequality (30) with the lower-RE condition to derive a lower bound on the left-hand side of the basic inequality (29). Combined with our earlier upper bound on the right-hand side, some algebra yields the claim.

We first upper-bound the right-hand side of inequality (29). Hölder’s inequality gives ⟨ν^,γ^−Γ^β∗⟩≤∥ν^∥1∥γ^−Γ^β∗∥∞\langle{\widehat{\nu}},\widehat{\gamma}-{\widehat{\Gamma}}\beta^{*}\rangle\leq\|{\widehat{\nu}}\|_{1}\|\widehat{\gamma}-{\widehat{\Gamma}}\beta^{*}\|_{\infty}. By the triangle inequality, we have

where inequality (i) follows from the deviation conditions (3.1.1). Combining the pieces, we conclude that

where we have exploited the sparsity of β∗\beta^{*} and applied the triangle inequality. Combining the pieces, we conclude that the right-hand side of inequality (29) is upper-bounded by

a bound that holds for any nonnegative choice of λn{\lambda_{n}}.

Proof of inequality (30)

We first consider the constrained program (4), with R=∥β∗∥1R=\|\beta^{*}\|_{1}, so ∥β^∥1=∥β∗+ν^∥1≤∥β∗∥1\|{\widehat{\beta}}\|_{1}=\|\beta^{*}+{\widehat{\nu}}\|_{1}\leq\|\beta^{*}\|_{1}. Combined with inequality (5.1), we conclude that ∥ν^Sc∥1≤∥ν^S∥1\|{\widehat{\nu}}_{{{S}^{c}}}\|_{1}\leq\|{\widehat{\nu}}_{S}\|_{1}. Consequently, we have the inequality ∥ν^∥1≤2∥ν^S∥1≤2k∥ν^∥2\|{\widehat{\nu}}\|_{1}\leq 2\|{\widehat{\nu}}_{S}\|_{1}\leq 2\sqrt{{k}}\|{\widehat{\nu}}\|_{2}, which is a slightly stronger form of the bound (30).

For the regularized estimator (7), we first note that our choice of λn{\lambda_{n}} guarantees that the term (33) is at most 3λn2∥ν^S∥1−λn2∥ν^Sc∥1\frac{3{\lambda_{n}}}{2}\|{\widehat{\nu}}_{S}\|_{1}-\frac{{\lambda_{n}}}{2}\|{\widehat{\nu}}_{S^{c}}\|_{1}. Returning to the basic inequality, we apply the lower-RE condition to lower-bound the left-hand side, thereby obtaining the inequality

by our choice of λn{\lambda_{n}}. Combining the pieces, we conclude that

and rearranging implies ∥ν^Sc∥1≤7∥ν^S∥1\|{\widehat{\nu}}_{{{S}^{c}}}\|_{1}\leq 7\|{\widehat{\nu}}_{S}\|_{1}, from which we conclude that ∥ν^∥1≤8k∥ν^∥2\|{\widehat{\nu}}\|_{1}\leq 8\sqrt{{k}}\|{\widehat{\nu}}\|_{2}, as claimed.

Lower bound on left-hand side

We now derive a lower bound on the left-hand side of inequality (29). Combining inequality (30) with the RE condition (12) gives

where the final step uses our assumption that kτ(n,p)≤α1128{k}{\tau}({n},{p})\leq\frac{{\alpha_{1}}}{128}.

Finally, combining bounds (33), (30) and (34) yields

giving inequality (19a). Using inequality (30) again gives inequality (19b).

2 Proof of Theorem 2

In order to apply Theorem 1 in their paper, we first need to compute the tolerance parameter ε2\varepsilon^{2} defined there; since β∗\beta^{*} is supported on the set SS with ∣S∣=k|S|={k} and the RE conditions hold with τ≍log⁡pn\tau\asymp\frac{\log{p}}{{n}}, we find that

where the final inequality makes use of the assumption that n≿klog⁡p{n}\succsim{k}\log{p}. Similarly, we may compute the contraction coefficient to be

so γ∈(0,1)\gamma\in(0,1) for n≿klog⁡pn\succsim{k}\log{p}.

combining the bounds yields ∥ΔSct∥1≤∥ΔSt∥1+∥β^−β∗∥1\|\Delta^{t}_{S^{c}}\|_{1}\leq\|\Delta^{t}_{S}\|_{1}+\|{\widehat{\beta}}-\beta^{*}\|_{1}. Then

Turning to the Lagrangian version, we exploit Theorem 2 in Agarwal et al. AgaEtal11 , with M\mathcal{M} corresponding to the subspace of all vectors with support contained within the support set of β∗\beta^{*}. With this choice, we have ψ(M)=k\psi(\mathcal{M})=\sqrt{k}, and the contraction coefficient γ{\gamma} takes the previous form (35), so that the assumption n≿klog⁡pn\succsim{k}\log{p} guarantees that γ∈(0,1){\gamma}\in(0,1). It remains to verify that the requirements are satisfied. From the conditions in our Theorem 2 and using the notation of Agarwal et al. AgaEtal11 , we have β(M)=O(log⁡pn)\beta(\mathcal{M})={\mathcal{O}}(\frac{\log{p}}{{n}}) and ρ‾=k\overline{\rho}=\sqrt{k}, and the condition n≿klog⁡p{n}\succsim{k}\log{p} implies that ξ(M)=O(1)\xi(\mathcal{M})={\mathcal{O}}(1). Putting together the pieces, we find that the compound tolerance parameter ε2\varepsilon^{2} satisfies the bound ε2=O(klog⁡pn∥β^−β∗∥22)=O(∥β^−β∗∥22)\varepsilon^{2}={\mathcal{O}}(\frac{{k}\log{p}}{{n}}\|{\widehat{\beta}}-\beta^{*}\|_{2}^{2})={\mathcal{O}}(\|{\widehat{\beta}}-\beta^{*}\|_{2}^{2}), so the claim follows.

Discussion

Future directions of research include studying more general types of dependencies or corruption in the covariates of regression, such as more general types of multiplicative noise, and performing sparse linear regression for corrupted data with additive noise when the noise covariance is unknown and replicates of the data may be unavailable. As pointed out by a reviewer, it would also be interesting to study the performance of our algorithms on data that are not sub-Gaussian, or even under model mismatch. In addition, one might consider other loss functions, where it is more difficult to correct the objective for corrupted covariates. Finally, it remains to be seen whether or not our techniques—used to show that certain nonconvex problems can solved to statistical precision—can be applied more broadly.

Acknowledgments

The authors thank Alekh Agarwal, Sahand Negahban, John Duchi and Alexandre Tsybakov for useful discussions and guidance. They are also grateful to the Associate Editor and anonymous referees for improvements on the paper.

Supplementary material for: High-dimensional regression with noisy and missing data: Provable guarantees with nonconvexity \slink[doi]10.1214/12-AOS1018SUPP \sdatatype.pdf \sfilenameaos1018_supp.pdf \sdescriptionDue to space constraints, we have relegated technical details of the remaining proofs to the supplement LohWai11 .

References