Support recovery without incoherence: A case for nonconvex regularization

We demonstrate that the primal-dual witness proof method may be used to establish variable selection consistency and $\ell_\infty$-bounds for sparse regression problems, even when the loss function and/or regularizer are nonconvex. Using this method, we derive two theorems concerning support recovery and $\ell_\infty$-guarantees for the regression estimator in a general setting. Our results provide rigorous theoretical justification for the use of nonconvex regularization: For certain nonconvex regularizers with vanishing derivative away from the origin, support recovery consistency may be guaranteed without requiring the typical incoherence conditions present in $\ell_1$-based methods. We then derive several corollaries that illustrate the wide applicability of our method to analyzing composite objective functions involving losses such as least squares, nonconvex modified least squares for errors-in variables linear regression, the negative log likelihood for generalized linear models, and the graphical Lasso. We conclude with empirical studies to corroborate our theoretical predictions.


Introduction
The last two decades have generated a significant body of work involving convex relaxations of nonconvex problems arising in high-dimensional sparse regression (e.g., see the papers [10,36,7,5,3,40] and references therein). In broad terms, the underlying goal is to identify a relatively sparse solution from among a set of candidates that also yields a good fit to the data. A hard sparsity constraint is most directly encoded in terms of the ℓ 0 -"norm," which counts the number of nonzero entries in a vector. However, this produces a nonconvex optimization problem that may be NP-hard to solve or even approximate [27,38]. As a result, much work has focused instead on a slightly different problem, where the ℓ 0 -based constraint is replaced by the convex ℓ 1 -norm. There is a relatively well-developed theoretical understanding of the conditions under which such ℓ 1 -relaxations produce good estimates of the underlying parameter vector (e.g., see the papers [3,37,25] and references therein).
Although the ℓ 1 -norm encourages sparsity in the solution, however, it differs from the ℓ 0 -norm in a crucial aspect: Whereas the ℓ 0 -norm is constant for any nonzero argument, the ℓ 1 -norm increases linearly with the absolute value of the argument. This linear increase leads to a bias in the resulting ℓ 1 -regularized solution, and noticeably affects the performance of the estimator in finite-sample settings [11,4,21]. Motivated by this deficiency of ℓ 1regularization, several authors have proposed alternative forms of nonconvex regularization, including the smoothly clipped absolute deviation (SCAD) penalty [11], minimax concave penalty (MCP) [44], and log-sum penalty (LSP) [6]. These regularizers may be viewed as a hybrid of ℓ 0 -and ℓ 1 -regularizers-they resemble the ℓ 1 -norm in a neighborhood of the origin, but become (asymptotically) constant at larger values. Although the nonconvexity of the regularizer causes the overall optimization problem to be nonconvex, numerous empirical studies have shown that gradient-based optimization methods, while only guaranteed to find local optima, often produce estimators with consistently smaller estimation error than the estimators produced by the convex ℓ 1 -penalty [11,16,46,4,21].
In recent years, several important advances have been made toward developing a theoretical framework for nonconvex regularization. Zhang and Zhang [43] provide results showing that global optima of nonconvex regularized least squares problems are statistically consistent for the true regression vector, leaving open the question of how to find such optima efficiently. Fan et al. [13] show that one step of a local linear approximation (LLA) algorithm, initialized at a Lasso solution, results in a local optimum of the nonconvex regularized least squares problem that satisfies oracle properties; Wang et al. [41] establish similar guarantees for the output of a particular path-following algorithm. Our own past work [20] provides a general set of sufficient conditions under which all stationary points of the nonconvex regularized problem are guaranteed to lie within statistical precision of the true parameter, which substantially simplifies the optimization problem to one of finding stationary points. Our work also establishes bounds on the ℓ 2 -norm and prediction error that agree with the well-known bounds for the convex ℓ 1 -regularizer, up to constant factors.
Despite these advances, however, an important question has remained: Are stationary points of such nonconvex problems also consistent for variable selection? In other words, does the support set of a stationary point agree with the support of the true regression vector, with high probability, and at what rate does the error probability tend to zero? In addition to providing a natural venue for establishing bounds on the ℓ ∞ -error, support recovery results furnish a much better understanding of when stationary points of nonconvex objectives are actually unique. For convex objectives, various standard proof techniques for variable selection consistency now exist, including approaches via Karush-Kuhn-Tucker optimality conditions and primal-dual witness arguments (e.g., [45,40,18]). However, these arguments, as previously stated, have relied crucially on convexity of both the loss and regularizer.
The first main contribution of our paper is to show how the primal-dual witness technique may be modified and extended to a certain class of nonconvex problems. Our proof hinges on the notion of generalized gradients from nonsmooth analysis [8], and optimization-theoretic results on norm-regularized, smooth, but possibly nonconvex functions [14]. Our main result is to establish sufficient conditions for variable selection consistency when both the loss and regularizer are allowed to be nonconvex, provided the loss function satisfies a form of restricted strong convexity and the regularizer satisfies suitable mild conditions. These assumptions are similar to the conditions required in our earlier work on ℓ 1 -and ℓ 2 -consistency [20], with an additional assumption on the minimum signal strength that allows us to derive stronger support recovery guarantees. Remarkably, our results demonstrate that for a certain class of regularizers-including the SCAD and MCP regularizers-we may dispense with the usual incoherence conditions required by ℓ 1 -based methods, and still guarantee support recovery consistency for all stationary points of the resulting nonconvex program. This provides a strong theoretical reason for why certain nonconvex regularizers might be favored over their convex counterparts. In addition, we establish that for the same class of nonconvex regularizers, the unique stationary point is in fact equal to the oracle solution.
Several other authors have mentioned the potential for nonconvex regularizers to deliver estimation and support recovery guarantees under weaker assumptions than the ℓ 1 -norm. The same line of work introducing nonconvex penalties such as the SCAD and MCP and developing subsequent theory [11,12,46,42,13] demonstrates that in the absence of incoherence conditions, nonconvex regularized problems possess local optima that are statistically consistent and satisfy an oracle property. Since nonconvex programs may have multiple solutions, however, these papers have focused on establishing theoretical guarantees for the output of specific optimization algorithms. More recently, Wang et al. [41] propose a path-following homotopy algorithm for obtaining solutions to nonconvex regularized M -estimators, and show that iterates of the homotopy algorithm converge at a linear rate to the oracle solution of the M -estimation problem. In contrast to theory of this type-applicable to a particular algorithm-the theory in our paper is purely statistical and does not concern iterates of a particular optimization algorithm. Indeed, the novelty of our theoretical results is that they establish support recovery consistency for all stationary points and, moreover, shed light on situations where such stationary points are actually unique. Finally, Pan and Zhang [29] provide results showing that under restricted eigenvalue assumptions on the design matrix that are weaker than the standard restricted eigenvalue conditions, a certain class of nonconvex regularizers yield estimates that are consistent in ℓ 2 -norm. They provide bounds on the sparsity of approximate global and approximate sparse (AGAS) solutions, a notion also studied in earlier work [43]. However, their theoretical development stops short of providing conditions for recovering the exact support of the underlying regression vector. The remainder of our paper is organized as follows: In Section 2, we provide basic background material on regularized M -estimators and set up notation for the paper. We also outline the primal-dual witness proof method. Section 3 is devoted to the statements of our main results concerning support recovery and ℓ ∞ -bounds, followed by corollaries that specialize our results to linear regression, generalized linear models, and the graphical Lasso. In each case, we contrast our conditions for nonconvex regularizers to those required by convex regularizers and discuss the implications of our significantly weaker assumptions. We provide proofs of our main results in Appendices A and B, with supporting results and more technical lemmas contained in later appendices. Finally, Section 4 contains convergence guarantees for the composite gradient descent algorithm and a variety of illustrative simulations that confirm our theoretical results.
Notation: For functions f (n) and g(n), we write f (n) g(n) to mean that f (n) ≤ cg(n) for some universal constant c ∈ (0, ∞), and similarly, f (n) g(n) when f (n) ≥ c ′ g(n) for some universal constant c ′ ∈ (0, ∞). We write f (n) ≍ g(n) when f (n) g(n) and f (n) g(n) hold simultaneously. For a vector v ∈ R p and a subset S ⊆ {1, . . . , p}, we write v S ∈ R S to denote the vector v restricted to S. For a matrix M , we write |||M ||| 2 and |||M ||| F to denote the spectral and Frobenius norms, respectively, and write |||M ||| max := max i,j |m ij | to denote the elementwise ℓ ∞ -norm of M . For a function h : R p → R, we write ∇h to denote a gradient or subgradient, if it exists. Finally, for q, r > 0, we write B q (r) to denote the ℓ q -ball of radius r centered around 0.

Problem formulation
In this section, we briefly review the theory of regularized M -estimators. We also outline the primal-dual witness proof technique that underlies our proofs of variable consistency for nonconvex problems.

Regularized M-estimators
The analysis of this paper applies to regularized M -estimators of the form β ∈ arg min where L n denotes the empirical loss function and ρ λ denotes the penalty function, both assumed to be continuous. In our framework, both of these functions are allowed to be nonconvex, but the theory applies a fortiori when only one function is convex. The prototypical example of a loss function to keep in mind is the least squares objective, L n (β) = 1 2n y−Xβ 2 2 . We include the side constraint, β 1 ≤ R, in order to ensure that a global minimum β exists. 1 For modeling purposes, we have also allowed for an additional constraint, β ∈ Ω, where Ω is an open convex set; note that we may take Ω = R p when this extra constraint is not needed.
The analysis of this paper is restricted to the class of coordinate-separable regularizers, meaning that ρ λ is expressible as the sum Here, we have engaged in some minor abuse of notation; the functions ρ λ : R → R appearing on the right-hand side of equation (2) are univariate functions acting upon each coordinate.
Our results are readily extended to the inhomogenous case, where different coordinates have different regularizers ρ j λ , but we restrict ourselves to the homogeneous case in order to simplify our discussion.
From a statistical perspective, the purpose of solving the optimization problem (1) is to estimate the vector β * ∈ R p that minimizes the expected loss function: where we assume that β * is unique and independent of the sample size. Our goal is to develop sufficient conditions under which a minimizer β of the composite objective (1) is a consistent estimator for β * . Consequently, we will always choose R ≥ β * 1 , which ensures that β * is a feasible point.

Class of regularizers
In this paper, we study the class of regularizers ρ λ : R → R that are amenable in the following sense.
(iii) The function t → ρ λ (t) t is nonincreasing on R + . 1 In the sequel, we will give examples of nonconvex loss functions for which the global minimum fails to exist without such a side constraint (cf. Section 2.3 below).
We say that ρ λ is (µ, γ)-amenable if, in addition: (vii) There is some scalar γ ∈ (0, ∞) such that ρ ′ λ (t) = 0, for all t ≥ γλ. Conditions (vi) and (vii) are also known as the selection and unbiasedness properties, respectively, and the MCP regularizer [42] described below minimizes the maximum concavity of ρ subject to (vi)-(vii). Note that the usual ℓ 1 -penalty ρ λ (t) = λ|t| is 0-amenable, but it is not (0, γ)-amenable, for any γ < ∞. The notion of µ-amenability was also used in our past work on ℓ 2 -bounds for nonconvex regularizers [20], with the exception of the selection property (vi). Since the goal of the current paper is to obtain stronger conclusions, in terms of variable selection and ℓ ∞ -bounds, we will also require ρ λ to satisfy the selection and unbiasedness properties.
Many popular regularizers are either µ-amenable or (µ, γ)-amenable. Let us consider a few examples to illustrate.

Nonconvexity and restricted strong convexity
Next, we consider some examples of the types of nonconvex loss functions treated in this paper. At a high level, we consider loss functions that are differentiable and satisfy a form of restricted strong convexity, as used in large body of past work on analysis of high-dimensional sparse M -estimators (e.g., [1,3,37,20,25]). In order to provide intuition before stating the formal definition, note that for any convex and differentiable function f : R p → R that is globally convex and locally strongly convex around a point β ∈ R p , there exists a constant α > 0 such that for all ∆ ∈ R p . The notion of restricted strong convexity (with respect to the ℓ 1 -norm) weakens this requirement by adding a tolerance term that penalizes non-sparse vectors. In particular, for positive constants {(α j , τ j )} 2 j=1 , we have the following definition: Restricted strong convexity: Given any pair of vectors β, ∆ ∈ R p , the loss function L n satisfies an (α, τ )-RSC condition, if: where (α 1 , α 2 ) are strictly positive constants, and (τ 1 , τ 2 ) are nonnegative constants.
As noted in inequality (7), any locally strongly convex function that is also globally convex satisfies the RSC condition with tolerance parameters τ 1 = τ 2 = 0. For τ 1 , τ 2 > 0, the RSC condition imposes strong curvature only in certain directions of p-dimensional space-namely, those nonzero directions ∆ ∈ R p for which the ratio ∆ 1 ∆ 2 is relatively small; i.e., less than a constant multiple of n log p . Note that for any k-sparse vector ∆, we have ∆ 1 ∆ 2 ≤ √ k, so that the RSC definition guarantees a form of strong convexity for all k-sparse vectors when n k log p.
A line of past work (e.g., [30,34,25,20]) shows that the RSC condition holds, with high probability, for many types of convex and nonconvex objectives arising in statistical estimation problems. We now consider a few illustrative examples.
Standard linear regression: Consider the standard linear regression model, in which we observe i.i.d. pairs (x i , y i ) ∈ R p × R, linked by the linear model and the goal is to estimate β * ∈ R p . A standard loss function in this case is the least squares function L n (β) = 1 2n y − Xβ 2 2 , where y ∈ R n is the vector of responses and X ∈ R n×p is the design matrix with x i ∈ R p as its i th row. In this special case, for any β, ∆ ∈ R p , we have Consequently, for the least squares loss function, the RSC condition is essentially equivalent to lower-bounding sparse restricted eigenvalues [37,3].
Linear regression with errors in covariates: Let us now turn to a simple extension of the standard linear regression model. Suppose that instead of observing the covariates x i ∈ R p directly, we observe the corrupted vectors z i = x i + w i , where w i ∈ R p is some type of noise vector. This setup is a particular instantiation of a more general errors-in-variables model for linear regression. The standard Lasso estimate (applied to the observed pairs {(z i , y i )} n i=1 ) is inconsistent in this setting.
As studied previously [19], it is natural to consider a corrected version of the Lasso, which we state in terms of the quadratic objective, Our past work [19] shows that as long as ( Γ, γ) are unbiased estimates of (Σ x , Σ x β * ), any global minimizer β of the appropriately regularized problem (1) is a consistent estimate for β * . In the additive corruption model described in the previous paragraph, a natural choice for the pair ( Γ, γ) is given by where the covariance matrix Σ w = Cov(w i ) is assumed to be known. However, in the highdimensional setting (n ≪ p), the random matrix Γ is not positive semidefinite, so the quadratic objective function (9) is nonconvex. (This is also a concrete instance where the objective function (1) requires the constraint β 1 ≤ R in order to be bounded below.) Nonetheless, our past work [19,20] shows that under certain tail conditions on the covariates and noise vectors, the loss function (9) does satisfy a form of restricted strong convexity.
Generalized linear models: Moving beyond standard linear regression, suppose the pairs {(x i , y i )} n i=1 are drawn from a generalized linear model (GLM). Recall that the conditional distribution for a GLM takes the form where σ > 0 is a scale parameter and ψ is the cumulant function. The loss function corresponding to the negative log likelihood is given by and it is easy to see that equation (11) reduces to equation (9) with the choice ψ(t) = t 2 2 . Using properties of exponential families, we may also check that equation (3) holds. Negahban et al. [25] show that a form of restricted strong convexity holds for a broad class of generalized linear models.
Graphical Lasso: Now suppose that we observe a sequence {x i } n i=1 ⊆ R p of p-dimensional random vectors with mean 0. Our goal is to estimate the inverse covariance matrix Θ * := Σ −1 x , assumed to be relatively sparse. In the Gaussian case, sparse inverse covariance matrices arise from imposing a Markovian structure on the random vector [35]. Letting Σ := X T X n denote the sample covariance matrix, consider the loss function which we refer to as the graphical Lasso loss. Taking derivatives, it is easy to check that Θ * = arg min Θ E[L n (Θ)], which verifies the population-level condition (3). As we show in Section 3.5, the graphical Lasso loss is locally strongly convex, so it satisfies the restricted strong convexity condition with τ 1 = τ 2 = 0.

Primal-dual witness proof technique
We now outline the main steps of the primal-dual witness (PDW) proof technique, which we will use to establish support recovery and ℓ ∞ -bounds for the program (1). Such a technique was previously applied only in situations where L n is convex and ρ λ is the ℓ 1 -penalty, but we show that this machinery may be extended via a careful analysis of local optima of normregularized functions based on generalized gradients. As stated in Theorem 1 below, the success of the PDW construction guarantees that stationary points of the nonconvex objective are consistent for variable selection consistency-in fact, they are unique. Recall that β ∈ R p is a stationary point of the optimization program (1) if we have ∇L n ( β) + ∇ρ λ ( β), β − β ≥ 0, for all β in the feasible region [2]. Due to the possible nondifferentiability of ρ λ at 0, we abuse notation slightly and denote ∇ρ λ ( β), β − β = lim t→0 + ∇ρ λ ( β + t(β − β)), β − β (see, e.g., Clarke [8] for a more comprehensive treatment of such generalized gradients). The set of stationary points includes all local/global minima of the program (1), as well as any interior local maxima.
The key steps of the PDW argument are as follows: (i) Optimize the restricted program β S ∈ arg min where we enforce the additional constraint that supp( β S ) ⊆ supp(β * ) := S. Establish that β S 1 < R; i.e., β S is an interior point of the feasible set.
(iii) Show that β is a local minimum of the full program (1), and moreover, all stationary points of the program (1) are supported on S.
Note that the output ( β, z) of the PDW construction depends implicitly on the choice of λ and R. Under the restricted strong convexity condition, the restricted problem (13) minimized in step (i) is actually a convex program. Hence, if β S 1 < R, the zero-subgradient condition (14) must hold at β S for the restricted problem (13). Note that when L n is convex and ρ λ is the ℓ 1 -penalty as in the conventional setting, the additional ℓ 1 -constraint in the programs (1) and (13) is omitted. If also Ω = R p , the vector β S is automatically a zero-subgradient point if it is a global minimum of the restricted program (13), which greatly simplifies the analysis. Our refined analysis shows that under suitable restrictions, global optimality still holds for β S and β, and the convexity of the restricted program therefore implies uniqueness.
In the sections to follow, we show how the primal-dual witness technique may be used to establish support recovery results for general nonconvex regularized M -estimators, and then derive sufficient conditions under which stationary points of the program (1) are in fact unique.

Main statistical results and consequences
We begin by stating our main theorems concerning support recovery and ℓ ∞ -bounds, and then specialize our analysis to particular settings of interest.

Main results
Our main statistical results concern stationary points of the regularized M -estimator (1), where the loss function satisfies the RSC condition (8) with parameters {(α j , τ j )} 2 j=1 , and the regularizer is µ-amenable with µ ∈ [0, α 1 ). Our first theorem concerns the success of the PDW construction described in Section 2.4. The theorem guarantees support recovery provided two conditions are met, the first involving an appropriate choice of λ and R, and the second involving strict dual feasibility of the dual vector z. Note that it is through validating the second condition that the incoherence assumption arises in the usual ℓ 1 -analysis, but we demonstrate in our corollaries to follow that strict dual feasibility may be guaranteed under weaker conditions when a (µ, γ)-amenable regularizer is used, instead. (See Appendix C for a technical discussion.) The proof of Theorem 1 is contained in Appendix A.
Theorem 1 (PDW construction for nonconvex functions). Suppose L n is an (α, τ )-RSC function and ρ λ is µ-amenable, for some µ ∈ [0, α 1 ). Suppose that: (a) The parameters (λ, R) satisfy the bounds For some δ ∈ 4Rτ 1 log p nλ , 1 , the dual vector z from the PDW construction satisfies the strict dual feasibility condition Then for any k-sparse vector β * , the program (1) with a sample size n ≥ 2τ 1 2α 1 −µ k log p has a unique stationary point, given by the primal output β of the PDW construction.
Remark: Of course, Theorem 1 is vacuous unless proper choices of λ, R, and δ exist. In the corollaries to follow, we show that ∇L n (β * ) ∞ ≤ c log p n , with high probability, in many settings of interest. In particular, we may choose λ ≍ log p n to satisfy inequality (15a) when the sample size satisfies n k log p. Note that R ≍ 1 λ then causes inequality (15b) to be satisfied under the same sample size scaling. Indeed, if the RSC parameters were known, we could simply take R = min α 2 8λ , α 2 τ 2 n log p and then focus on tuning λ. Finally, note that the inequality 4Rτ 1 log p nλ ≤ 1 is satisfied as long as R ≤ nλ 4τ 1 log p , which is guaranteed by the preceding choice of (λ, R) and the scaling n k log p. In this way, the existence of an appropriate choice of δ is guaranteed. 2 Note also that our results require the assumption µ < 2α 1 , where a smaller gap of (2α 1 −µ) translates into a larger sample size requirement. This consideration may motivate an advantage of using the LSP regularizer over a regularizer such as SCAD or MCP; as discussed in Section 2.2, the SCAD and MCP regularizers have µ equal to a constant value, whereas µ = λ 2 → 0 for the LSP. On the other hand, the LSP is not (µ, γ)-amenable, which as discussed later, allows us to remove the incoherence condition for SCAD and MCP when establishing strict dual feasibility (16). Indeed, the MCP is designed so that µ is minimal subject to unbiasedness and selection of the regularizer [44]. This suggests that for more incoherent designs, the LSP may be preferred for variable selection, whereas for less incoherent designs, SCAD or MCP may be better. In the simulations of Section 4, however, the LSP regularizer only performs negligibly better than the ℓ 1 -penalty in situations where the incoherence condition holds and the same regularization parameter λ is chosen.
Finally, we note that although the conditions of Theorem 1 are already relatively mild, they are nonetheless sufficient conditions. Indeed, as confirmed experimentally, there are many situations where the condition µ < 2α 1 does not hold, yet the stationary points of the program (1) still appear to be supported on S and/or unique. Two feasible explanations for this phenomenon are the following: First, it is possible that in cases where 2α 1 ≤ µ, the composite objective function is not convex over the entire feasible set, yet β is still sufficiently close to β * and the positive definite condition holds at the point β. Examining the proof of Theorem 1, we may see that equation (17), together with strict dual feasibility, is sufficient for establishing that β is a local minimum of the program (1); Lemma 3 in Appendix A then implies that all stationary points are supported on S. Nonetheless, since the restricted program (13) is no longer guaranteed to be convex, multiple stationary points may exist. Second, we note that although strict convexity of the objective and a zero-subgradient condition are certainly sufficient conditions to guarantee a unique global minimum, a function may have diverse regions of convexity/concavity and still possess a unique global minimum. In such cases, we may still be lucky in simulation studies and obtain the global optimum.
Our second general theorem provides control on the ℓ ∞ -error between any stationary point and β * , and shows that the local/global optimum of a nonconvex regularized program agrees with the oracle result when the regularizer is (µ, γ)-amenable. We define the oracle estimator according to In other words, the oracle estimator is the unpenalized estimator obtained from minimizing L n over the true support set S. As shown in the proof, under the assumed RSC conditions, the restricted function L n | S is strictly convex and β O S is uniquely defined. With this notation, we have the following result: Theorem 2. Suppose the assumptions of Theorem 1 are satisfied and strict dual feasibility (16) holds. Then the unique stationary point β of the program (1) has the following properties: then β agrees with the oracle estimator β O , and we have the tighter bound The proof of Theorem 2 is provided in Appendix B.
Remark: Theorem 2 underscores the strength of (µ, γ)-amenable regularizers. Indeed, with the addition of a beta-min condition (19a), which provides a lower bound on the minimum signal strength, the unbiasedness property allows us to remove the second term in inequality (18) and obtain a faster oracle rate (19b). As described in greater detail in the corollaries below, we may show that the right-hand expression in inequality (19b) is bounded by O log p n , with high probability, provided the spectrum of ∇ 2 L n (β * ) is bounded appropriately.
We now unpack the consequences of Thoerems 1 and 2 for several concrete examples.

Ordinary least squares linear regression
Our first application focuses on the setting of ordinary least squares, together with the nonconvex regularizers introduced in Section 2.2: SCAD, MCP, and LSP. We compare the consequences of Theorems 1 and 2 for each of these regularizers with the corresponding results for the convex ℓ 1 -penalty. Our theory demonstrates a clear advantage of using nonconvex regularizers such as SCAD and MCP that are (µ, γ)-amenable; whereas support recovery based on ℓ 1 -based methods is known to require fairly stringent incoherence conditions, our corollaries show that methods based on nonconvex regularizers will guarantee support recovery even in the absence of incoherence conditions. The ρ λ -regularized form of least squares regression may be written in the form β ∈ arg min Note that the Hessian of the loss function is given by ∇ 2 L n (β) = X T X n . While the sample covariance matrix is always positive semidefinite, it has rank at most n. Hence, in highdimensional settings where n < p, the Hessian of the loss function has at least p − n zero eigenvalues, implying that any nonconvex regularizer ρ λ makes the overall program (20) nonconvex.
In analyzing the family of estimators (20), we assume throughout that n ≥ c 0 k log p, for a sufficiently large constant c 0 . By known information-theoretic results [39], this type of lower bound is required for any method to recover the support of a k-sparse signal, hence is not a limiting restriction. With this setup, we have the following result, proved in Appendix D.1: Corollary 1. Suppose X and ǫ are sub-Gaussian, and regularization parameters (λ, R) are chosen such that β * 1 ≤ R 2 and c ℓ log p n ≤ λ ≤ cu R , for some constants c ℓ and c u . Also suppose the sample covariance matrix Γ = X T X n satisfies the condition (a) Suppose ρ λ is µ-amenable, with µ < λ min (Σ x ), and Γ also satisfies the incoherence condition Then with probability at least 1 − c 1 exp(−c 2 min{k, log p}), the nonconvex objective (20) has a unique stationary point β (corresponding to the global optimum). Furthermore, supp( β) ⊆ supp(β * ), and Then with probability at least 1 − c 1 exp(−c 2 min{k, log p}), the nonconvex objective (20) has a unique stationary point β given by the oracle estimator β O , and Note that if we also have the beta-min condition β * min ≥ 2 c 3 log p n + c ∞ λ in part (a), then β is still a sign-consistent estimate of β * ; however, the guaranteed bound (23) is looser than the oracle bound (24) derived in part (b).
The proof of Corollary 1 is provided in Appendix D.1. Here, we make some comments about its consequences. Regularizers satisfying the conditions of part (b) include the SCAD and MCP penalties. Recall that for the SCAD penalty, we have µ = 1 a−1 ; and for the MCP, we have µ = 1 b . Hence, the lower-eigenvalue condition translates into 1 a−1 < λ min (Σ x ) and 1 b < λ min (Σ x ), respectively. The LSP penalty is an example of a regularizer that satisfies the conditions of part (a), but not part (b): with this choice, we have µ = λ 2 , so the condition µ < λ min (Σ x ) is satisfied asymptotically whenever λ min (Σ x ) is bounded below by a constant. A version of part (a) also holds for the ℓ 1 -penalty, as shown in past work [40]. A valuable consequence of Corollary 1 is that it establishes conditions under which stationary points are unique and variable selection consistency holds, when using certain nonconvex regularizers. The distinguishing point between parts (a) and (b) of the corollary is that using (µ, γ)-amenable regularizers allow us to do away with an incoherence assumption (22) and guarantee that the unique stationary point is in fact equal to the oracle estimator.
Furthermore, a great deal of past work on nonconvex regularizers [11,44,42,43] has focused on the ordinary least squares regression objective (20); hence, it is instructive to interpret the results of Corollary 1 in light of this existing work. Zhang [42] shows that the two-step MC+ estimator (beginning with a global optimum of the program (20) with the MCP regularizer) is guaranteed to be consistent for variable selection, under only a sparse eigenvalue assumption on the design matrix. Our result shows that the global optimum obtained in the MCP step is actually already guaranteed to be consistent for variable selection, provided we have only slightly stronger assumptions about lower-and upper-eigenvalue bounds on the design matrix. In another related paper, Wainwright [39] establishes necessary conditions for support recovery in a linear regression setting when the covariates are drawn from a Gaussian distribution. As remarked in that paper, the necessary conditions only require eigenvalue bounds on the design matrix, in contrast to the more stringent incoherence conditions appearing in necessary and sufficient conditions for the success of the Lasso [40,45]. Using standard matrix concentration results for sub-Gaussian variables, it may be shown that the inequalities (21) and (22) hold, with high probability, when the population-level bounds are satisfied: However, the second inequality (25) is a fairly strong assumption on the covariance matrix Σ x , and as we explore in the simulations of Section 4 below, simple covariance matrices such as the class (43) defined in Section 4.2 fail to satisfy the latter condition. Hence, Corollary 1 shows clear advantage of using the SCAD or MCP regularizers over the ℓ 1 -penalty or LSP when Σ x is not incoherent.

Linear regression with corrupted covariates
We now shift our focus to an application where the loss function itself is nonconvex. In particular, we analyze the situation when the general loss function is defined according to equation (9). To simplify our discussion, we only state an explicit corollary for the case when ρ λ is the convex ℓ 1 -penalty; the most general case, involving a nonconvex quadratic form and a nonconvex regularizer, is simply a hybrid of the analysis below and the arguments of the previous section. Our goal is to illustrate the applicability of the primal-dual witness technique for nonconvex loss functions.
Let us recall the problem of linear regression with corrupted covariates, as previously introduced in Section 2.3. The pairs {(x i , y i )} n i=1 are drawn according to the standard linear model y i = x T i β * + ǫ i . While the response vector y = {y i } n i=1 is assumed to be observed, suppose we observe only the corrupted versions z i = x i + w i of the covariates. Based on the observed variables {(z i , y i )} n i=1 , we may then compute the quantities and estimate β * based on the following nonconvex program: Also suppose n ≥ k 2 and n ≥ c 0 k log p, for a sufficiently large constant c 0 .
then with probability at least 1 − c 1 exp(−c 2 min{k, log p}), the nonconvex objective (27) has a unique stationary point β (corresponding to the global optimum) such that supp( β) ⊆ supp(β * ), and Note that if in addition, we have a lower bound of the form β * min ≥ 2 c 3 log p n + c ∞ λ , then we are guaranteed that β is sign-consistent for β * .
Corollary 2 may be understood as an extension of part (a) of Corollary 1: it shows how the primal-dual witness technique may be used even in a setting where the loss function is nonconvex. Under the same incoherence assumption (28) and the familiar sample size scaling n ≥ c 0 k log p of the usual Lasso, stationary points of the modified (nonconvex) Lasso program (27) are also support-recovery consistent. Corollary 2 also implies the rather surprising result that, although the objective (27) is indeed nonconvex whenever n < p and Σ w ≻ 0, it nonetheless has a unique stationary point that is in fact equal to the global optimum. This result further clarifies the simulation results appearing in Loh and Wainwright [19]. Indeed, those simulations are performed with the setting Γ = I p , so the incoherence condition (28) holds, with high probability, with η close to 0. A careful inspection of the plots in Figure  2 of Loh and Wainwright [19] confirms the theoretical conclusion of Corollary 2; more detailed simulations for non-identity assignments of Γ appear in Section 4 below. The proof of Corollary 2 is provided in Appendix D.2.

Generalized linear models
To further illustrate the power of nonconvex regularizers, we now move to the case where the loss function is the negative log likelihood of a generalized linear model. We show that the incoherence condition may again be removed if the regularizer ρ λ is (µ, γ)-amenable (as is the case for the SCAD and MCP regularizers).
For (10), we take Ω = R p and construct the composite objective β ∈ arg min We impose the following conditions on the covariates and the link function: (i) The covariates are uniformly bounded as x i ∞ ≤ M , for all i = 1, . . . , n.
The conditions of Assumption 1, although somewhat stringent, are nonetheless satisfied in various settings of interest. In particular, for logistic regression, we have ψ(t) = log(1+exp(t)), so and we may verify that the boundedness conditions in Assumption 1(ii) are satisfied with κ 2 = 0.25 and κ 3 = 0.1. Also note that the uniform bound on ψ ′′′ is used implicitly in the proof for support recovery consistency in the logistic regression analysis of Ravikumar et al. [31], whereas the uniform bound on ψ ′′ also appears in the conditions for ℓ 1 -and ℓ 2 -consistency in other past work [25,20]. The uniform boundedness condition in Assumption 1(i) is somewhat less desirable: although it always holds for categorical data, it does not hold for Gaussian covariates. We suspect that is possible to relax this constraint, but since our main goal is to illustrate the more general theory, we keep it here. In what follows, let Then with probability at least 1 − c 1 exp(−c 2 min{k, log p}), the nonconvex objective (30) has a unique stationary point β given by the oracle estimator β O , and It is worthwhile to compare Corollary 3 with the analysis of ℓ 1 -regularized logistic regression given in Ravikumar et al. [31] (see Theorem 1 in their paper). Both results require that the sample size is lower-bounded as n ≥ c 0 k 3 log p, but Ravikumar et al. [31] also require Q * to satisfy the incoherence condition As noted in their paper and by other authors, the incoherence condition (32) is difficult to interpret and verify for general GLMs. In contrast, Corollary 3 shows that by using a properly chosen nonconvex regularizer, this incoherence requirement may be removed entirely. In addition, Corollary 3 is attractive in its generality, since it applies to more than just the logistic case with an ℓ 1 -penalty and extends to various nonconvex problems where the uniqueness of stationary points is not evident a priori. The proof of Corollary 3 is contained in Appendix D.3.

Graphical Lasso
Finally, we discuss the consequences of our theorems for the graphical Lasso, as previously described in Section 2.3. Recall that for the graphical Lasso, the observations consist of a collection {x i } n i=1 of p-dimensional vectors, and the goal is to recover the support of the inverse covariance matrix Θ * = (Cov(X)) −1 . The analysis here is different and more subtle, because we seek a high-dimensional result in which the sample size scales only with the row sparsity of the inverse covariance matrix, as opposed to the total number of of nonzero parameters (which grows linearly with the matrix dimension, for any connected graph).
In order to prove such a result, we consider the constrained estimator where S p ++ denotes the convex cone of symmetric, strictly positive definite matrices, and we impose the spectral norm bound |||Θ||| 2 ≤ κ on the estimate. 3 A more standard choice would be to use the ℓ 1 -norm bound i,j |Θ ij | ≤ R as the sideconstraint, as done in our past work on Frobenius norm bounds for the graphical Lasso with nonconvex regularizers [20]. However, as we will show here, the formulation (33) actually leads to variable selection consistency results under the milder scaling n d 2 log p, rather than the scaling n s log p obtained in our analysis on Frobenius norm bounds [20]. Here, we use d to denote the maximum number of nonzeros in any row/column of Θ * , and s := | supp(Θ * )| to denote the total number of nonzero entries.
A few remarks are in order. First, when ρ λ is the convex ℓ 1 -penalty and κ = ∞, the program (33) is identical to the standard graphical Lasso (e.g., [9,32,15,33]). However, since we are interested in scenarios where ρ λ is allowed to be nonconvex, we include an additional spectral norm constraint governed by κ. As a technical comment, note that our original assumptions required Ω to be an open subset of R p . In the analysis to follow, we handle the symmetry constraint Θ = Θ T by treating the program (33) as an optimization problem over the space R corresponding to positive definite matrices. Doing so makes the program (33) consistent with the framework laid out earlier in our paper. In this case, the oracle estimator is defined by With this setup, we have the following guarantee: Corollary 4. Given a sample size n ≥ c 0 d 2 log p, suppose the x i 's are drawn from a sub-Gaussian distribution, and the regularizer ρ λ is (µ, γ)-amenable. Also suppose that Then with probability at least 1 − c 1 exp(−c 2 log p), the program (33) with κ = 2 µ has a unique stationary point Θ given by the oracle estimator Θ O , and Moreover, as shown in the proof in Appendix D.4, the support containment condition and elementwise bound (35) also imply bounds on the Frobenius and spectral norms of the error, namely We reiterate that Corollary 4 does not involve any restrictive incoherence assumptions on the matrix Θ * . As remarked in past work [22,32], such incoherence conditions for the graphical Lasso are very restrictive-much more so than the corresponding incoherence conditions for graph recovery using neighborhood regression [23,45,40]. Thus, in this setting, our corollary illustrates another distinct advantage of nonconvex regularization.

Simulations
In this section, we report the results of various simulations that we ran in order to verify our theoretical results.

Optimization algorithm
We begin by describing the algorithm we use to optimize the program (1). We rewrite the program as β ∈ arg min and apply the composite gradient descent algorithm due to Nesterov [26]. The updates of the composite gradient procedure are given by β t+1 ∈ arg min where 1 η is the stepsize. In the particular simulations of this section, we take Ω = R p . Then the iterates (38) have the convenient closed-form expression The following proposition guarantees the computational efficiency of the general composite gradient descent algorithm (38). For ease of analysis, we assume that β * 2 ≤ 1, and the Taylor error T (β 1 , β 2 ) := L n (β 1 )−L n (β 2 )− ∇L n (β 2 ), β 1 −β 2 satisfies the following restricted strong convexity condition, for all β 2 ∈ B 2 (3) ∩ B 1 (R): as well as the restricted smoothness condition We also assume for simplicity that q λ is convex, as is the case for all the regularizers studied in this paper. In the following statement, we let β be the unique global optimum of the program (37). Also denote φ(β) := L n (β) + ρ λ (β). Proposition 1. Suppose L n satisfies the RSC (40) and RSM (41) conditions, and assume the regularizer ρ λ is µ-amenable with µ 2 < α := min{α 1 , α 2 }, and q λ is convex. Let the scalars (R, λ) be chosen to satisfy the bounds R log p n ≤ c and λ ≥ 4 max ∇L n (β * ) ∞ , τ log p n .
Since we are interested in ℓ ∞ -error bounds, we state a simple corollary to Proposition 1 that ensures convergence of the iterates in ℓ ∞ -norm, up to accuracy O log p n , assuming statistical consistency of the global optimum and the scaling n k 2 log p: Corollary 5. Suppose, in addition to the assumptions of the previous proposition, that β − β * 2 ≤ c k log p n , and the sample size is lower-bounded as n ≥ c 0 k 2 log p. Then the iterates of the composite gradient descent algorithm (38) satisfy the ℓ ∞ -bound The proof of Corollary 5 is a simple consequence of Proposition 1 and is supplied in Appendix E.2. Note that the bound β − β * 2 ≤ c k log p n holds, with high probability, as a consequence of Lemma 9 in Appendix F.2. Corollary 5 has a natural consequence for support recovery: Suppose the estimate β satisfies an ℓ ∞ -bound of the form β − β * ∞ ≤ c 3 log p n , as guaranteed by Theorem 2. Combining the ℓ ∞ -bound (42) on iterate β t with the triangle inequality, we are then guaranteed that so the composite gradient descent algorithm converges to a vector with the correct support,

Classes of matrices
Next, we describe two classes of matrices to be used in our simulations. The first class consists of matrices that do not satisfy the incoherence conditions, although the maximum and minimum eigenvalues are bounded by constants. We define Hence, M 1 (θ) is a matrix with 1's on the diagonal, θ's in the first k positions of the (k + 1) st row and column, and 0's everywhere else. The following lemma, proved in Appendix E.3, provides the incoherence parameter and eigenvalue bounds for M 1 (θ) as a function of θ: Lemma 1. With the shorthand notation Γ = M 1 (θ) and S = {1, . . . , k}, the incoherence parameter is given by Γ S c S Γ −1 SS ∞ = kθ, and the minimum and maximum eigenvalues are given by λ min (Γ) = 1 − θ √ k and λ max (Γ) = 1 + θ √ k.
In particular, if θ ∈ 1 k , 1 √ k , Lemma 1 ensures that M 1 (θ) has bounded eigenvalues but does not satisfy the incoherence condition.

Experimental results
We ran experiments with the loss function coming from (a) ordinary least squares linear regression, (b) least squares linear regression with corrupted covariates, and (c) logistic regression. For all our simulations, we used the regularization parameters R = 1.1 β * 1 and λ = log p n , and we set the SCAD and MCP parameters to be a = 2.5 and b = 1.5, respectively. Note that although the covariates in our simulations for logistic regression do not satisfy the boundedness Assumption 1(i) imposed in our corollary, the generated plots still agree qualitatively with our predicted theoretical results.
In our first set of simulations, we show that using the SCAD or MCP regularizer in situations where the design matrix does not satisfy incoherence conditions still results in an estimator that is variable selection consistent. We generated i.i.d. covariates x i ∼ N (0, Σ x ), where Σ x = M 1 (θ) was obtained from the family of non-incoherent matrices (43), with θ = 2.5 k . We chose k ≈ √ p and β * = 1 √ k , · · · , 1 √ k , 0, · · · , 0 , the unit vector with the first k components equal to 1 √ k , and generated response variables according to the linear model y . In addition, we generated corrupted covariates z i = x i + w i , where w i ∼ N (0, (0.2) 2 ), and w i ⊥ ⊥ x i . We then ran the composite gradient descent algorithm with updates given by equation (39), where the loss function is given by equation (9), and ( Γ, γ) are defined as in equation (26). Figure 1 shows the results of our simulations for the problem sizes p = 128, 256, and 256. In panel (a), we see that the probability of correct support recovery transitions sharply from 0 to 1 as the sample size increases and ρ λ is the SCAD or MCP regularizer. In contrast, the probability of recovering the correct support remains at 0 when ρ λ is the ℓ 1 -penalty or LSP-by the structure of Σ x , regularization with the ℓ 1 -penalty or LSP results in an estimator β that puts nonzero weight on the (k + 1) st coordinate, as well. Note that we have rescaled the horizontal axis according to n k log p in order to match the scaling prescribed by our theory; the three sets of curves for each regularizer roughly align, as predicted by Theorem 1. Panel (b) confirms that the ℓ ∞ -error β − β * ∞ decreases to 0 when using the SCAD and MCP, as predicted by Theorem 2. Finally, we plot the ℓ 2 -error β − β * 2 of the SCAD and MCP regularizers alongside the ℓ 2 -error for the ℓ 1 -penalty and LSP in panel (c). Although the ℓ 2 -error is noticeably smaller for SCAD and MCP than for the ℓ 1 -penalty and LSP, as noted by previous authors [11,4,21], all four regularizers are nonetheless consistent in ℓ 2 -error, since a lower-eigenvalue bound on the covariance matrix of the design is sufficient for ℓ 2 -consistency [20]. For our choice of regularization parameters, where the same value of λ is shared between the ℓ 1 -penalty and LSP, the two sets of curves for the ℓ 1 -penalty and LSP nearly agree; as shown in Candes et al. [6], the relative improvement of the LSP in comparison to the ℓ 1 -penalty may vary widely depending on the regularization parameter.
In our second set of simulations, we explore the uniqueness of stationary points of the composite objectives. We focus on settings where the loss function comes from either linear regression with ordinary least squares, or logistic regression. Our theory guarantees that stationary points are unique when µ < 2α 1 , but when 2α 1 ≤ µ, multiple stationary points may emerge. In fact, when 2α 1 ≤ µ, convergence of the composite gradient descent algorithm and consistent support recovery are no longer guaranteed. In practice, we observe that multiple initializations of the composite gradient descent algorithm still appear to converge to a single  stationary point with the correct support, when µ is slightly larger than 2α 1 ; however, when that condition is violated more severely, the composite gradient descent algorithm indeed terminates at several distinct stationary points. Figure 2 shows the result of multiple runs of the composite gradient descent algorithm with different regularizers in the cleanly-observed linear regression setting. We generated observations x i ∼ N (0, Σ x ), with Σ x = M 2 (θ) coming from the family of spiked identity models (44), for θ = 0.7 and 0.8, and independent additive noise, ǫ i ∼ N (0, (0.1) 2 ). We set the problem dimensions to be p = 128, k ≈ √ p, and n ≈ 20k log p, and generated β * to have k nonzero values ± 1 √ k with equal probability for each sign. When using the SCAD or MCP regularizers (panels (b) and (d)), distinct stationary points emerge and the recovered support is incorrect, since 2α 1 < µ. In contrast, the ℓ 1penalty and LSP still continue to produce unique stationary points with the correct support (panels (a) and (b)). Observe from the plots in Figure 2 that the error β t − β * 2 decreases at a rate that is linear on a log scale, as predicted by Theorem 3 of Loh and Wainwright [20], until it reaches the threshold of statistical accuracy. Further note the significant increase in overall precision from using SCAD or MCP, as seen by comparing the vertical axes in panels (a) & (b) and panels (c) & (d).
Finally, we present a third set of simulations, analogous to the second, with the OLS loss function replaced by the maximum likelihood loss function for logistic regression: We generated x i ∼ N (0, σ 2 x I), with σ x ∈ {1, 3}, and set the problem dimensions to be p = 128, k ≈ √ p, and n ≈ 10k 3 log p. We generated β * to have k nonzero values ± 1 √ k with equal probability for each sign, and generated response variables y i ∈ {0, 1} according to . Figure 3 shows the results of our simulations. Panels (a)-(d) plot the log ℓ 2 -error as a function of iteration number, when σ x = 1. Note that in this case, an empirical evaluation shows that λ min (∇ 2 L(β * )) ≈ 0.14, so we expect α 1 ≈ 0.14 and µ ≮ 2α 1 . As in the plots of Figure 2, multiple stationary points emerge in panels (c) and (d) when ρ λ is the SCAD or MCP regularizer; in contrast, we see from panels (a) and (b) that all 15 runs of composite gradient descent converge to the same stationary point when ρ λ is the ℓ 1 -penalty or LSP. In panels (e) and (f), we repeat the simulations with σ x = 3. In this case, λ min (∇ 2 L(β * )) ≈ 0.25, and we see from our plots that although the condition µ < 2α 1 is still violated, the larger value of α 1 is enough to make the stationary points under SCAD or MCP regularization unique. We may again observe the geometric rate of convergence of the ℓ 2 -error β t − β * 2 in each plot, up to a certain small threshold. The improved performance from using the SCAD and MCP regularizers may be observed empirically by comparing the vertical axes in the panels of Figure 3.   x I), with σ x = 1 in plots (a)-(d), and σ x = 3 in plots (e)-(f). In panels (c) and (d), the composite gradient descent algorithm settles into multiple distinct stationary points, which exist because 2α 1 < µ for the SCAD and MCP. However, when the covariates have a larger covariance, the SCAD and MCP regularizers produce unique stationary points, as observed in panels (e) and (f).

Discussion
We have developed an extended framework for analyzing a variety of nonconvex problems via the primal-dual witness proof technique. Our results apply to composite optimization programs where both the loss and regularizer function are allowed to be nonconvex, and our analysis significantly generalizes the machinery previously established to study convex objective functions. As a consequence, we have provided a powerful reason for using nonconvex regularizers such as the SCAD and MCP rather than the convex ℓ 1 -penalty: In addition to being consistent in ℓ 2 -error, the nonconvex regularizers actually produce an overall estimator that is consistent for support recovery when the design matrix is non-incoherent and the usual ℓ 1 -regularized program fails in recovering the correct support. We have also established a similar strong result for the graphical Lasso objective function with nonconvex regularizers, which eliminates the need for complicated incoherence conditions on the inverse covariance matrix.
Future directions of research include devising theoretical guarantees when the condition µ < 2α 1 is only mildly violated, since the condition does not appear to be strictly necessary based on our simulations, and establishing a rigorous justification for why the SCAD and MCP regularizers perform appreciably better than the ℓ 1 -penalty even in terms of ℓ 2 -error, in situations where the assumptions are not strong enough for an oracle result to apply. It is also an open question as to how generally the restricted strong convexity condition may hold for various other nonconvex loss functions of interest, or whether a local RSC condition is sufficient to guarantee good behavior with proper initializations of a gradient descent algorithm. Finally, it would be useful to be able to compute the RSC constants (α 1 , α 2 ) empirically from data, so as to assign a nonconvex regularizer with the proper amount of curvature.

Acknowledgments
The work of PL was partly supported from a Hertz Foundation Fellowship, an NSF Graduate Research Fellowship and a PD Fellowship while studying at Berkeley. MJW and PL were also partially supported by grants NSF grant DMS-1107000, NSF grant CIF-31712-23800, and Air Force Office of Scientific Research Grant AFOSR-FA9550-14-1-0016.

A Proof of Theorem 1
In this Appendix, we provide the proof of Theorem 1. We begin with the main body of the argument, with the proofs of some more technical lemmas deferred to later subsections.

A.1 Main part of proof
We follow the outline of the primal-dual witness construction described in Section 2.4. For step (i) of the construction, we use Lemma 9 in Appendix F, where we simply replace p by k and L n by (L n ) S , which is the function L n restricted to R S . It follows that as long as Here, the final inequality follows by the lower bound in inequality (15b). We conclude that β S must be in the interior of the feasible region. Moving to step (ii) of the PDW construction, we define the shifted objective function L n (β) := L n (β) − q λ (β). Since β S is an interior point, it must be a zero-subgradient point for the restricted program (13), so ∇(L n ) S ( β S ) + λ z S = 0, where z S ∈ ∂ β S 1 is the dual vector.
By the chain rule, this implies that ∇L n ( β) S + λ z S = 0, where β := ( β S , 0 S c ). Accordingly, we may define the subvector z S c ∈ R S c such that where z := ( z S , z S c ) is the extended subgradient. Under the assumption (16), this completes step (ii) of the construction. For step (iii), we first establish that β is a local minimum for the program (1) by verifying the sufficient conditions of Lemma 10 in Appendix F, with f = L n , g = q λ , and (x * , v * , w * , µ * ) = ( β, z, z, 0). Note that Lemma 8(b) from Appendix F.1 ensures the concavity and differentiability of g(x) − µ 2 x 2 2 . Since µ * = 0, condition (100a) is trivially satisfied. Furthermore, condition (100b) holds by equation (45). Hence, it remains to verify the condition (100c).
We first show that G * ⊆ R S . Supposing the contrary, consider a vector ν ∈ G * such that supp(ν) S. Fixing some index j ∈ S c such that ν j = 0, by the definition of G * , we have However, if z denotes the vector z with entry j replaced by sign(s j ) ∈ {−1, 1}, we clearly still have z ∈ ∂ β 1 , but ν T (∇L n ( β) + λ z) > ν T (∇L n ( β) + λ z) = 0, where the strict inequality holds because z S c ∞ < 1, by our assumption. We have thus obtained a contradiction to equation (46); consequently, our initial assumption was false, and we may conclude that G * ⊆ R S .
In particular, since G * ⊆ S and supp( β) ⊆ S, Lemma 2 immediately implies condition (100c) of Lemma 10. We conclude that β is indeed a local minimum of the program (1).
The following lemma, proved in Appendix A.3, show that all stationary points of the program (1) are supported on S: Turning to the uniqueness assertion, note that since all stationary points are supported in S, any stationary point β of the program (1) must satisfy β = ( β S , 0 S c ), where β S is a stationary point of the restricted program (13). By Lemma 2, the restricted program is strictly convex. Hence, the vector β S , and consequently also β, is unique.

A.2 Proof of Lemma 2
We begin by establishing the bound ∇ 2 L n (β) SS α 1 − τ 1 k log p n I, for all β ∈ R p . Equivalently, When this lower bound holds, we are guaranteed that so L n (β) − µ 2 β 2 2 is strictly convex on R S under the prescribed sample size. In order to prove the bound (47), consider a fixed v ∈ R p such that supp(v) ⊆ S and v 2 = 1. For this vector, we have Furthermore, by the RSC assumption (8) and the fact that v 2 = 1, we have Combining this bound with equations (48) and (49) then gives the desired inequality (47).
Finally, we note the decomposition L n (β) = L n (β) − µ 2 β 2 2 + µ 2 β 2 2 − q λ (β) , showing that L n is the sum of a strictly convex and convex function over R S . Hence, L n | S is strictly convex, as claimed. The strict convexity of L n (β) + λ β 1 over R S follows immediately.
Now, applying the RSC condition (8b), we have which implies that By inequality (50), we also have where z ∈ ∂ β 1 . From the zero-subgradient condition (14), we have ∇L n ( β)+λ z, β− β = 0. Combining with inequality (53) then yields Rearranging, we have where the second inequality comes from the fact that z, β ≤ z ∞ · β 1 ≤ β 1 , and the third inequality comes from the bound (52). Finally, we need a lemma showing that ν lies in a cone set: Proof. From inequality (54) together with the RSC bound (52), we have Note that since supp( β) ⊆ S, we have Combining inequalities (56a), (56b), and (56c) then yields Under the assumption on λ, we have τ 1 log p n Putting together the pieces, we then have which establishes the claim.
By the zero-subgradient condition (14) and the invertibility of ∇ 2 L n ( β) SS , we then have which is inequality (18).
Lemma 5 implies that ∇q λ ( β S ) = λ z S . Hence, the zero-subgradient condition (14) reduces to ∇L n ( β S ) S = 0. Since L n is strictly convex on R S by Lemma 2, this zero-gradient condition implies that β S is the unique global minimum of (L n ) S , so we have the equalities β S = β O S and β = β O , as claimed. Finally, inequality (58) simplifies to inequality (19b), using the equality ∇q λ ( β S ) = λ z S .

C Establishing strict dual feasibility
In this Appendix, we derive conditions that allow us to establish the strict dual feasibility conditions required in order to apply Theorem 1. We use these derivations to prove Corollaries 1-3.
Defining Q := 1 0 ∇ 2 L n β * + t( β − β * ) dt, we then have In block form, this means By simple algebraic manipulations, we then have Note that by the selection property (vi), we have ∇q λ ( β S c ) = ∇q λ (0 S c ) = 0 S c , so the first term in equation (60) vanishes, giving If the unbiasedness property (vii) also holds, we have the following result: Under the conditions of Theorem 1, suppose ρ λ is (µ, γ)-amenable. Also suppose inequality (19a) holds. Then strict dual feasibility holds provided Proof. Combining the result of Lemma 5 in Appendix B with equation (61), we see that as claimed.
However, for regularizers that do not satisfy the unbiasedness property (vii), such as the LSP and the usual ℓ 1 -norm, we impose slightly stronger conditions to ensure strict dual feasibility. The notion of an incoherence condition is taken from the Lasso literature [45,40,24]: Based on this condition, we have the following result: Proposition 3. Under the conditions of Theorem 1, suppose in addition that Assumption 2 holds. Then strict dual feasibility holds provided Proof. Lemma 8(a) implies that λ z S − ∇q λ ( β S ) ∞ = ∇ρ λ ( β S ) ∞ ≤ λ. Hence, equation (61), together with the triangle inequality, implies where we have used inequalities (63a) and (63b) and the incoherence condition (Assumption 2).
Comparing Propositions 2 and 3 reveals a theoretical benefit of using unbiased penalties such as SCAD or MCP over LSP or the ℓ 1 -norm: Whereas the bound (62b) is essentially the same as the bound (63b), Proposition 3 imposes an additional incoherence condition on the covariates.

D Proofs of corollaries in Section 3
In this section, we provide detailed proofs of the corollaries to Theorems 1 and 2 appearing in Section 3.

D.1.1 Proof of part (a)
We now use the machinery developed in Appendix C. We use a slightly more direct analysis than the result stated in Proposition 3, which allows us to produce tighter bounds. Denoting ( Γ, γ) = X T X n , X T y n , we have so equation (61) takes the form using the fact that λ z s − ∇q λ ( β S ) ∞ = ∇ρ λ ( β S ) ∞ ≤ λ, by Lemma 8. We now write where Π := I − X S (X T S X S ) −1 X T S is an orthogonal projection matrix. Note that for t > 0, we have , since ǫ is a sub-Gaussian vector with parameter σ 2 ǫ , and Π is a projection matrix. Taking another expectation, we have so e T j X T S c Πǫ n is sub-Gaussian with parameter cσ 2 ǫ σ 2 x n . Using sub-Gaussian tail bounds and a union bound, we conclude from equation (66) that with probability at least 1−c exp(−c ′ log p), provided n k log p. Hence, strict dual feasibility holds, w.h.p., provided λ > C 1−η log p n , and the variable selection consistency property follows by Theorem 1.
Turning to ℓ ∞ -bounds, we have For j ∈ S, we may write e T j X T S X S n −1 X T S ǫ n with probability at least 1 − c 1 exp(−c 2 k). Conditioned on X, the variables v T j ǫ are sub-Gaussian with parameter σ 2 ǫ v j 2

. Denoting
A := |v T j ǫ| ≤ σ ǫ v j 2 log p, ∀j ∈ S and applying a union bound, we then have as well. We conclude that with probability at least 1 − c ′′ 1 exp(−c ′′ 2 min{k, log p}). Combining this with the assumption (21), we conclude by part (a) of Theorem 2 that the desired ℓ ∞ -bound holds.

D.1.2 Proof of part (b)
The proof of this corollary is nearly identical to the proof of part (a), except that equation (65) does not have the extra term η; by the assumption of (µ, γ)-amenability, Lemma 5 implies that λ z S − ∇q λ ( β S ) = 0. Hence, equation (61) takes the form Observe that the matrix concentration arguments from the proof of part (a) also imply that The remainder of the proof follows by part (b) of Theorem 2.

D.2 Proof of Corollary 2
The expressions for ∇L n and ∇ 2 L n are as in equation (64), with ( Γ, γ) defined in equation (26). Furthermore, as remarked in the proof of Corollary 1, Corollary 1 of our previous work [20] implies that the RSC condition (8) holds w.h.p., when α 1 = α 2 = 1 2 λ min (Σ x ) and τ 1 , τ 2 ≍ 1. We may verify using standard techniques that when the x i 's and ǫ i 's are sub-Gaussian and n k log p, we have the bound ∇L n (β * ) ∞ log p n , with probability at least 1 − c 1 exp(−c 2 log p) [19]. Hence, if λ log p n , the lower bound in inequality (15a) is satisfied w.h.p. We may then check that for the choice of regularization parameters λ ≍ log p n and R ≍ 1 λ , and under the scaling n k log p, both bounds in condition (15a) hold. We now use Proposition 3 to verify condition (a) in Theorem 1, establishing strict dual feasibility (16). From the bound (63b), and with the scaling λ ≍ log p n , it suffices to show that w.h.p. To this end, we have the following proposition, which we state in some generality to indicate its applicability to other types of corrupted linear models [19]: Also suppose γ satisfies the deviation bound and suppose |||Γ||| 2 and Γ −1 2 are uniformly bounded. Under the scaling n k · max{log p, k} and λ ≍ log p n , inequality (67) holds, with probability at least 1 − c exp(c ′ min{k, log p}). Proof. First note that inequality (68) implies the following deviation bounds: Inequality (70a) follows from from a discretization argument over the ℓ 2 -ball in R S , and inequality (70b) is similar. Inequality (70c) follows by taking w = Ae j and a union bound over 1 ≤ j ≤ p. Furthermore, by Lemma 11, inequality (70a) implies that as well.
Note that by inequalities (70c) and (69), we have with probability at least 1 − c log(−c ′ log p). Furthermore, the triangle inequality gives and the first term is bounded above by C log p n , by inequality (72). For the second term, we write where we have again used inequality (72) in the second inequality. Furthermore, The terms in inequality (75) are bounded as Combining the bounds with inequalities (70b) and (71), we conclude from inequality (75) that with probability at least 1 − c exp(−c ′ min{k, log p}), under the scaling n max{k, log p}. Plugging inequality (76) into inequality (74) and combining with inequality (72), we conclude that with probability at least 1 − c exp(−c ′ min{k, log p}). The desired result then follows under the scaling n k · max{k, log p}.
It is easy to check that the bounds (68) and (69) hold when ( Γ, γ) are defined by equation (26) and X, W , and ǫ are sub-Gaussian. Hence, strict dual feasibility holds by Propositions 3 and 4, and the PDW technique succeeds.
Turning to ℓ ∞ -bounds, note that We have the following result: Proposition 5. Under the same conditions as Proposition 4, we have with probability at least 1 − c exp(−c ′ min{k, log p}).
Proof. We decompose where the last inequality holds with probability at least 1 − c exp(−c ′ min k, log p), and we have used the bounds (71) and (72) to establish the last inequality. Inequality (77) then follows under the scaling n ≥ k 2 .

D.3 Proof of Corollary 3
Note that Corollary 2 of our previous work [20] establishes the RSC condition (8) when ψ ′′ is bounded and the x i 's are again sub-Gaussian. We again proceed via the framework and terminology of Appendix C, particularly Proposition 2. Taking derivatives of the loss function (11), we have To verify inequality (62a), we may use straightforward sub-Gaussian concentration techniques (cf. the proof of Corollary 2 in Loh and Wainwright [20]). For inequality (62b), note that by the mean value theorem, where β t lies on the line segment between β * and β * + t( β − β * ).
where x i,S denotes the vector x i restricted to S. Inequality (78) follows by two applications of the Cauchy-Schwarz inequality and the fact that supp( β) ⊆ S. Furthermore, by Lemma 9 in Appendix F, we have with probability at least 1 − c exp(−c ′ log p). By Assumption 1(i), we also have In particular, taking a supremum over unit vectors v, w ∈ R S , we have with probability at least 1 − c exp(−c ′ k), where the second inequality holds by a standard spectral norm bound on the sample covariance matrix.
average of a product of sub-Gaussian random variables ψ ′′ (x T i β * )v T x i and w T x i , since the x i 's are sub-Gaussian and ψ ′′ is uniformly bounded by assumption. Defining L(β * ) := E [L n (β * )], a standard discretization argument of k-dimensional unit sphere yields with probability at least 1 − c exp(−c ′ k). Hence, by inequalities (79) and (80), we have . Applying Lemma 11 in Appendix F yields as well. A similar argument shows that with probability at least 1− c exp(−c ′ min{k, log p}). Putting together the pieces, we find that with high probability. Returning to the expression (62b), we have the bound Since the eigenspectrum of ∇ 2 L(β * ) is bounded by assumption, standard techniques (cf. the proofs of Corollary 2 in Loh and Wainwright [20] and Lemma 6 in Negahban et al. [25]) guarantee that T 1 ≤ c log p n , with probability at least 1 − c 1 exp(−c 2 log p). Turning to the second term, we have Now, and by the triangle inequality, where By the bounds (81) Hence, by Proposition 2, strict dual feasibility holds under the scaling n k 3 log p.
Turning to ℓ ∞ -bounds, we have using inequality (81) and the same sub-Gaussian concentration result used to bound the term T 1 above. Hence, under the scaling n k 3 log p. Finally, note that inequality (81) implies Theorem 2 then implies the desired result.

D.4 Proof of Corollary 4
In this section, we provide the proof of Corollary 4. Corollary 3 of our previous work [20] establishes the RSC condition (8), with α 1 and α 2 scaling as |||Θ * ||| 2 and τ 1 = 0. Our proof strategy deviates mildly from the framework described in Section 2.4, in that we provide a slightly different way of constructing a matrix Θ S such that supp( Θ S ) ⊆ S and Θ S is a zerosubgradient point of the restricted program (13). However, subject to this minor adjustment, the remainder of the primal-dual witness technique proceeds as before.
We begin with a simple lemma establishing the convexity of the program (33): Lemma 6. Suppose ρ λ is µ-amenable. Then for the choice of parameter κ = 2 µ , the objective function in the program (33) is strictly convex over the constraint set.
In particular, if there exists a zero-subgradient point of the composite objective function L n (Θ) + ρ λ (Θ) within the feasible set, it must be the unique global minimum. Our strategy is to construct such a zero-subgradient point. Let S := {(j, k) | j ≤ k, Θ * jk = 0} denote the support of Θ * , where we remove redundant elements. We define the map F : R |S| → R |S| according to where ∆ ∈ R p×p is the symmetric matrix agreeing with ∆ S on S and having 0's elsewhere, and Γ * := Θ * −1 ⊗ Θ * −1 = Σ * ⊗ Σ * . We analyze the behavior of F over the ℓ ∞ -ball B ∞ (r) of radius r to be specified later. In particular, note that for ∆ S ∈ B ∞ (r), we have since ∆ has at most d nonzero entries per row. Hence, when dr < λ min (Θ * ), we are guaranteed that the matrix (Θ * +∆) is invertible, making F a continuous map. We show that F (B ∞ (r)) ⊆ B ∞ (r), so by Brouwer's fixed point theorem [28], the function must have a fixed point, which we denote by ∆ * S ∈ B ∞ (r). Defining the constants κ Γ := (Γ * SS ) −1 ∞ and κ Σ := |||Σ * ||| ∞ , this insight is summarized in the following lemma: and suppose the sample size satisfies n κ 2 Γ log p. Then with probability at least 1−c exp(−c ′ log p), there exists Θ ∈ R p×p such that Furthermore, if ρ λ is (µ, γ)-amenable and min (j,k)∈S |Θ * jk | ≥ γλ + r, we have ρ ′ λ ( Θ jk ) = 0, for all (j, k) ∈ S.
Proof. We first establish that F (B ∞ (r)) ⊆ B ∞ (r). Consider ∆ S ∈ B ∞ (r). We have For the first term, we have w.h.p., by the sub-Gaussian assumption on the x i 's. Furthermore, by the matrix expansion we have By the triangle inequality, we then have By Hölder's inequality, we have Using inequality (87) and plugging back into inequality (93), we then have where the last inequality follows from our assumption (88). Combining inequalities (90), (91), and (94), and using the assumption that 2κ Γ κ 3 Σ dr 2 ≤ r 2 , we find that F (vec(∆ S )) ∞ ≤ r, as desired. Applying Brouwer's fixed point theorem then yields the required fixed point ∆ * S . Defining Θ := Θ * + ∆ * , the third equality in line (89) follows. The operator norm bound in inequality (89) follows from the same argument as in inequality (87). Finally, by the triangle inequality and the assumption that ρ λ is (µ, γ)-amenable, we have ρ ′ λ ( Θ jk ) = ρ ′ λ (Θ * jk ) = 0.
Equipped with Lemma 7, we now prove that Θ is the (unique) global optimum of the program (33). Note that it suffices to show that Θ is a zero-subgradient point of the objective function in the program (33); since the objective is strictly convex by Lemma 6, Θ must be the global minimum. The same manipulations used in Appendix C and Proposition 2 reveal that Θ is a zero-subgradient point satisfying dual feasibility, if inequalities (62a) and (62b) hold. In fact, we may take δ = 0, since the convexity of the problem immediately implies uniqueness of a stationary point if it exists. We treat ∇ 2 L n (Θ) as a p 2 × p 2 matrix operating on the p-dimensional vector vec(Θ). For the graphical Lasso, we have ∇L n (Θ * ) = Σ− Σ * , so we need to show that Σ − Σ * max ≤ λ 2 . It is straightforward to see that by the sub-Gaussianity of the x i 's, setting λ = c log p n for a suitably large constant is sufficient to satisfy inequality (62a). Turning to inequality (62b), note that Furthermore, for t ∈ [0, 1], we have using Lemma 7. Lemma 11 in Appendix F then implies that Hence, by Lemma 13 in Appendix F, we have implying the bounds and Applying Lemma 11 in Appendix F yields Next we define the pair (R, Ξ) according to Note that we have R ≤ |||Ξ| || ∞ · (∇L n (Θ * )) S max |||Ξ||| ∞ · log p n , and moreover, Combined with our earlier upper bound on R, we conclude that R log p n , under the scaling n d 2 log p.
Next, observe that Hence, the primal-dual witness technique succeeds. Note that by inequality (95), we also have By Theorem 2, we conclude that Θ is the unique global minimum of the program (33) with the desired properties. Finally, let us prove the claimed bounds on the Frobenius and spectral norms. Note that The bounds then follow from our earlier bound on the elementwise ℓ ∞ -norm.

E Proofs for Section 4
In this section, we provide details of proofs for the results in Section 4.

E.2 Proof of Corollary 5
We set δ ≍ k log p n · β − β * 2 in Proposition 1. Then where the last inequality follows by the assumption of statistical consistency for β. Under the scaling n k 2 log p, the desired result follows.

E.3 Proof of Lemma 1
The computation of the incoherence parameter is straightforward. For the spectral properties, note that v T Γv = 1 + 2θv k+1 k j=1 v j for any vector v ∈ R p . Consequently, we have We may write where we have used the fact that w 1 ≤ √ k w 2 . It is easy to see that the final expression is maximized when α = 1 √ 2 , so λ max (Γ) ≤ θ √ k. Equality is achieved for the vector (v 1 , . . . , v k , v k+1 ) =

F Some useful auxiliary results
Finally, we provide some useful auxiliary results, which we employ in the proofs of our main theorems.

F.1 Properties of amenable regularizers
The following lemma is useful in various parts of our analysis. Part (a) is based on Lemma 4 of Loh and Wainwright [20].

F.3 Sufficient conditions for local minima
The following lemma is a minor extension of results from Fletcher and Watson [14]. It applies to functions f ∈ C 2 and g ∈ C 1 , such that g(x) − κ 2 x 2 2 is concave, for some κ ≥ 0.
Proof. The proof of this lemma essentially follows the proof of Theorem 3 of Fletcher and Watson [14], except it allows for a composite function h = f − g in the objective that is not in C 2 . Nonetheless, we include a full proof for clarty and completeness. Suppose for the sake of contradiction that x * is not an isolated local minimum. Then there exists a sequence of feasible points {x (k) } → x * with φ(x (k) ) ≤ φ(x * ), where φ(x) := h(x) + λ x 1 .
Let s (k) := x (k) −x * x (k) −x 2 ; then {s (k) } is a set of feasible directions. Since {s (k) } ⊆ B 2 (1), the set must possess a point of accumulation s ∈ B 2 (1), and we may extract a convergent subsequence. Relabeling the points as necessary, we assume that {s (k) } → s. We show that s ∈ G * .
Since the feasible region is closed, s is also a feasible direction at x * . In particular, if x * 1 = R, we must have sup Together with equation (100a), this implies that Note that by equation (100b), we have s T (∇h(x * ) + λv * ) = −µ * s T w * ≥ 0, where the inequality follows from the fact that if µ * = 0, we have x * 1 = R, by equation (100a), Hence, s T w * ≤ 0, by inequality (101).
By the definition of the subgradient, we have

F.4 Some matrix-theoretic lemmas
Here, we collect some useful lemmas on matrices and Kronecker products.
Lemma 11. Let A, B ∈ R p×p be invertible. For any matrix norm |||·|||, we have In particular, if A −1 |||A − B||| ≤ 1/2, then Proof. We use the matrix expansion (92), where ∆ = B − A. By the triangle inequality and multiplicativity of the matrix norm, we then have as claimed.
By symmetry, the same bound holds with the roles of A and B reversed, from which the claim follows.