Regularized M-estimators with nonconvexity: Statistical and algorithmic theory for local optima

Introduction

Although recent years have brought about a flurry of work on optimization of convex functions, optimizing nonconvex functions is in general computationally intractable (Nesterov and Nemirovskii 1987; Vavasis 1995). Nonconvex functions may possess local optima that are not global optima, and iterative methods such as gradient or coordinate descent may terminate undesirably in local optima. Unfortunately, standard statistical results for nonconvex $M$-estimators often only provide guarantees for global optima. This leads to a significant gap between theory and practice, since computing global optima—or even near-global optima—in an efficient manner may be extremely difficult in practice. Nonetheless, empirical studies have shown that local optima of various nonconvex $M$-estimators arising in statistical problems appear to be well-behaved (e.g., Breheny and Huang 2011). This type of observation is the starting point of our work.

A key insight is that nonconvex functions occurring in statistics are not constructed adversarially, so that “good behavior” might be expected in practice. Our recent work (Loh and Wainwright 2012) confirmed this intuition for one specific case: a modified version of the Lasso applicable to errors-in-variables regression. Although the Hessian of the modified objective has many negative eigenvalues in the high-dimensional setting, the objective function resembles a strongly convex function when restricted to a cone set that includes the stationary points of the objective. This allows us to establish bounds on the statistical and optimization error.

Our current paper is framed in a more general setting, and we focus on various $M$-estimators coupled with (nonconvex) regularizers of interest. On the statistical side, we establish bounds on the distance between any local optimum of the empirical objective and the unique minimizer of the population risk. Although the nonconvex functions may possess multiple local optima (as demonstrated in simulations), our theoretical results show that all local optima are essentially as good as a global optima from a statistical perspective. The results presented here subsume our previous work (Loh and Wainwright 2012), and our present proof techniques are much more direct.

Our theory also sheds new light on a recent line of work involving the nonconvex SCAD and MCP regularizers (Fan and Li 2001; Breheny and Huang 2011; Zhang 2010; Zhang and Zhang 2012). Various methods previously proposed for nonconvex optimization include local quadratic approximation (LQA) (Fan and Li 2001), minorization-maximization (MM) (Hunter and Li 2005), local linear approximation (LLA) (Zou and Li 2008), and coordinate descent (Breheny and Huang 2011; Mazumder et al. 2011). However, these methods may terminate in local optima, which were not previously known to be well-behaved. In a recent paper, Zhang and Zhang 2012 provided statistical guarantees for global optima of least-squares linear regression with nonconvex penalties and showed that gradient descent starting from a Lasso solution would terminate in specific local minima. Fan et al. 2014 also showed that if the LLA algorithm is initialized at a Lasso optimum satisfying certain properties, the two-stage procedure produces an oracle solution for various nonconvex penalties. Finally, Chen and Gu 2014 showed that specific local optima of nonconvex regularized least-squares problems are stable, so optimization algorithms initialized sufficiently close by will converge to the same optima. See the survey paper (Zhang and Zhang 2012) for a more complete overview of related work.

In contrast, our paper is the first to establish appropriate regularity conditions under which all stationary points (including both local and global optima) lie within a small ball of the population-level minimum. Thus, standard first-order methods such as projected and composite gradient descent (Nesterov 2007) will converge to stationary points that lie within statistical error of the truth, eliminating the need for specially designed optimization algorithms that converge to specific local optima. Our work provides an important contribution to the growing literature on the tradeoff between statistical accuracy and optimization efficiency in high-dimensional problems, establishing that certain types of nonconvex $M$-estimators arising in statistical problems possess stationary points that both enjoy strong statistical guarantees and may be located efficiently. For a higher-level description of contemporary problems involving statistical and optimization tradeoffs, see Wainwright 2014 and the references cited therein.

Panel (b) exhibits the same behavior for a problem in which both the cost function (a corrected form of least-squares suitable for missing data, described in Loh and Wainwright 2013a) and the regularizer (the MCP function, described in Zhang 2010) are nonconvex. Nonetheless, as guaranteed by our theory, we still see the same qualitative behavior of the statistical and optimization error. Moreover, our theory also predicts the geometric convergence rates that are apparent in these plots. More precisely, under the same sufficient conditions for statistical consistency, we show that a modified form of composite gradient descent only requires $\log(1/\epsilon_{\mbox{\tiny{stat}}})$ steps to achieve a solution that is accurate up to the statistical precision $\epsilon_{\mbox{\tiny{stat}}}$, which is the rate expected for strongly convex functions. Furthermore, our techniques are more generally applicable than the methods proposed by previous authors and are not restricted to least-squares or even convex loss functions.

Problem Formulation

In this section, we develop some general theory for regularized $M$-estimators. We begin by establishing our notation and basic assumptions, before turning to the class of nonconvex regularizers and nonconvex loss functions to be covered in this paper.

To this end, we consider a regularized $M$-estimator of the form

2 Nonconvex Regularizers

On the nonnegative real line, the function $\rho_{\lambda}$ is nondecreasing.

For $t>0$, the function $t\mapsto\frac{\rho_{\lambda}(t)}{t}$ is nonincreasing in $t$.

The function $\rho_{\lambda}$ is differentiable for all $t\neq 0$ and subdifferentiable at $t=0$, with $\lim_{t\rightarrow 0^{+}}\rho_{\lambda}^{\prime}(t)=\lambda L$.

There exists $\mu>0$ such that $\rho_{\lambda,\mu}(t):=\rho_{\lambda}(t)+\frac{\mu}{2}t^{2}$ is convex.

SCAD penalty: This penalty, due to Fan and Li 2001, takes the form

where $a>2$ is a fixed parameter. As verified in Lemma 6 of Appendix A.2, the SCAD penalty satisfies the conditions of Assumption 1 with $L=1$ and $\mu=\frac{1}{a-1}$.

MCP regularizer: This penalty, due to Zhang 2010, takes the form

where $b>0$ is a fixed parameter. As verified in Lemma 7 in Appendix A.2, the MCP regularizer satisfies the conditions of Assumption 1 with $L=1$ and $\mu=\frac{1}{b}$.

3 Nonconvex Loss Functions and Restricted Strong Convexity

Throughout this paper, we require the loss function $\mathcal{L}_{n}$ to be differentiable, but we do not require it to be convex. Instead, we impose a weaker condition known as restricted strong convexity (RSC). Such conditions have been discussed in previous literature (Negahban et al. 2012; Agarwal et al. 2012), and involve a lower bound on the remainder in the first-order Taylor expansion of $\mathcal{L}_{n}$. In particular, our main statistical result is based on the following RSC condition:

where the $\alpha_{j}$’s are strictly positive constants and the $\tau_{j}$’s are nonnegative constants.

To understand this condition, note that if $\mathcal{L}_{n}$ were actually strongly convex, then both these RSC inequalities would hold with $\alpha_{1}=\alpha_{2}>0$ and $\tau_{1}=\tau_{2}=0$. However, in the high-dimensional setting ($p\gg n)$, the empirical loss $\mathcal{L}_{n}$ will not in general be strongly convex or even convex, but the RSC condition may still hold with strictly positive $(\alpha_{j},\tau_{j})$. In fact, if $\mathcal{L}_{n}$ is convex (but not strongly convex), the left-hand expression in (4b) is always nonnegative, so (4a) and (4b) hold trivially for $\frac{\|\Delta\|_{1}}{\|\Delta\|_{2}}\geq\sqrt{\frac{\alpha_{1}n}{\tau_{1}\log p}}$ and $\frac{\|\Delta\|_{1}}{\|\Delta\|_{2}}\geq\frac{\alpha_{2}}{\tau_{2}}\sqrt{\frac{n}{\log p}}$, respectively. Hence, the RSC inequalities only enforce a type of strong convexity condition over a cone of the form $\left\{\frac{\|\Delta\|_{1}}{\|\Delta\|_{2}}\leq c\sqrt{\frac{n}{\log p}}\right\}$.

It is important to note that the class of functions satisfying RSC conditions of this type is much larger than the class of convex functions; for instance, our own past work (Loh and Wainwright 2012) exhibits a large family of nonconvex quadratic functions that satisfy the condition (see Section 3.2 below for further discussion). Furthermore, note that we have stated two separate RSC inequalities (4b) for different ranges of $\|\Delta\|_{2}$, unlike in past work (Negahban et al. 2012; Agarwal et al. 2012; Loh and Wainwright 2012). As illustrated in the corollaries of Sections 3.3 and 3.4 below, an equality of the first type (4a) will only hold locally over $\Delta$ when we have more complicated types of loss functions that are only quadratic around a neighborhood of the origin. As proved in Appendix B.1, however, (4b) is implied by (4a) in cases when $\mathcal{L}_{n}$ is convex, which sustains our theoretical conclusions even under the weaker RSC conditions (4b). Further note that by the inequality

Finally, we clarify that whereas Negahban et al. 2012 define an RSC condition with respect to a fixed subset $S\subseteq\{1,\dots,p\}$, we follow the setup of Agarwal et al. 2012 and Loh and Wainwright 2012 and essentially require an RSC condition of the type defined in Negahban et al. 2012 to hold uniformly over all subsets $S$ of size $k$. Although the results on statistical consistency may be established under the weaker RSC assumption with $S:=\operatorname{supp}(\beta^{*})$, a uniform RSC condition is preferred because the true support set is not known a priori. The uniform RSC condition may be shown to hold w.h.p. in the sub-Gaussian settings we consider here (cf. Sections 3.2—3.4 below); in fact, the proofs contained in Negahban et al. 2012 establish a uniform RSC condition, as well.

Statistical Guarantees and Consequences

When $\widetilde{\beta}$ lies in the interior of the constraint set, this condition reduces to the usual zero-subgradient condition:

Such vectors $\widetilde{\beta}$ satisfying the condition (5) are also known as stationary points (Bertsekas 1999); note that the set of stationary points also includes interior local maxima. Hence, although some of the discussion below is stated in terms of “local minima,” the results hold for interior local maxima, as well.

Suppose the regularizer $\rho_{\lambda}$ satisfies Assumption 1, the empirical loss $\mathcal{L}_{n}$ satisfies the RSC conditions (4b) with $\frac{3}{4}\mu<\alpha_{1}$, and $\beta^{*}$ is feasible for the objective. Consider any choice of $\lambda$ such that

and suppose $n\geq\frac{16R^{2}\max(\tau_{1}^{2},\tau_{2}^{2})}{\alpha_{2}^{2}}\log p$. Then any vector $\widetilde{\beta}$ satisfying the first-order necessary conditions (5) satisfies the error bounds

Our next theorem provides a bound on a measure of the prediction error, as defined by the quantity

When the empirical loss $\mathcal{L}_{n}$ is a convex function, this measure is always nonnegative, and in various special cases, it has a form that is readily interpretable. For instance, in the case of the least-squares objective function $\mathcal{L}_{n}(\beta)=\frac{1}{2n}\|y-X\beta\|_{2}^{2}$, we have

corresponding to the usual measure of (fixed design) prediction error for a linear regression problem (cf. Corollary 1 below). More generally, when the loss function is the negative log likelihood for a generalized linear model with cumulant function $\psi$, the error measure (8) is equivalent to the symmetrized Bregman divergence defined by $\psi$. (See Section 3.3 for further details.)

Under the same conditions as Theorem 1, the error measure (8) is bounded as

This result shows that the prediction error (8) behaves similarly to the squared Euclidean norm between $\widetilde{\beta}$ and $\beta^{*}$.

We return to the proofs of Theorems 1 and 2 in Section 3.5. First, we develop various consequences of these theorems for various nonconvex loss functions and regularizers of interest. The main technical challenge is to establish that the RSC conditions (4b) hold with high probability for appropriate choices of positive constants $\{(\alpha_{j},\tau_{j})\}_{j=1}^{2}$.

2 Corrected Linear Regression

We begin by considering the case of high-dimensional linear regression with systematically corrupted observations. Recall that in the framework of ordinary linear regression, we have the linear model

We use the population and empirical loss functions

where $(\widehat{\Gamma},\widehat{\gamma})$ are estimators for $(\Sigma_{x},\Sigma_{x}\beta^{*})$ that depend only on $\{(z_{i},y_{i})\}_{i=1}^{n}$. It is easy to see that $\beta^{*}=\arg\min_{\beta}\mathcal{L}(\beta)$. From the formulation (1), the corrected linear regression estimator is given by

We now state a concrete corollary in the case of additive noise (model (a) above). In this case, as discussed in Loh and Wainwright 2012, an appropriate choice of the pair $(\widehat{\Gamma},\widehat{\gamma})$ is given by

Here, we assume the noise covariance $\Sigma_{w}$ is known or may be estimated from replicates of the data. Such an assumption also appears in canonical errors-in-variables literature (Carroll et al. 1995), but it is an open question how to devise a corrected estimator when an estimate of $\Sigma_{w}$ is not readily available. If we assume a sub-Gaussian model on the covariates and errors (i.e., $x_{i}$, $w_{i}$, and $\epsilon_{i}$ are sub-Gaussian with parameters $\sigma_{x}^{2}$, $\sigma_{w}^{2}$, and $\sigma_{\epsilon}^{2}$, respectively), the contribution of the error covariances may be summarized in the error term

which appears as a prefactor in the deviation bounds and estimation/prediction error bounds for the subsequent estimators (cf. Lemma 2 in Loh and Wainwright 2012). We make this dependence explicit in the statement of the corollary for high-dimensional errors-in-variables regression below. Note in particular that $\varphi$ scales up with both $\sigma_{\epsilon}$ and $\sigma_{w}$. Hence, even when $\sigma_{\epsilon}=0$, corresponding to no additive error, we will have $\varphi\neq 0$ due to errors in the covariates; whereas when $\sigma_{w}=0$, corresponding to cleanly observed covariates, we will still have $\varphi\neq 0$ due to the additional additive error introduced by the $\epsilon_{i}$’s, agreeing with canonical results for the Lasso (Bickel et al. 2009).

In the high-dimensional setting ($p\gg n$), the matrix $\widehat{\Gamma}$ in (13) is always negative definite: the matrix $\frac{Z^{T}Z}{n}$ has rank at most $n$, and the positive definite matrix $\Sigma_{w}$ is then subtracted to obtain $\widehat{\Gamma}$. Consequently, the empirical loss function $\mathcal{L}_{n}$ previously defined (11) is nonconvex. Other choices of $\widehat{\Gamma}$ are applicable to missing data (model (b)), and also lead to nonconvex programs (see Loh and Wainwright 2012 for further details).

Suppose we have i.i.d. observations $\{(z_{i},y_{i})\}_{i=1}^{n}$ from a corrupted linear model with additive noise, where the covariates and error terms are sub-Gaussian. Let $\varphi$ be defined as in (14) with respect to the sub-Gaussian parameters. Suppose $(\lambda,R)$ are chosen such that $\beta^{*}$ is feasible and

Also suppose $\frac{3}{4}\mu<\frac{1}{2}\lambda_{\min}(\Sigma_{x})$. Then given a sample size $n\geq C\,\max\{R^{2},k\}\log p$, any stationary point $\widetilde{\beta}$ of the nonconvex program (12) satisfies the estimation error bounds

with probability at least $1-c_{1}\exp(-c_{2}\log p)$, where $\|\beta^{*}\|_{0}=k$.

Furthermore, our theory provides a theoretical motivation for why the usual choice of $a=3.7$ for linear regression with the SCAD penalty (Fan and Li 2001) is reasonable. Indeed, as discussed in Section 2.2, we have

in that case. Since $x_{i}\sim N(0,I)$ in the SCAD simulations, we have $\frac{3}{4}\mu<\frac{1}{2}\lambda_{\min}(\Sigma_{x})$ for the choice $a=3.7$. For further comments regarding the parameter $a$ in the SCAD penalty, see the discussion concerning Figure 3 in Section 5.

3 Generalized Linear Models

Moving beyond linear regression, we now consider the case where observations are drawn from a generalized linear model (GLM). Recall that a GLM is characterized by the conditional distribution

where $\sigma>0$ is a scale parameter and $\psi$ is the cumulant function, By standard properties of exponential families (McCullagh and Nelder 1989; Lehmann and Casella 1998), we have

The population loss corresponding to the negative log likelihood is then given by

giving rise to the population-level and empirical gradients

Since we are optimizing over $\beta$, we will rescale the loss functions and assume $c(\sigma)=1$. We may check that if $\beta^{*}$ is the true parameter of the GLM, then $\nabla\mathcal{L}(\beta^{*})=0$; furthermore,

We will assume that $\beta^{*}$ is sparse and optimize the penalized maximum likelihood program

We then have the following corollary, proved in Appendix B.3:

Suppose we have i.i.d. observations $\{(x_{i},y_{i})\}_{i=1}^{n}$ from a GLM, where the $x_{i}$’s are sub-Gaussian. Suppose $(\lambda,R)$ are chosen such that $\beta^{*}$ is feasible and

Then given a sample size $n\geq CR^{2}\log p$, any stationary point $\widetilde{\beta}$ of the nonconvex program (15) satisfies

with probability at least $1-c_{1}\exp(-c_{2}\log p)$, where $\|\beta^{*}\|_{0}=k$. Here, $\alpha_{1}$ is a constant depending on $\|\beta^{*}\|_{2}$, $\psi$, $\lambda_{\min}(\Sigma_{x})$, and the sub-Gaussian parameter of the $x_{i}$’s, and we assume $\mu<2\alpha_{1}$.

Although $\mathcal{L}_{n}$ is convex in this case, the overall program may not be convex if the regularizer $\rho_{\lambda}$ is nonconvex, giving rise to multiple local optima. For instance, see the simulations of Figure 4 in Section 5 for a demonstration of such local optima. In past work, Breheny and Huang 2011 studied logistic regression with SCAD and MCP regularizers, but did not provide any theoretical results on the quality of the local optima. In this context, Corollary 2 shows that their coordinate descent algorithms are guaranteed to converge to a stationary point $\widetilde{\beta}$ within close proximity of the true parameter $\beta^{*}$.

In the statement of Corollary 2, we choose not to write out the form of $\alpha_{1}$ explicitly as in Corollary 1, since it is rather complicated. As explained in the proof of Corollary 2 in Appendix B.3, the precise form of $\alpha_{1}$ may be traced back to Proposition 2 of Negahban et al. 2012.

4 Graphical Lasso

Finally, we specialize our results to the case of the graphical Lasso. Given $p$-dimensional observations $\{x_{i}\}_{i=1}^{n}$, the goal is to estimate the structure of the underlying (sparse) graphical model. Recall that the population and empirical losses for the graphical Lasso are given by

where $\widehat{\Sigma}$ is an empirical estimate for the covariance matrix $\Sigma=\operatorname{Cov}(x_{i})$. The objective function for the graphical Lasso is then given by

As suggested by Loh and Wainwright 2013a, the graphical Lasso easily accommodates systematically corrupted observations, with the only modification being the form of the sample covariance matrix $\widehat{\Sigma}$. Just as in Corollary 1, the magnitude and form of corruption would occur as a prefactor in the deviation condition captured in (17) below; for instance, in the case of $\widehat{\Sigma}=\frac{Z^{T}Z}{n}-\Sigma_{w}$, corresponding to additive noise in the $x_{i}$’s, the bound (17) would involve a prefactor of $\sigma_{z}^{2}$ rather than $\sigma_{x}^{2}$, where $\sigma_{z}^{2}$ and $\sigma_{x}^{2}$ are the sub-Gaussian parameters of $z_{i}$ and $x_{i}$, respectively.

Further note that the program (16) is always useful for obtaining a consistent estimate of a sparse inverse covariance matrix, regardless of whether the $x_{i}$’s are drawn from a distribution for which $\Theta^{*}$ is relevant in estimating the edges of the underlying graph. Note that other variants of the graphical Lasso exist in which only off-diagonal entries of $\Theta$ are penalized, and similar results for statistical consistency hold in that case. Here, we assume that all entries are penalized equally in order to simplify our arguments. The same framework is considered by Fan et al. 2009.

We have the following result, proved in Appendix B.4. The statement of the corollary is purely deterministic, but in cases of interest (say, sub-Gaussian observations), the deviation condition (17) holds with probability at least $1-c_{1}\exp(-c_{2}\log p)$, translating into the Frobenius norm bound (18) holding with the same probability.

Suppose we have an estimate $\widehat{\Sigma}$ of the covariance matrix $\Sigma$ based on (possibly corrupted) observations $\{x_{i}\}_{i=1}^{n}$, such that

Also suppose $\Theta^{*}$ has at most $s$ nonzero entries. Suppose $(\lambda,R)$ are chosen such that $\Theta^{*}$ is feasible and

Suppose $\frac{3}{4}\mu<\left(\left|\!\left|\!\left|{\Theta^{*}}\right|\!\right|\!\right|_{2}+1\right)^{-2}$. Then with a sample size $n>Cs\log p$, for a sufficiently large constant $C>0$, any stationary point $\widetilde{\Theta}$ of the nonconvex program (16) satisfies

5 Proof of Theorems 1 and 2

We now turn to the proofs of our two main theorems.

Proof of Theorem 1: Introducing the shorthand $\widetilde{\nu}:=\widetilde{\beta}-\beta^{*}$, we begin by proving that $\|\widetilde{\nu}\|_{2}\leq 1$. If not, then (4b) gives the lower bound

Since $\beta^{*}$ is feasible, we may take $\beta=\beta^{*}$ in (5), and combining with (19) yields

By Hölder’s inequality, followed by the triangle inequality, we also have

where inequality (i) follows since $\|\nabla\mathcal{L}_{n}(\beta^{*})\|_{\infty}\leq\frac{\lambda L}{2}$ by the bound (6), and $\|\nabla\rho_{\lambda}(\widetilde{\beta})\|_{\infty}\leq\lambda L$ by Lemma 4 in Appendix A.1. Combining this upper bound with (20) and rearranging then yields

By our choice of $\lambda$ from (6) and the assumed lower bound on the sample size $n$, the right hand side is at most $1$, so $\|\widetilde{\nu}\|_{2}\leq 1$, as claimed.

Consequently, we may apply (4a), yielding the lower bound

Since the function $\rho_{\lambda,\mu}(\beta):=\rho_{\lambda}(\beta)+\frac{\mu}{2}\|\beta\|_{2}^{2}$ is convex by assumption, we have

Combining (21) with (5) and (22), we obtain

Rearranging and using Hölder’s inequality, we then have

Combining this with (3.5) and (52) in Lemma 4 in Appendix A.1, as well as the subadditivity of $\rho_{\lambda}$, we then have

In particular, we have $3\rho_{\lambda}(\beta^{*})-\rho_{\lambda}(\widetilde{\beta})\geq 0$, so we may apply Lemma 5 in Appendix A.1 to conclude that

where $A$ denotes the index set of the $k$ largest elements of $\widetilde{\beta}-\beta^{*}$ in magnitude. In particular, we have the cone condition

Substituting (25) into (24), we then have

Proof of Theorem 2: In order to establish (9), note that combining the first-order condition (5) with the upper bound (22), we have

Furthermore, as noted earlier, Lemma 4 in Appendix A.1 implies that

Optimization Algorithms

We now describe how a version of composite gradient descent (Nesterov 2007) may be applied to efficiently optimize the nonconvex program (1), and show that it enjoys a linear rate of convergence under suitable conditions. In this section, we focus exclusively on a version of the optimization problem with the side function

Note that this choice of $g_{\lambda,\mu}$ is convex by Assumption 1. We may then write the program (1) as

In this way, the objective function decomposes nicely into a sum of a differentiable but nonconvex function and a possibly nonsmooth but convex penalty. Applied to the representation (29) of the objective function, the composite gradient descent procedure of Nesterov 2007 produces a sequence of iterates $\{\beta^{t}\}_{t=0}^{\infty}$ via the updates

where $\frac{1}{\eta}$ is the stepsize. As discussed in Section 4.2, these updates may be computed in a relatively straightforward manner.

The main result of this section is to establish that the algorithm defined by the iterates (30) converges very quickly to a $\delta$-neighborhood of any global optimum, for all tolerances $\delta$ that are of the same order (or larger) than the statistical error.

We begin by setting up the notation and assumptions underlying our result. The Taylor error around the vector $\beta_{2}$ in the direction $\beta_{1}-\beta_{2}$ is given by

We analogously define the Taylor error $\makebox[0.0pt][l]{\hskip 2.05pt\hskip 0.0pt\rule[8.12498pt]{4.91359pt}{0.43057pt}}{{\mathcal{T}}}$ for the modified loss function $\makebox[0.0pt][l]{\hskip 2.05pt\hskip 0.0pt\rule[8.12498pt]{3.98999pt}{0.43057pt}}{\mathcal{L}}_{n}$, and note that

a condition referred to as restricted smoothness in past work (Agarwal et al. 2012). Throughout this section, we assume $2\alpha_{i}>\mu$ for all $i$, where $\mu$ is the coefficient ensuring the convexity of the function $g_{\lambda,\mu}$ from (28). Furthermore, we define $\alpha=\min\{\alpha_{1},\alpha_{2}\}$ and $\tau=\max\{\tau_{1},\tau_{2},\tau_{3}\}$.

Under the stated scaling on the sample size, we are guaranteed that $\kappa\in(0,1)$, so it is a contraction factor. Roughly speaking, we show that the squared optimization error will fall below $\delta^{2}$ within $T\asymp\frac{\log(1/\delta^{2})}{\log(1/\kappa)}$ iterations. More precisely, our theorem guarantees $\delta$-accuracy for all iterations larger than

where $\phi(\beta):=\mathcal{L}_{n}(\beta)+\rho_{\lambda}(\beta)$ denotes the composite objective function. As clarified in the theorem statement, the squared tolerance $\delta^{2}$ is not allowed to be arbitrarily small, which would contradict the fact that the composite gradient method may converge to a stationary point. However, our theory allows $\delta^{2}$ to be of the same order as the squared statistical error $\epsilon_{\mbox{\tiny{stat}}}^{2}=\|\widehat{\beta}-\beta^{*}\|_{2}^{2}$, the distance between a fixed global optimum and the target parameter $\beta^{*}$. From a statistical perspective, there is no point in optimizing beyond this tolerance.

With this setup, we now turn to a precise statement of our main optimization-theoretic result. As with Theorems 1 and 2, the statement of Theorem 3 is entirely deterministic.

Suppose the empirical loss $\mathcal{L}_{n}$ satisfies the RSC/RSM conditions (33b) and (34), and suppose the regularizer $\rho_{\lambda}$ satisfies Assumption 1. Suppose $\widehat{\beta}$ is any global minimum of the program (29), with regularization parameters chosen such that

Suppose $\mu<2\alpha$. Then for any stepsize parameter $\eta\geq\max\{2\alpha_{3}-\mu,\,\mu\}$ and tolerance $\delta^{2}\geq\frac{c\epsilon_{\mbox{\tiny{stat}}}^{2}}{1-\kappa}\cdot\frac{k\log p}{n}$, we have

Remark: Note that for the optimal choice of tolerance parameter $\delta\asymp\frac{k\log p}{n}\epsilon_{\mbox{\tiny{stat}}}$, the error bound appearing in (37) takes the form $\frac{c\epsilon_{\mbox{\tiny{stat}}}^{2}}{2\alpha-\mu}\cdot\frac{k\log p}{n}$, meaning that successive iterates of the composite gradient descent algorithm are guaranteed to converge to a region within statistical accuracy of the true global optimum $\widehat{\beta}$. Concretely, if the sample size satisfies $n\succsim Ck\log p$ and the regularization parameters are chosen appropriately, Theorem 1 guarantees that $\epsilon_{\mbox{\tiny{stat}}}={\mathcal{O}}\left(\sqrt{\frac{k\log p}{n}}\right)$ with high probability. Combined with Theorem 3, we then conclude that

for all iterations $t\geq T(\epsilon_{\mbox{\tiny{stat}}})$.

As would be expected, the (restricted) curvature $\alpha$ of the loss function and nonconvexity parameter $\mu$ of the penalty function enter into the bound via the denominator $2\alpha-\mu$. Indeed, the bound is tighter when the loss function possesses more curvature or the penalty function is closer to being convex, agreeing with intuition. Similar to our discussion in the remark following Theorem 2, the requirement $\mu<2\alpha$ is certainly necessary for our proof technique, but it is possible that composite gradient descent still produces good results when this condition is violated. See Section 5 for simulations in scenarios involving mild and severe violations of this condition.

Finally, note that the parameter $\eta$ must be sufficiently large (or equivalently, the stepsize must be sufficiently small) in order for the composite gradient descent algorithm to be well-behaved. See Nesterov 2007 for a discussion of how the stepsize may be chosen via an iterative search when the problem parameters are unknown.

In the case of corrected linear regression (Corollary 1), Lemma 13 of Loh and Wainwright 2012 establishes the RSC/RSM conditions for various statistical models. The following proposition shows that the conditions (33b) and (34) hold in GLMs when the $x_{i}$’s are drawn i.i.d. from a zero-mean sub-Gaussian distribution with parameter $\sigma_{x}^{2}$ and covariance matrix $\Sigma=\cov(x_{i})$. As usual, we assume a sample size $n\geq c\,k\log p$, for a sufficiently large constant $c>0$. Recall the definition of the Taylor error ${\mathcal{T}}(\beta_{1},\beta_{2})$ from (31).

with probability at least $1-c_{1}\exp(-c_{2}n)$. With the bound $\|\psi^{\prime\prime}\|_{\infty}\leq\alpha_{u}$, we also have

with probability at least $1-c_{1}\exp(-c_{2}n)$.

For the proof of Proposition 1, see Appendix D.

2 Form of Updates

If $g_{\lambda,\mu}(\widehat{\beta})\leq R$, define $\beta^{t+1}=\widehat{\beta}$.

Otherwise, if $g_{\lambda,\mu}(\widehat{\beta})>R$, optimize the constrained program

We derive the correctness of this procedure in Appendix C.1. For many nonconvex regularizers $\rho_{\lambda}$ of interest, the unconstrained program (40) has a convenient closed-form solution: For the SCAD penalty (2), the program (40) has simple closed-form solution given by

For the MCP (3), the optimum of the program (40) takes the form

and the operations are taken componentwise. See Appendix C.2 for the derivation of these closed-form updates.

3 Proof of Theorem 3

We provide the outline of the proof here, with more technical results deferred to Appendix C. In broad terms, our proof is inspired by a result of Agarwal et al. 2012, but requires various modifications in order to be applied to the much larger family of nonconvex regularizers considered here.

Our first lemma shows that the optimization error $\beta^{t}-\widehat{\beta}$ lies in an approximate cone set:

Under the conditions of Theorem 3, suppose there exists a pair $(\bar{\eta},T)$ such that

Then for any iteration $t\geq T$, we have

Our second lemma shows that as long as the composite gradient descent algorithm is initialized with a solution $\beta^{0}$ within a constant radius of a global optimum $\widehat{\beta}$, all successive iterates also lie within the same ball:

Under the conditions of Theorem 3, and with an initial vector $\beta^{0}$ such that $\|\beta^{0}-\widehat{\beta}\|_{2}\leq 3$, we have

In particular, suppose we initialize the composite gradient procedure with a vector $\beta^{0}$ such that $\|\beta^{0}\|_{2}\leq\frac{3}{2}$. Then by the triangle inequality,

where we have assumed our scaling of $n$ guarantees $\|\widehat{\beta}-\beta^{*}\|_{2}\leq 1/2$.

Finally, recalling our earlier definition (35) of $\kappa$, the third lemma combines the results of Lemmas 1 and 2 to establish a bound on the value of the objective function that decays exponentially with $t$:

Under the same conditions of Lemma 2, suppose in addition that (44) holds and $\frac{32k\tau\log p}{n}\leq\frac{2\alpha-\mu}{4}$. Then for any $t\geq T$, we have

where $\makebox[0.0pt][l]{\hskip 1.29167pt\hskip 0.0pt\rule[5.59721pt]{2.62898pt}{0.43057pt}}{\epsilon}:=8\sqrt{k}\epsilon_{\mbox{\tiny{stat}}}$, $\epsilon:=2\cdot\min\left(\frac{2\bar{\eta}}{\lambda L},R\right)$, the quantities $\kappa$ and $\varphi$ are defined according to (35), and

The remainder of the proof follows an argument used in Agarwal et al. 2012, so we only provide a high-level sketch. We first prove the following inequality:

In the first iteration, we apply Lemma 3 with $\bar{\eta}_{0}=\phi(\beta^{0})-\phi(\widehat{\beta})$ to obtain

Let $\bar{\eta}_{1}:=\frac{2\xi}{1-\kappa}(4R^{2}+\makebox[0.0pt][l]{\hskip 1.29167pt\hskip 0.0pt\rule[5.59721pt]{2.62898pt}{0.43057pt}}{\epsilon}^{2})$, and note that for $T_{1}:=\Bigg\lceil\frac{\log(2\bar{\eta}_{0}/\bar{\eta}_{1})}{\log(1/\kappa)}\Bigg\rceil$, we have

Finally, by (84b) in the proof of Lemma 3 in Appendix C.5 and the relative scaling of $(n,p,k)$, we have

Simulations

Linear regression: In the case of linear regression, we simulated covariates corrupted by additive noise according to the mechanism described in Section 3.2, giving the estimator

We generated i.i.d. samples $x_{i}\sim N(0,I)$ and set $\Sigma_{w}=(0.2)^{2}I$, and generated additive noise $\epsilon_{i}\sim N(0,(0.1)^{2})$.

Logistic regression: In the case of logistic regression, we also generated i.i.d. samples $x_{i}\sim N(0,I)$. Since $\psi(t)=\log(1+\exp(t))$, the program (15) becomes

We optimized the programs (50) and (51) using the composite gradient updates (30). In order to compute the updates, we used the three-step procedure described in Section 4.2, together with the updates for SCAD and MCP given by (42) and (43). Note that the updates for the Lasso penalty may be generated more simply and efficiently as discussed in Agarwal et al. 2012.

Figure 4 provides analogous results to Figure 3 in the case of logistic regression, using $p=64,k=\lfloor\sqrt{p}\rfloor$, and $n=\lfloor 20k\log p\rfloor$. The plot shows solution trajectories for 20 different initializations of composite gradient descent. Again, we see that the log optimization error decreases at a linear rate up to the level of statistical error, as predicted by Theorem 3. Furthermore, the Lasso penalty yields a unique global optimum $\widehat{\beta}$, since the program (51) is convex, as we observe in panel (a). In contrast, the nonconvex program based on the SCAD penalty produces multiple local optima, whereas the MCP yields a relatively large number of local optima. Note that empirically, all local optima appear to lie within the small ball around $\beta^{*}$ defined in Theorem 1. However, if we use $\lambda_{\min}(\nabla^{2}\mathcal{L}_{n}(\beta^{*}))$ as a surrogate for $\alpha_{1}$, we see that $2\alpha_{1}<\mu$ in the case of the SCAD or MCP regularizers, which is not covered by our theory.

Discussion

We have analyzed theoretical properties of local optima of regularized $M$-estimators, where both the loss and penalty function are allowed to be nonconvex. Our results are the first to establish that all stationary points of such nonconvex problems are close to the truth, implying that any optimization method guaranteed to converge to a stationary point will provide statistically consistent solutions. We show concretely that a variant of composite gradient descent may be used to obtain near-global optima in linear time, and verify our theoretical results with simulations.

Appendix A Properties of Regularizers

In this section, we establish properties of some nonconvex regularizers covered by our theory (Appendix A.1) and verify that specific regularizers satisfy Assumption 1 (Appendix A.2). The properties given in Appendix A.1 are used in the proof of Theorem 1.

We begin with some general properties of regularizers that satisfy Assumption 1.

Under conditions (i)–(ii) of Assumption 1, conditions (iii) and (iv) together imply that $\rho_{\lambda}$ is $\lambda L$-Lipschitz as a function of $t$. In particular, all subgradients and derivatives of $\rho_{\lambda}$ are bounded in magnitude by $\lambda L$.

Under the conditions of Assumption 1, we have

(a): Suppose $0\leq t_{1}\leq t_{2}$. Then

by condition (iii). Applying (iii) once more, we have

where the last equality comes from condition (iv). Hence,

A similar argument applies to the cases when one (or both) of $t_{1}$ and $t_{2}$ are negative.

(b): Clearly, it suffices to verify the inequality for the scalar case:

The inequality is trivial for $t=0$. For $t>0$, the convexity of the right-hand expression implies that for any $s\in(0,t)$, we have

Taking a limit as $s\rightarrow 0^{+}$ then yields the desired inequality. The case $t<0$ follows by symmetry. ∎

where $\nu:=\beta-\beta^{*}$ and $A$ is the index set of the $k$ largest elements of $\nu$ in magnitude.

We first establish (53). Define $f(t):=\frac{t}{\rho_{\lambda}(t)}$ for $t>0$. By our assumptions on $\rho_{\lambda}$, the function $f$ is nondecreasing in $|t|$, so

Again using the nondecreasing property of $f$, we have

where the last equality follows from condition (iv) of Assumption 1. Combining this result with (55) and (56) yields

We now turn to the proof of the bound (54). Letting $S:=\operatorname{supp}(\beta^{*})$ denote the support of $\beta^{*}$, the triangle inequality and subadditivity of $\rho$ (see the remark following Assumption 1; cf. Lemma 1 of Chen and Gu 2014) imply that

A.2 Verification for Specific Regularizers

We now verify that Assumption 1 is satisfied by the SCAD and MCP regularizers. (The properties are trivial to verify for the Lasso penalty.)

The SCAD regularizer (2) with parameter $a$ satisfies the conditions of Assumption 1 with $L=1$ and $\mu=\frac{1}{a-1}$.

Conditions (i)–(iii) were already verified in Zhang and Zhang 2012. Furthermore, we may easily compute the derivative of the SCAD regularizer to be

and any point in the interval $[-\lambda,\lambda]$ is a valid subgradient at $t=0$, so condition (iv) is satisfied for any $L\geq 1$. Furthermore, we have $\frac{\partial^{2}}{\partial t^{2}}\rho_{\lambda}(t)\geq\frac{-1}{a-1}$, so $\rho_{\lambda,\mu}$ is convex whenever $\mu\geq\frac{1}{a-1}$, giving condition (v). ∎

The MCP regularizer (3) with parameter $b$ satisfies the conditions of Assumption 1 with1 $L=1$ and $\mu=\frac{1}{b}$.

Again, the conditions (i)–(iii) are already verified in Zhang and Zhang 2012. We may compute the derivative of the MCP regularizer to be

with subgradient $\lambda[-1,+1]$ at $t=0$, so condition (iv) is again satisfied for any $L\geq 1$. Taking another derivative, we have $\frac{\partial^{2}}{\partial t^{2}}\rho_{\lambda}(t)\geq\frac{-1}{b}$, so condition (v) of Assumption 1 holds with $\mu=\frac{1}{b}$. ∎

Appendix B Proofs of Corollaries in Section 3

In this section, we provide proofs of the corollaries to Theorem 1 stated in Section 3. Throughout this section, we use the convenient shorthand notation

We begin with two lemmas that will be useful for establishing the RSC conditions (4b) in the special case where $\mathcal{L}_{n}$ is convex. We assume throughout that $\|\Delta\|_{1}\leq 2R$, since $\beta^{*}$ and $\beta^{*}+\Delta$ lie in the feasible set.

Suppose $\mathcal{L}_{n}$ is convex. If condition (4a) holds and $n\geq 4R^{2}\tau_{1}^{2}\log p$, then

Taking $t=\frac{1}{\|\Delta\|_{2}}\in(0,1]$ and applying condition (4a) to the rescaled vector $\frac{\Delta}{\|\Delta\|_{2}}$ then yields

where the third inequality uses the assumption on the relative scaling of $(n,p)$ and the fact that $\|\Delta\|_{2}\geq 1$. ∎

again using the assumption on the scaling of $(n,p)$. ∎

B.2 Proof of Corollary 1

Note that $\mathcal{E}_{n}(\Delta)=\Delta^{T}\widehat{\Gamma}\Delta$, so in particular,

Applying Lemma 12 in Loh and Wainwright 2012 with $s=\frac{n}{\log p}$ to bound the second term, we have

As shown in previous work (Loh and Wainwright 2012), both of these terms are upper-bounded by $c^{\prime}\,\varphi\sqrt{\frac{\log p}{n}}$ with high probability. Consequently, the claim in the corollary follows by applying Theorem 1.

B.3 Proof of Corollary 2

Applying the mean value theorem, we find that

where $t_{i}\in$. From (the proof of) Proposition 2 in Negahban et al. 2012, we then have

with probability at least $1-c_{1}\exp(-c_{2}n)$, for an appropriate choice of $\alpha_{1}$. Note that by the arithmetic mean-geometric mean inequality,

which establishes (4a). Inequality (4b) then follows via Lemma 8 in Appendix B.1.

It remains to show that there are universal constants $(c,c_{1},c_{2})$ such that

For each $1\leq i\leq n$ and $1\leq j\leq p$, define the random variable $V_{ij}:=(\psi^{\prime}(x_{i}^{T}\beta^{*})-y_{i})x_{ij}$. Our goal is to bound $\max_{j=1,\ldots,p}|\frac{1}{n}\sum_{i=1}^{n}V_{ij}|$. Note that

using the fact that $\psi$ is the cumulant generating function for the underlying exponential family. Thus, by a Taylor series expansion, there is some $v_{i}\in$ such that

B.4 Proof of Corollary 3

We first verify condition (4a) in the case where $|\!|\!|\Delta|\!|\!|_{{F}}\leq 1$. A straightforward calculation yields

for some $t\in$. By standard properties of the Kronecker product (Horn and Johnson 1990), we have

using the fact that $\left|\!\left|\!\left|{\Delta}\right|\!\right|\!\right|_{2}\leq\left|\!\left|\!\left|{\Delta}\right|\!\right|\!\right|_{F}\leq 1$. Plugging back into (65) yields

so (4a) holds with $\alpha_{1}=\left(\left|\!\left|\!\left|{\Theta^{*}}\right|\!\right|\!\right|_{2}+1\right)^{-2}$ and $\tau_{1}=0$. Lemma 9 then implies (4b) with $\alpha_{2}=\left(\left|\!\left|\!\left|{\Theta^{*}}\right|\!\right|\!\right|_{2}+1\right)^{-2}$. Finally, we need to establish that the given choice of $\lambda$ satisfies the requirement (6) of Theorem 1. By the assumed deviation condition (17), we have

Applying Theorem 1 then implies the desired result.

Appendix C Auxiliary Optimization-Theoretic Results

In this section, we provide proofs of the supporting lemmas used in Section 4.

We begin by deriving the correctness of the three-step procedure given in Section 4.2. Let $\widehat{\beta}$ be the unconstrained optimum of the program (40). If $g_{\lambda,\mu}(\widehat{\beta})\leq R$, we clearly have the update given in step (2). Suppose instead that $g_{\lambda,\mu}(\widehat{\beta})>R$. Then since the program (30) is convex, the iterate $\beta^{t+1}$ must lie on the boundary of the feasible set; i.e.,

By Lagrangian duality, the program (30) is also equivalent to

In fact, since the projection will shrink the vector to the boundary of the constraint set, (66) forces $R^{\prime}=R$. This yields the update (41) appearing in step (3).

C.2 Derivation of Updates for SCAD and MCP

We now derive the explicit form of the updates (42) and (43) for the SCAD and MCP regularizers, respectively. We may rewrite the unconstrained program (40) as

Since the program in the last line of equation (C.2) decomposes by coordinate, it suffices to solve the scalar optimization problem

We first consider the case when $\rho$ is the SCAD penalty. The solution $\widehat{x}$ of the program (68) in the case when $\nu=1$ is given in Fan and Li 2001; the expression (42) for the more general case comes from writing out the subgradient of the objective as

using the equation for the SCAD derivative (57), and setting the subgradient equal to zero.

Similarly, when $\rho$ is the MCP parametrized by $(b,\lambda)$, the subgradient of the objective takes the form

using the expression for the MCP derivative (58), leading to the closed-form solution given in (43). This agrees with the expression provided in Breheny and Huang 2011 for the special case when $\nu=1$.

C.3 Proof of Lemma 1

We first show that if $\lambda\geq\frac{8}{L}\cdot\|\nabla\mathcal{L}_{n}(\beta^{*})\|_{\infty}$, then for any feasible $\beta$ such that

Defining the error vector $\Delta:=\beta-\beta^{*}$, (69) implies

so subtracting $\langle\nabla\mathcal{L}_{n}(\beta^{*}),\,\Delta\rangle$ from both sides gives

We divide the argument into two cases. First suppose $\|\Delta\|_{2}\leq 3$. Note that if $\bar{\eta}\geq\frac{\lambda L}{4}\|\Delta\|_{1}$, the claim (70) is trivially true; so assume $\bar{\eta}\leq\frac{\lambda L}{4}\|\Delta\|_{1}$. Then the RSC condition (33a), together with (71), implies that

Rearranging and using the assumption $\lambda L\geq 16R\tau_{1}\frac{\log p}{n}$, along with Lemma 4 in Appendix A.1, we then have

by Lemma 5 in Appendix A.1. Furthermore, note that the bound (C.3) also implies that

In the case when $\|\Delta\|_{2}\geq 3$, the RSC condition (33b) gives

In particular, if $\rho_{\lambda}(\beta^{*})-\rho_{\lambda}(\beta^{*}+\Delta)\leq 0$, we have

a contradiction. Hence, using Lemma 5 in Appendix A.1, we have

Note that under the scaling $\lambda L\geq 4\tau_{2}\sqrt{\frac{\log p}{n}}$, the bound (C.3) also implies (74). Combining (74) and (76), we then have

Using the trivial bound $\|\Delta\|_{1}\leq 2R$, we obtain the claim (70).

We now apply the implication (69) to the vectors $\widehat{\beta}$ and $\beta^{t}$. Note that by optimality of $\widehat{\beta}$, we have

C.4 Proof of Lemma 2

Our proof proceeds via induction on the iteration number $t$. Note that the base case $t=0$ holds by assumption. Hence, it remains to show that if $\|\beta^{t}-\widehat{\beta}\|_{2}\leq 3$ for some integer $t\geq 1$, then $\|\beta^{t+1}-\widehat{\beta}\|_{2}\leq 3$, as well.

We assume for the sake of a contradiction that $\|\beta^{t+1}-\widehat{\beta}\|_{2}>3$. By the RSC condition (33b) and the relation (32), we have

Furthermore, by convexity of $g:=g_{\lambda,\mu}$, we have

Multiplying by $\lambda$ and summing with (77) then yields

Together with the first-order optimality condition $\langle\nabla\phi(\widehat{\beta}),\,\beta^{t+1}-\widehat{\beta}\rangle\geq 0$, we then have

Since $\|\widehat{\beta}-\beta^{t}\|_{2}\leq 3$ by the induction hypothesis, applying the RSC condition (33a) to the pair $(\widehat{\beta},\beta^{t})$ also gives

Finally, the RSM condition (34) on the pair $(\beta^{t+1},\beta^{t})$ gives

since $\frac{\eta}{2}\geq\alpha_{3}-\frac{\mu}{2}$ by assumption, and $\|\beta^{t+1}-\beta^{t}\|_{1}\leq 2R$. It is easy to check that the update (30) may be written equivalently as

and the optimality of $\beta^{t+1}$ then yields

Summing up (C.4), (81), and (83), we then have

Combining this last inequality with (79), we have

since $\|\beta^{t}-\widehat{\beta}\|_{2}\leq 3$ by the induction hypothesis and $\|\beta^{t+1}-\widehat{\beta}\|_{2}>3$ by assumption, and using the fact that $\eta\geq\mu$. It follows that

where the final inequality holds whenever $2R\tau\sqrt{\frac{\log p}{n}}+\frac{8R^{2}\tau\log p}{n}\leq 3\left(\alpha-\frac{3\mu}{2}\right)$. Rearranging gives $\|\beta^{t+1}-\widehat{\beta}\|_{2}\leq 3$, providing the desired contradiction.

C.5 Proof of Lemma 3

We prove this result later, taking it as given for the moment.

the objective function minimized over the constraint set $\{g(\beta)\leq R\}$ at iteration $t$. For any $\gamma\in$, the vector $\beta_{\gamma}:=\gamma\widehat{\beta}+(1-\gamma)\beta^{t}$ belongs to the constraint set, as well. Consequently, by the optimality of $\beta^{t+1}$ and feasibility of $\beta_{\gamma}$, we have

where inequality (i) incorporates the fact that

since $\alpha_{3}-\mu\leq\frac{\eta}{2}$ by assumption, and adding $\lambda g(\beta^{t+1})$ to both sides gives

where we have defined $\Delta^{t}:=\beta^{t}-\widehat{\beta}$. Combined with (C.5), we therefore have

Now introduce the shorthand $\delta_{t}:=\phi(\beta^{t})-\phi(\widehat{\beta})$ and $\upsilon(k,p,n)=k\tau\frac{\log p}{n}$. By applying (84b) and subtracting $\phi(\widehat{\beta})$ from both sides of (87), we have

Choosing $\gamma=\frac{2\alpha-\mu}{4\eta}\in(0,1)$ yields

or $\delta_{t+1}\leq\kappa\delta_{t}+\xi(\epsilon+\makebox[0.0pt][l]{\hskip 1.29167pt\hskip 0.0pt\rule[5.59721pt]{2.62898pt}{0.43057pt}}{\epsilon})^{2}$, where $\kappa$ and $\xi$ were previously defined in (35) and (46), respectively. Finally, iterating the procedure yields

The only remaining step is to prove the auxiliary lemma.

Proof of Lemma 10: By the RSC condition (33a) and the assumption (45), we have

Furthermore, by convexity of $g$, we have

and the first-order optimality condition for $\widehat{\beta}$ gives

Applying Lemma 1 to bound the term $\|\widehat{\beta}-\beta^{t}\|_{1}^{2}$ and using the assumption $\frac{ck\tau\log p}{n}\leq\frac{2\alpha-\mu}{4}$ yields the bound (84b). On the other hand, applying Lemma 1 directly to (89) with $\beta^{t}$ and $\widehat{\beta}$ switched gives

Appendix D Verifying RSC/RSM Conditions

In this Appendix, we provide a proof of Proposition 1, which verifies the RSC (33b) and RSM (34) conditions for GLMs.

Using the notation for GLMs in Section 3.3, we introduce the shorthand $\Delta:=\beta_{1}-\beta_{2}$ and observe that, by the mean value theorem, we have

for some $t_{i}\in$. The $t_{i}$’s are i.i.d. random variables, with each $t_{i}$ depending only on the random vector $x_{i}$.

Proof of bound (39): The proof of this upper bound is relatively straightforward given earlier results (Loh and Wainwright 2013a). From the Taylor series expansion (92) and the boundedness assumption $\|\psi^{\prime\prime}\|_{\infty}\leq\alpha_{u}$, we have

By known results on restricted eigenvalues for ordinary linear regression (cf. Lemma 13 in Loh and Wainwright 2012), we also have

with probability at least $1-c_{1}\exp(-c_{2}n)$. Combining the two inequalities yields the desired result.

Proof of bounds (38b): The proof of the RSC bound is much more involved, and we provide only high-level details here, deferring the bulk of the technical analysis to later in the appendix. We define

where $T$ is a suitably chosen constant depending only on $\lambda_{\min}(\Sigma)$ and the sub-Gaussian parameter $\sigma_{x}$. (In particular, see (98) below, and take $T=3\tau$.) The core of the proof is based on the following lemma, proved in Section D.2:

With probability at least $1-c_{1}\exp(-c_{2}n)$, we have

Taking Lemma 11 as given, we now complete the proof of the RSC condition (38b). By the arithmetic mean-geometric mean inequality, we have

where $\Delta:=\beta_{1}-\beta_{2}$. Rearranging yields

D.2 Proof of Lemma 11

For a truncation level $\tau^{\prime}>0$ to be chosen, define the functions

By construction, $\varphi_{\tau^{\prime}}$ is $\tau^{\prime}$-Lipschitz and

In addition, we define the trapezoidal function

Taking $T\geq 3\tau^{\prime}$ so that $T\geq\tau^{\prime}\|\Delta\|_{2}$ (since $\|\Delta\|_{2}\leq 3$ by assumption), and defining

where the first equality is the expansion (92) and the second inequality uses the bound (96).

By construction of $\varphi$, each summand in the expression for $Z(\delta)$ is sandwiched as

Consequently, applying the bounded differences inequality yields

Furthermore, by Lemmas 12 and 13 in Appendix E, we have

where the $g_{i}$’s are i.i.d. standard Gaussians. Conditioned on $\{x_{i}\}_{i=1}^{n}$, define the Gaussian processes

and note that for pairs $(\beta,\Delta)$ and $(\widetilde{\beta},\widetilde{\Delta})$, we have

since $\varphi_{\tau^{\prime}\|\Delta\|_{2}}\leq\frac{\tau^{{}^{\prime}2}\|\Delta\|_{2}^{2}}{4}$ and $\gamma_{T}$ is $\frac{2}{T}$-Lipschitz. Similarly, using the homogeneity property

and the fact that $\varphi_{\tau^{\prime}\|\Delta\|_{2}}$ is $\tau^{\prime}\|\Delta\|_{2}$-Lipschitz, we have

where the $\widehat{g}_{i}$’s and $\widetilde{g}_{i}$’s are independent standard Gaussians, it follows that

Applying Lemma 14 in Appendix E, we then have

Note further (cf. p.77 of Ledoux and Talagrand 1991) that

by Lemma 16 in Appendix E. Combining (101), (102), (103), (104), and (D.2), we then obtain

Finally, combining (99), (100), and (106), we see that under the scaling $R\sqrt{\frac{\log p}{n}}\precsim 1$, we have

It remains to extend this bound to one that is uniform in the ratio $\frac{\|\Delta\|_{1}}{\|\Delta\|_{2}}$, which we do via a peeling argument (Alexander 1987; van de Geer 2000). Consider the inequality

Since $\frac{\|\Delta\|_{1}}{\|\Delta\|_{2}}\geq 1$, we have

over the region of interest. For each integer $m\geq 1$, define the set

where $\mu=c^{\prime}\tau^{\prime}\sigma_{x}\sqrt{\frac{\log p}{n}}$. By a union bound, we then have

where the index $m$ ranges up to $M:=\Big\lceil\log\left(c\sqrt{\frac{n}{\log p}}\right)\Big\rceil$ over the relevant region (111). By the definition (109) of $f$, we have

where inequality (i) applies the tail bound (110). It follows that

Multiplying through by $\|\Delta\|_{2}^{2}$ then yields the desired result.

Appendix E Auxiliary Results

In this section, we provide some auxiliary results that are useful for our proofs. The first lemma concerns symmetrization and desymmetrization of empirical processes via Rademacher random variables:

Let $\{Z_{i}\}_{i=1}^{n}$ be independent zero-mean stochastic processes. Then

We also have a useful lemma that bounds the Gaussian complexity in terms of the Rademacher complexity:

Let $Z_{1},\dots,Z_{n}$ be independent stochastic processes. Then

where the $\epsilon_{i}$’s are Rademacher variables and the $g_{i}$’s are standard normal.

We next state a version of the Sudakov-Fernique comparison inequality:

Given a countable index set $T$, let $\{X(t),t\in T\}$ and $\{Y(t),t\in T\}$ be centered Gaussian processes such that

We also have a lemma about maxima of products of sub-Gaussian variables:

Conditioned on $\{X_{i}\}_{i=1}^{n}$, for each $j=1,\ldots,p$, the variable $\left|\frac{1}{n}\sum_{i=1}^{n}g_{i}X_{ij}\right|$ is zero-mean and sub-Gaussian with parameter bounded by $\frac{\sigma_{x}}{n}\,\sqrt{\sum_{i=1}^{n}X_{ij}^{2}}$. Hence, by Lemma 15, we have

Now fix some $t\geq\sqrt{2\sigma_{x}^{2}}$. Since the $\{Z_{j}\}_{j=1}^{p}$ are all nonnegative, we have

where the final inequality follows from the bound (113) with $u=s^{2}-2\sigma_{x}^{2}$, valid as long as $s^{2}\geq t^{2}\geq 2\sigma_{x}^{2}$. Integrating, we have the bound

In order to handle the case when $\rho_{\lambda}$ has points where neither a gradient nor subderivative exists, we assume the existence of a function $\widetilde{\rho}_{\lambda}$ (possibly defined according to the particular local optimum $\widetilde{\beta}$ of interest), such that the following conditions hold:

The function $\widetilde{\rho}_{\lambda}$ is differentiable/subdifferentiable everywhere, and $\|\nabla\widetilde{\rho}_{\lambda}(\widetilde{\beta})\|_{\infty}\leq\lambda L$.

The equality $\widetilde{\rho}_{\lambda}(\widetilde{\beta})=\rho_{\lambda}(\widetilde{\beta})$ holds.

There exists $\mu_{1}\geq 0$ such that $\widetilde{\rho}_{\lambda}(\beta)+\frac{\mu_{1}}{2}\|\beta\|_{2}^{2}$ is convex.

For some index set $A$ with $|A|\leq k$ and some parameter $\mu_{2}\geq 0$, we have

In addition, we assume conditions (i)–(iii) of Assumption 1 in Section 2.2 above.

Under the conditions of Assumption 2, we have the following variant of Theorems 1 and 2:

Suppose $\mathcal{L}_{n}$ satisfies the RSC conditions (4b), and the functions $\rho_{\lambda}$ and $\widetilde{\rho}_{\lambda}$ satisfy Assumption 1 and Assumption 2, respectively. Suppose $\lambda$ is chosen according to the bound (6) and $n\geq\frac{16R^{2}\max(\tau_{1}^{2},\tau_{2}^{2})}{\alpha_{2}^{2}}\log p$. Then for any stationary point $\widetilde{\beta}$ of the program (1), we have

The proof is essentially the same as the proofs of Theorems 1 and 2, so we only mention a few key modifications here. First note that any local minimum $\widetilde{\beta}$ of the program (1) is a local minimum of $\mathcal{L}_{n}+\widetilde{\rho}_{\lambda}$, since

locally for all $\beta$ in the constraint set, where the first inequality comes from the fact that $\widetilde{\beta}$ is a local minimum of $\mathcal{L}_{n}+\rho_{\lambda}$, and the second inequality holds because $\widetilde{\rho}_{\lambda}$ upper-bounds $\rho_{\lambda}$. Hence, the first-order condition (5) still holds with $\rho_{\lambda}$ replaced by $\widetilde{\rho}_{\lambda}$. Consequently, (20) holds, as well.

Next, note that (22) holds as before, with $\rho_{\lambda}$ replaced by $\widetilde{\rho}_{\lambda}$ and $\mu$ replaced by $\mu_{1}$. By condition (v) on $\widetilde{\rho}_{\lambda}$, we then have () with $\mu$ replaced by $\mu_{1}+\mu_{2}$. The remainder of the proof is exactly as before. ∎

For a fixed local optimum $\widetilde{\beta}$, note that we have $\widetilde{\rho}_{\lambda}(\beta)=\sum_{j\in T}\lambda|\widetilde{\beta}_{j}|+\sum_{j\in T^{c}}\frac{\lambda^{2}c}{2}$, where $T:=\left\{j\,\mid|\widetilde{\beta}_{j}|\leq\frac{\lambda c}{2}\right\}$. Clearly, $\widetilde{\rho}_{\lambda}$ is a convex upper bound on $\rho_{\lambda}$, with $\widetilde{\rho}_{\lambda}(\widetilde{\beta})=\rho_{\lambda}(\widetilde{\beta})$. Furthermore, by the convexity of $\widetilde{\rho}_{\lambda}$, we have

using decomposability of $\widetilde{\rho}$. For $j\in T$, we have

whereas for $j\notin T$, we have $\widetilde{\rho}_{\lambda}(\beta^{*}_{j})-\widetilde{\rho}_{\lambda}(\widetilde{\beta}_{j})=0\leq\lambda|\widetilde{\nu}_{j}|$. Combined with the bound (115), we obtain

which is condition (v) of Assumption 2 on $\widetilde{\rho}_{\lambda}$ with $L=1$, $A=S$, and $\mu_{2}=\frac{1}{c}$. The remaining conditions are easy to verify (see also Zhang and Zhang 2012). ∎

References