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 diverges with the dimension 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 . 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 , 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 to grow and possibly exceed the sample size . Of course, consistent estimation when is impossible unless the model is endowed with additional structure—for instance, sparsity in the parameter vector . Consequently, we study the class of models where has at most nonzero parameters, where is also allowed to increase to infinity with and .
2 M𝑀M-estimators for noisy and missing covariates
As long as the constraint radius is at least , the unique solution to this convex program is . Of course, this program is an idealization, since in practice we may not know the covariance matrix , and we certainly do not know —after all, is the quantity we are trying to estimate!
or alternatively, the regularized version
where 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 serve as unbiased estimators for .
Suppose we observe , where is a random matrix independent of , with rows drawn i.i.d. from a zero-mean distribution with known covariance . We consider the pair
where denotes elementwise division. A small calculation shows that these are unbiased estimators of and , respectively. The estimators (10) have been studied in past work XuYou07 , but only under classical scaling ().
As a special case of the estimators (10), suppose the entries of are independent random variables. Then the observed matrix corresponds to a missing-data matrix, where each element of the th column has probability of being missing. In this case, the estimators (10) become
is the parameter vector containing the ’s, and 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 satisfies a lower restricted eigenvalue condition with curvature and tolerance 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 satisfies an upper restricted eigenvalue condition with smoothness and tolerance 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 . In application to the constrained version, the method of projected gradient descent generates a sequence of iterates 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 satisfy the deviation bound (16), and the matrix satisfies the lower-RE condition (12) with parameters such that
Then for any vector with sparsity at most , there is a universal positive constant such that any global optimum of the Lagrangian program (7) with any satisfies the bounds
The same bounds (without ) also apply to the constrained program (4) with radius choice .
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 and of the program (4), Theorem 1 combined with the triangle inequality guarantees that
Finally, as noted by a reviewer, the constraint in the program (4) is rather restrictive, since 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 to be larger than , and the conditions on the regularizer 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 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 . In both cases, we assume that , as is required for statistical consistency in Theorem 1.
Under the conditions of Theorem 1: {longlist}[(a)]
For any global optimum of the constrained program (4), there are universal positive constants and a contraction coefficient , independent of , such that the gradient descent iterates (14) satisfy the bounds
Letting denote the objective function of Lagrangian program (7) with global optimum , and applying composite gradient updates (15), there are universal positive constants and a contraction coefficient , independent of , such that
where .
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 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 . We then reapplied the projected gradient method to the same problem instance times, each time with a random starting point, and measured the error between the iterates and the first estimate (optimization error), and the error between the iterates and the truth (statistical error). Within each panel, the blue traces show the optimization error over 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 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 , 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 , we form the pair such that and . In the proofs of Section 5, we will analyze the case where and show that the result of Corollary 1 still holds when 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 , we define , and assume .
Remarks. Suppose is a Gaussian random matrix and for all . In this case, the ratio is the condition number of . Then
a quantity that depends on both the conditioning of , and the fraction 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 is estimated by .
2.3 Bounds for dependent data
Turning to the case of dependent data, we consider the setting where the rows of are drawn from a stationary vector autoregressive (VAR) process according to
Note that we may extend the cases of dependent data to situations when and 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 is a vector of i.i.d. Gaussians and for each . If we define , we can verify that . Our algorithm, described below, forms estimates and for each , then combines the estimates to obtain an estimate .
In the additive noise case, we observe the matrix . From the equations (26), we obtain . Note that is a vector of i.i.d. Gaussians, and since , we have . Hence, our results on covariates with additive noise allow us to recover from . We can verify that this reduces to solving the program (4) or (7) with the pair , where .
When is a missing-data version of , we similarly estimate the vectors 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 , we use the pair
where {\widehat{\Sigma}}=\frac{1}{n}Z^{T}Z\mbox{\,{}:\ominus{}\,}M, and is defined as in Example 3.
To obtain the estimate , we therefore propose the following procedure, based on the estimators and .
(1) Perform linear regressions of the variables upon the remaining variables , using the program (4) or (7) with the estimators , to obtain estimates of .
(2) Estimate the scalars using the quantity , based on the estimator . Form with and .
(3) Set , where 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 :
Suppose the columns of the matrix are -sparse, and suppose the condition number is nonzero and finite. Suppose we have
and suppose we have the following additional deviation condition on :
Finally, suppose the lower-RE condition holds uniformly over the matrices with the scaling (18). Then under the estimation procedure of Algorithm 3.1, there exists a universal constant 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 ; and the lower-RE condition holds over all matrices by the same technique used to establish the lower-RE condition for . 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 when , for the appropriate values of and .
Simulations
In order to verify this theoretical prediction, we plotted versus the rescaled error . As shown by Figure 4(a), the curve is roughly constant, as predicted by the theory.
The plot of versus the rescaled error 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 denote the loss function to be minimized. This definition captures both the estimator (4) with and the estimator (7) with the choice of given in the theorem statement. For either estimator, we are guaranteed that is feasible and is optimal for the program, so . Indeed, in the regularized case, the -sparsity of implies that . Defining the error vector 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 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 . 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 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 .
Proof of inequality (30)
We first consider the constrained program (4), with , so . Combined with inequality (5.1), we conclude that . Consequently, we have the inequality , which is a slightly stronger form of the bound (30).
For the regularized estimator (7), we first note that our choice of guarantees that the term (33) is at most . 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 . Combining the pieces, we conclude that
and rearranging implies , from which we conclude that , 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 .
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 defined there; since is supported on the set with and the RE conditions hold with , we find that
where the final inequality makes use of the assumption that . Similarly, we may compute the contraction coefficient to be
so for .
combining the bounds yields . Then
Turning to the Lagrangian version, we exploit Theorem 2 in Agarwal et al. AgaEtal11 , with corresponding to the subspace of all vectors with support contained within the support set of . With this choice, we have , and the contraction coefficient takes the previous form (35), so that the assumption guarantees that . 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 and , and the condition implies that . Putting together the pieces, we find that the compound tolerance parameter satisfies the bound , 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 .