HIGH-DIMENSIONAL MINIMAX CONCAVE PENALIZED GENERALIZED LINEAR MODELS
Year: 2019

Abstract:
Global solutions to high-dimensional sparse estimation problems with a folded concave penalty (FCP) have been shown to be statistically desirable but are strongly NP-hard to compute which implies the non-existence of pseudopolynomial time global optimization schemes in the worst case. This paper shows that, with high probability, a global solution to generalized linear models with minimax concave penalty (MCP), a specific type of FCP, coincides with a stationary point characterized by the significant subspace second order necessary conditions (SONC). Given that the desired SONC solution admits a fully polynomialtime randomized approximation scheme (FPRAS), we are able to demonstrate the existence of an FPRAS for this strongly NP-hard problem. We further demonstrate two versions of the FPRAS for generating the desired SONC solutions. One follows the paradigm of an interior point trust region algorithm and the other is the well-studied local linear approximation (LLA). Our analysis thus provides new techniques for global optimization of certain NP-Hard problems and new insights on the effectiveness of LLA.

heading: 1 INTRODUCTION
text: This paper concerns global optimization of a folded concave penalized formulation of highdimensional learning generalized linear models, which belongs to statistical/machine learning problems such that the number of dimensions (or number of fitting parameters) p is (much) larger than the number of samples n. This type of problem has recently become very common in a variety of engineering and scientific applications (Fan et al., 2014a; Fan & Li, 2006) including computational biology, speech recognition and image processing (Huang et al., 2016; Algamal & Lee, 2015; Xing et al., 2001; Palmer & Ostendorf, 2001; Bagaria et al., 2012; Paul et al., 2013). Globally minimal solutions to such a nonconvex learning formulation have been shown effective to guarantee desirable statistical performance in order to address high dimensionality (Zhang & Zhang, 2012). Nonetheless, generating a global solution admits no pseudo polynomial-time algorithm, unless “P = NP”; Indeed, global optimality is shown strongly NP-hard to achieve by (Chen et al., 2017; Huo & Chen, 2010) while Chen et al. (2014) show similar results for several related problems in regularized minimization. In contrast to the existing pessimistic result, we derive herein a fully polynomial-time randomized approximation scheme (FPRAS, as defined in 4) that theoretically ensures global minimality to 2 with high probability.
Specifically, consider a high-dimensional generalized linear model (GLM) as follows. Let X = (x1, ..., xn)
ᵀ be the n×p design matrix with xi = (xi1, ..., xip)ᵀ, i = 1, ..., n, and Y = (y1, ..., yn)ᵀ be the n-dimensional response vector. We will assume the design matrix X is fixed, while the mean of the response is given by E[yi] = ψ′(x>i β
true) for some known link function ψ : Θ→ <, where Θ ⊆ < and βtrue = (βtrue1 , ..., βtruep ) is the unknown vector of true parameters of the model. Such a setup can be seen as a generalization of linear regression models with the link function allowing for nonlinear transformations that enable a more flexible approach to model estimation. The high-dimensional regression problem is to estimate βtrue given knowledge of X , Y , and ψ in the undesirable scenario where p n > 0. To that end, traditional statistical learning schemes
often resort to the following formulation:
L(β) = n∑ i=1 `(yi, xi, β) = 1 n n∑ i=1 [ψ(xᵀi β)− yix ᵀ i β], (1)
which, according to traditional statistical theories, would result in over fitting in general under the high-dimensional setting.
To resolve over fitting, modern statistical theories favor a modified formulation as below:
min β Q(β) := L(β) + p∑ j=1 Pλ(|βj |)  , (2) where Pλ(| · |) is sparsity-inducing regularization term that penalizes any nonzero dimensions in the minimizer, and λ > 0 is a tuning parameter. Under the assumption that the true regression parameter βtrue is sparse, a global optimizer to equation 2 has been shown effective to address over fitting for many choices of specific regularization functions Pλ. Indeed, one of the most successful choices of Pλ is the much studied Lasso-based regularized (Tibshirani, 1996), aka, the `1(-norm) penalty, which has been demonstrated to entail desirable statistical properties (Bickel et al., 2009; Negahban et al., 2012). Another admirable property of the Lasso is that, especially when applied to least squares linear regression, it yields an extremely tractable problem via a variety of algorithms (Friedman et al., 2008; 2010). However, per (Zhao & Yu, 2006; Fan & Li, 2001), Lasso is not selection consistent without a strong irrepresentable condition and may sometimes introduce nontrivial estimation bias.
As a popular alternative to Lasso, the folded concave penalty (FCP) was first introduced by (Fan & Li, 2001). There are mainstream examples of FCP functions, including the SCAD by (Fan & Li, 2001) and MCP by (Zhang, 2010). This paper will focus on MCP, defined as Pλ(|t|) =∫ |t|
0 (aλ−s)+
a ds for some fixed parameter a > 0. In contrast to the Lasso, the FCP regulaizations achieve variable selection consistency non-contingent on the irrepresentable condition and is demonstrated to be unbiased (Fan & Li, 2001). Furthermore, (Zhang & Zhang, 2012) demonstrated that the global solution to the FCP-regularized formulation leads to desirable recovery of the oracle solution.
Nonetheless, FCP problems are significantly harder to solve than Lasso, the new penalty term moves the problem outside the realm of convex optimization, (Chen et al., 2017) even showed that any estimation problem with convex loss and folded concave regularization to be strongly NP-hard, ruling out the possibility of a pseudo-polynomial-time global optimization algorithm. (Liu et al., 2016) were seemingly the first to propose a global approach to the problem called MIPGO which reformulates the problem into a mixed integer program. Yet, the worst-case complexity of MIPGO is in exponential time.
Perhaps for this reason, current literature tends to focus on local algorithms for the FCP-regularized learning problems. The local quadratic approximation algorithm by (Fan & Li, 2001) is an example of a majorization minimization algorithm, an approach which is also related to the local linear approximation (LLA) algorithm proposed by (Zou & Li, 2008). LLA was further explored by (Fan et al., 2014b) showing the oracle property can be obtained with high probability despite the local approach. (Mazumder et al., 2011; Fan & Lv, 2011) demonstrate coordinate optimization approaches for FCP while (Wang et al., 2014) used an approximate regularization path-following algorithm to obtain the optimal convergence rate to statistically desirable local solution. (Wang et al., 2013) analyzed the CCCP algorithm and showed under certain conditions that it asymptotically finds the oracle estimator. (Liu et al., 2017) took an algorithm agnostic approach by analyzing local solutions satisfying second order KKT conditions and showed desirable statistical properties like recovering the oracle solution and sparisty. These results discussed above primarily relate to FCP-regularized linear regression, a special case of GLM where ψ is specifically the identity function. For analyses which encompass GLM’s with FCP regularizers, (Fan & Lv, 2011) showed that GLM’s, even in ultra high dimensional variable selection problems, have oracle properties when using FCP regularization and demonstrated a coordinate optimization algorithm for finding local solutions. In the area of Mestimators, which is a further generalization of our estimation method beyond even GLMs, (Loh & Wainwright, 2013; Loh et al., 2017a) showed that under certain conditions all local solutions must be within statistical precision of the true parameter and its support while (Loh et al., 2017b) demonstrate a two-step algorithm involving composite gradient descent to find a local solution. Bian &
Chen (2014) demonstrate a optimality conditions for a class of nonconvex optimization problems using nonlipscitz regularization.
From the numerous results pertaining to local solution schemes above, our research question is why local solutions are repetitively successful. In other words, are there certain geometric properties of the learning formulation equation 2 with FCP that allow all local schemes to perform well independent of the specific designs of the algorithmic procedures? Our answer to this question is affirmative; we show herein that all local solutions within an efficiently achievable sub-level set are actually globally optimal. Those local solutions are characterized by the significant subspace second-order necessary conditions (S3ONC) provable admit FPRAS’s. The S3ONC are weaker conditions than the standard second-order KKT conditions. As per this result, all S3ONC-guaranteeing algorithms (which include a large spectrum of local algorithms) belong to the class of FPRAS’s for global optimization of the strongly NP-hard FCP-based learning problem. We subsequently develop theories for two specific algorithms of this type: one gradient-based method and the other is the same as the LLA.
It is worth noting that (Zhang, 2010) provides conditions to establish the uniqueness of local solutions to FCP-based learning. When local solutions are unique, then any local optimization algorithms would ensure global optimality. However, a few critical assumptions are necessary to achieve the uniqueness result and, furthermore, many report numerical experiments, e.g., those in (Fan et al., 2014b; Liu et al., 2017; 2016; Fan & Li, 2001) indicate the non-uniqueness of local solutions, instead. In contrast, our results in this paper imposes only standard assumptions commonly shared by a flexible set of high-dimensional GLMs and are applicable even if the local solutions are non-unique. To our knowledge, this is the first geometric proof that global solutions coincide with pseudo-polynomial-time computable local solutions in an FCP-based regression formulation with high probability. The resulting algorithms are the first few FPRAS’s to this problem.
Two works with notable relations to our own are Liu et al. (2017) and Liu & Ye (2019). The first applies a similar analytical framework to linear regression problems however our generalization to GLMs adds significant flexibility and it was unknown for their result that the oracle solution implies global optimality since it was only as of Zhang & Zhang (2012) that global optimality was known to potentially imply the oracle solution. Further, it is nontrivial to extend their existing result to global optimal results. (Liu & Ye, 2019) on the other hand is a more general setup than our own though the tradeoff is that our rate is better and S3ONC solutions do not ensure global optimality to the in-sample training error for their setup.
The rest of this paper is organized as follows. Section 2 goes through specific problem assumptions and explains the S3ONC. Section 3 contains our main result for global optimality and uses it to make additional claims for LLA. Section 4 contains numerical results to verify our theoretical findings. Section 5 provides concluding remarks.
In this paper we will use ‖·‖0 to denote the number of nonzero entries, | · | to denote the `1-norm if the argument is a vector, or cardinality if the argument is a set, ‖·‖ to denote the `2-norm, ‖·‖max to denote the maximum norm and ‖·‖min to denote the absolute value of the entry with the smallest magnitude. (·)+ is used equivalently to max(0, ·). For any vector v, vQ is intended as (vj : j ∈ Q). For any set Q, we denote the complement as Qc. In particular, let S be the true support set, that is, S := {j : βtruej 6= 0} and its complement is Sc. We occasionally use the term the “oracle solution” to refer to the solution βoracle defined as βoracle ∈ arg min
β: βj=0, ∀j /∈S L(β). The oracle solution
is a hypothetical solution which assumes prior knowledge of the true support set S and thus can be considered a theoretical benchmark.

heading: 2 SETUPS, PRELIMINARIES, AND ASSUMPTIONS

heading: 2.1 SETUPS AND ASSUMPTIONS
text: Our analysis will focus on sparse GLMs that have a fixed design matrix and satisfy the following assumptions:
(A1) Assume that
(i) bu ≥ ψ′′(xᵀi β) ≥ bl > 0 for all x ᵀ i β ∈ Θ;
(ii) There exists K > 0 such that the design matrix satisfies 1n ‖Xj‖ 2 < K for all j ∈ [p].
Let the tuning parameter a in Pλ satisfy K < (bua)−1. (A2) The vector of residuals W ∈ <n such that W := y − E[y|X] is subgaussian(σ) which
means it satisfies that P [|〈W, v〉| ≥ t] ≤ 2exp(−t2/2σ2), for all v : ‖v‖ = 1 and any t > 0;
(A3) There exists a sequence {rd ≥ 0 : d = 1, 2, ..., p} such that the following are satisfied: (i) For any d1, d2 : 1 ≤ d1 ≤ d2 ≤ p, we have rd1 ≥ rd2 ; (ii) There exists some p̃∗ : 2|S| ≤ p̃∗ ≤ p such that rp̃∗ > 0; (iii) For all d : 1 ≤ d ≤ p, it holds that n−1 ‖Xβ‖2 ≥ rd ‖β‖2 for any β ∈ <p : ‖β‖0 ≤ d. Remark 1. Part (i) of (A1) states that our link function is both strongly convex and continuously differentiable; that is, the gradient being Lipschitz continuous. Many types of traditional GLM problems satisfy this constraint including those for normal (linear regression), categorical (logistic regression), binomial, gamma and Poisson distributions, although in some cases the mild assumption on the boundedness of Θ has to be made. Even though the original domain of the link function can be unbounded, one may still consider its bounded subset given that it contains the vector of true parameters. Part (ii) of (A1) can be assumed without loss of generality by normalizing the design matrix columns. Remark 2. (A2) is a common assumption in the literature, for example (Negahban et al., 2012) (Wang et al., 2013). Remark 3. Both (i) and (A2) are satisfied by a number of GLM setups, one example is linear regression. In such a setup, the response Y takes a gaussian distribution, while the gradient of the link function (encoding E[Y |X]) is the identity. Note, it is difficult to have one without the other, since the loss formulation 1 is simply a log-likelihood maximization applied to a distribution within the class of exponential dispersion models (Jørgensen, 1987). Another classic example is logistic regression which is used for a Bernoulli or binomial distributed Y along with a logit link function. It should be noted that we treat the matrix X as fixed, so its generative distribution is not important to the analysis outside of whether it satisfies the assumptions and constraints mentioned. In the numerical experiments in Section 4, we use i.i.d. gaussian generation method since, as discussed for Definition 1, it means our design matrices will satisfy (A3) with high probability. Remark 4. Assumption (A3) can be understood to be a lower bound on the eigenvalues for principal sub-matrices of XᵀX of dimension d × d for all d ∈ [p]. For every d : d ≤ p̃∗, the lower bounds are positive, meaning that the smallest eigenvalues of the d× d principal sub-matrices are assumed positive.
According to (Liu et al., 2017), Assumption (A3), for certain parameters, is provably a weaker condition than the restricted eigenvalue (RE) condition, as defined in Definition 1 below and first introduced by (Bickel et al., 2009) as a plausible assumption to allow for the desired recovery quality of Lasso. The RE condition is a common assumption in the high-dimensional learning literature, such as (Zhang et al., 2014) and (Fan et al., 2014b). Definition 1. (RE condition (Zhang et al., 2014)) The matrix X ∈ <n×p is said to satisfy the RE condition if, for some re > 0, it holds that 1n ‖Xδ‖ 2 ≥ re ‖δ‖2 for all δ ∈ ⋃ |Ŝ|=s C(Ŝ) where C(Ŝ) := {δ := (δi) ∈ <p : |δŜc | ≤ 3|δŜ |}, δŜc := (δj : j ∈ Ŝc), and δŜ := (δj : j ∈ Ŝ). Furthermore, the largest possible re is said to be the restricted eigenvalue constant of X .
Random design matrices with with i.i.d. rows generated following subgaussian distributions as in (A2) have been shown to satisfy the RE condition with high probability (Zhou, 2009) while proposition 1 in Loh & Wainwright (2013) includes a proof that with high probability, restricted strong convexity (RSC) is satisfied for a setup equivalent to our own. Note that satisfaction RSC implies the RE condition above. Thus, within our setup, (A3) is also satisfied with high probability for our setting.
2.2 PRELIMINARIES ON S3ONC
Our results focus on the S3ONC solutions, which has been formerly introduced by (Liu et al., 2017) in the special case of high-dimensional linear regression as a relaxation of the standard second-order
KKT conditions. The definition of S3ONC depends on the notion of first order necessary conditions (FONC) as below. Definition 2 (FONC). A solution β∗ satisfies the first order necessary conditions (FONC) if
0 ∈ 1/n n∑ i=1 [ψ′(xᵀi β ∗)− yi]xi + ( P ′λ(|β∗j |)∂(|β∗j |), 1 ≤ j ≤ p ) (3)
where ∂(| · |) denotes the subdifferential of | · |. Definition 3 (S3ONC). A solution β∗ satisfies the significant subspace second-order necessary condition (S3ONC) if it satisfies FONC and for all j ∈ {j : β∗j 6= 0},
∂2Q(β) (∂βj)2 ∣∣∣∣ β=β∗ ≥ 0 (4)
if the second derivative exists. Remark 5. The S3ONC can be intuited as the second order necessary condition applied only to the dimensions where βj 6= 0, i.e., the significant dimensions. Since the S3ONC is weaker than the standard second-order KKT conditions, any algorithm that guarantees the second-order KKT conditions can be used to obtain an S3ONC solution, by requiring a more stringent optimality condition, may be slower than necessary. One specifically S3ONC guaranteeing approach, presented in (Liu & Ye, 2019), utilizes an interior point trust region algorithm in order to guarantee an S3ONC solution in polynomial time. This is the scheme which will be used later in Section 4.

heading: 3 MAIN RESULTS
text: We now present our theoretical results for global optimization of FCP penalized GLMs. All proofs can be found in the appendix. We will make use of a short-hand notation:
βLasso ∈ arg minL(β) + λ|β|. (5) Theorem 1. Suppose assumptions (A1), (A2), and (A3) with any p̃∗ : p̃∗ ≥ 2|S|. Let β∗ be an arbitrary S3ONC solution to equation 2 with Pλ specified as the MCP. Assume that Q(β∗) ≤ Q(βtrue) + Γ for an arbitrary Γ ≥ 0. (i) Let the sub-optimality gap satisfies Γ < Pλ(aλ) − σ2 bln ( p̃∗ + 2 √ p̃∗t+ 2t ) ; (ii) choose Pλ(aλ) > σ 2 2nbl (1 + 2 √ t′ + 2t′) + σ2 n |S|(1+2 √ t′+2t′)+Γbl
bl(p̃∗−2|S|+1) and (iii) assume that the minimal signal strength satisfy
∥∥βtrueS ∥∥min > √ 8σ2
rp̃b2l n
( p̃∗ + 2 √ p̃∗t+ 2t ) + 8
rp̃bl min
{ λ2
rp̃ |S|, Pλ(aλ)|S|+ Γ } then the following two statements hold
(a) β∗ is an oracle solution with probability at least 1− exp(−t+ p̃∗ ln( pep̃∗ ))− exp(−(p̃ ∗ +
1)(t′ − ln p)) · 1−exp(−(p−p̃ ∗)(t′−ln p))
1−exp(−t′+ln p) .
(b) β∗ is both an oracle solution and an globally optimal solution to equation 2 with probability 1− 2 exp(−t+ p̃∗ ln( pep̃∗ ))− 2 exp(−(p̃ ∗ + 1)(t′ − ln p)) · 1−exp(−(p−p̃ ∗)(t′−ln p)) 1−exp(−t′+ln p) .
Remark 6. Theorem 1 (especially in the second statement) is perhaps the first result that establishes a set of conditions for any S3ONC solution to be globally optimal with high probability. Further, this result is algorithm independent which allows for greater flexibility compared to most existing results as in (Loh et al., 2017b) and (Fan & Lv, 2011) which rely on a specific algorithm choice. Remark 7. The second part follows quite easily from the first due to the uniqueness of βoracle as well as the fact that βopt must also be an S3ONC solution. Thus by applying the first part of the Theorem to βopt we are able to show that both our arbitrary β∗ and βopt coincide with the unique βoracle.
Remark 8. The above constraints on Γ, Pλ(aλ) and ‖βtrueS ‖min may initially seem disparate but can all be converted to constraints on the sample size n as is shown in Corollary 1 below. This is possible because Γ can be bounded by some function of n−γ for some γ > 0. Given that, it can be seen that the lesser side of inequalities (i),(ii) and (iii) go to 0 as n grows. Further discussion of how this is achieved for Corollary 1 can be found in Remark 12. Corollary 1. Assume ln p ≥ 1, bl ≤ 1, and s ≥ 1. Let β∗ be an S3ONC solution to equation 2. Let assumptions (A1), (A2), and the RE condition as defined in Definition 1 hold. Assume that Q(β∗) ≤ Q(βLasso) almost surely, where βLasso is the optimal solution to the Lasso problem with penalty coefficient λLasso = σ √ lnp
n1−γ/2 where γ ∈ [0, 1] is an arbitrary scalar. Let λ = σre √ ln p nγ/2
and a ∈ [.8, 1). There exist problem independent constants C1 > 0,C2 > 0 and C2 > 0 such that if
n > max C1bl , [ C2 s bl ] 2 1−γ , [ C3
sσ2 ln p
‖βtrueS ‖ 2 min b 2 l r 4 e ]2/γ (6) Then β∗ is the global solution to 2 with probability at least 1 − C4 exp(−C5snγ/2 ln p) − C6 exp(−C7bunγ/2 ln(p)) for problem independent constants C4, C5, C6 and C7 Remark 9. Corollary 1 indicates that for γ > 0, the global optimal solution coincides with computable S3ONC solution with overwhelming probability given that the sample size meets certain requirements. It should specifically be noted that the relationship between n and p require only ln p nγ/2
= O(1), which ensures the applicability to the high-dimensional setting even if n p. Remark 10. (Liu & Ye, 2019) has derived a gradient-based algorithm that provably ensures an S3ONC solution at pseudo-polynomial-time complexity. When n is properly large, this pseudopolynomial-time algorithm enables a straightforward design of an FPRAS for generating the global optimal solution as follows.
FPRAS: A pseudo-polynomial-time algorithm that generates global optimum at high probability
Step 1. Initialize the parameters δ, λ, a, â, k = 0 and βLasso by solving equation 5 Step 2. If Case 1: |βkj | ∈ (0, aλ) for some j = 1, ..., p, then choose an arbitrary ι ∈ {j : |βkj ∈
(0, aλ)} and solve
βk+1ι ∈ arg min β [∇L(βk)]ι · β + Pλ(|β|)
s.t. (β − βkι )2 ≤ δ2
and let βk+1j = β k j ∀j 6= ι. Go to Step 3.
Else Case 2: If —βkj /∈ (0, aλ) for all j = 1, ..., p then for all j = 1, ..., p: — If βkj = 0 then β k+1 j = â · [ |[∇L(βk)]j | − λ ] + · sign(−[∇L(βk)]j). — If |βkj | ≥ aλ, then β k+1 j = β k j − â · [∇L(βk)]j . Go to Step 3.
Step 3. Algorithm stops if |βkj | /∈ (0, aλ) and ∥∥βkj − βk+1j ∥∥ < δ. Otherwise, let k := k + 1 and go
to Step 2. Remark 11. Here, the above algorithm has iteration complexity of O (( Q(βLasso) − Q(βopt) ) ·
max{(1/(2a) − bu/2)−1, 2b−1u , (1/a − bu/2)−1} · 1/δ2 )
for any γ-accuracy S3ONC solution. In this iteration complexity, all the quantities are verifiably upper bounded by a polynomial function of dimensionality p and the desired accuracy 1/γ. Furthermore, βLasso is a solution to a convex problem, which be generated within polynomial time and the per–iteration problem admits a closed form, whose complexity is strongly in polynomial-time. Therefor, this algorithm is an FPRAS in generating an S3ONC (global) solution where the term FPRAS is defined in Definition 4. Definition 4. (FPRAS) Let f : L → <+0 be a function representing the optimal solution to a problem. Let A be a probabilistic algorithm, which takes as an input an instance of the problem x ∈ L, and a parameter > 0. We call A a fully polynomial randomized approximation scheme, if it has the following properties
• the algorithm returns a value within the required precision with probability at least 1/2.
P [ |A(x, )− f(x)| ≤ f(x) ] > 1/2
• the running time of A is polynomial in both the size of the input x and the precision 1/
Note that the specific probability bound is fairly arbitrary as long as it is a constant probability bounded away from 1/2 (Vazirani, 2013).
Remark 12. We are able to remove Γ from the result by bounding the performance difference between βtrue and βlasso using similar techniques as in Bickel et al. (2009). In order to use this bound for our S3ONC solution, we require that Q(β∗) ≤ Q(βLasso). However, this can generally be obtained by initializing any S3ONC guaranteeing algorithm with βLasso in a similar fashion to (Fan et al., 2014b) for LLA. The FPRAS above follows the same initialization scheme.
Remark 13. The above specification of values for a, λ and λLasso can be thought of as examples rather than strict requirements. A closer examination of the proof for Corollary 1 will reveal that the values for λ and λLasso can be chosen in a much more flexible fashion, though the corresponding values of C1 through C7 may be different for different combinations of λ and λLasso.
The techniques used in the proof of Theorem 1 can be used to provide insights into other optimization schemes. As an example, we can apply the same analysis to the state-of-the-art FCP-based algorithm, LLA, using the framework in (Fan et al., 2014b) as a starting point.
LLA: local linear approximation. Step 1. Set k = 0. Initialize the algorithm with β0 = βLasso, where βLasso is generated by solving equation 5. Let N be the maximal iteration number. Step 2. For all k = 1, ..., N , solve the following convex program to generate βk+1:
βk+1 ∈ arg minL(β) + ∑ j∈[p] P ′λ(|βkj |) · |βj |,
where P ′λ is the first derivative of Pλ. Let k := k + 1.
We can show that in fact the LLA is another FPRAS that achieves the global optimal solution. The proof of this can be found in the appendix.
Corollary 2. For problem equation 2. If ‖βtrueS ‖min > (a + 1)λ, λ > max{ 3λ Lassos1/2
blre ,
4σ √ s+2 √ st1+2t1
bl(an/bu)1/2 ,
2σ √ s+2 √ st2+2t2
bl √ nre } and the RE condition in Definition 1 holds, the following holds.
(a) The LLA algorithm initialized with βLasso converges to the oracle solution in two iterations with probability 1− φ0 − φ1 − φ2, where
φ0 := P ( ∥∥βLasso − βtrue∥∥ max > λ) ≤ 2p exp(−(λ Lasso)2nbua 8σ2 ),
φ1 := P ( ∥∥∥OScp̃`n(βoracle)∥∥∥max ≥ λ) ≤ (pes )s exp(−t1) + 2 exp(−λ2abun8σ2 ),
φ2 := P ( ∥∥βoracleS ∥∥min ≤ aλ) ≤ (pes )s exp(−t2),
(b) If in addition (A1) and (A2) holds, while the parameters of (a, λ) satisfy that
Pλ(aλ) > σ2
2nbl (1 + 2
√ t4 + 2t4) + σ2 n |S|(1+2 √ t4+2t4)bl
bl(p̃∗−2|S|+1) and the and Pλ(aλ) > σ2 bln ( p̃∗ + 2 √ p̃∗t3 + 2t3 ) and let the minimal signal strength satisfy ‖βtrueS ‖min >√
8σ2
rp̃b2l n
( p̃∗ + 2 √ p̃∗t3 + 2t3 ) + 8rp̃bl min{ λ2 rp̃ |S|, Pλ(aλ)|S|} then the LLA algorithm ini-
tialized by βLasso converges to the global solution in two iterations with probability at least 1− φ0 − φ1 − φ2 − φ3
where
φ3 := P (B oracle 6= Bopt) ≤ exp(−t3 + p̃∗ ln(
pe p̃∗ ))
+ exp(−(p̃∗ + 1)(t4 − ln p)) · 1− exp(−(p− p̃∗)(t4 − ln p))
1− exp(−t4 + ln p) , (7)
and t1, t2, t3, t4 > 0 are arbitrary constants. Remark 14. Since each iteration of the LLA solves a convex program, which can be done within polynomial-time. When n is properly large, the above theorem then indicates that the LLA is another FPRAS in globally optimizing the FCP-based nonconvex formulation.

heading: 4 NUMERICAL EXPERIMENTS

heading: 4.1 EXPERIMENTAL SETUP
text: We focus our tests on sparse logistic regression. Our problem and data are implemented in a similar way as (Fan et al., 2014b). We construct βtrue as below: Firstly, βtrueS is constructed randomly by choosing 10 elements of β and choosing the magnitude of each to be a uniform value within [1, 2]. Each value is chosen to be negative with probability .5. Then, the remaining entries βtrueSc are set to be 0. The design matrix X ∈ <n×p is constructed by generating n iterations of xi ∼ Np(0,Σ) where Σ = (.5|j−j
′|)p×p. We then generate Y using a Bernoulli distribution where P (yi = 1) = (1 + e−x ᵀ i β true
)−1. With this data, we train a logistic regression model by invoking Algorithm 1 in solving equation 2 with MCP for S3ONC solutions initialized with Lasso implemented in Python 3. The tuning parameters λ and a are obtained by cross validation following (Fan et al., 2014b).
We would like to ascertain whether our FCP classifier, obtained using S3ONC methods, is actually the global optimal solution. We do this by taking each element of the FCP classifier and perturbing it to find a new potential solution. Each element’s perturbation is independent and generated by a N(0, 1/p1/2))-random variable. We then check if this perturbed classifier has better FCP regularized performance on the training data than the FCP classifier. If not, we repeat until either a better solution is found, or until 2000 perturbations have been tried.
Additionally, we compare our solution’s statistical performance to those of other popular regularization methods. Using the data generation method above, we obtain two sets of data, both with 100 samples. One set is for training the model, and the other is the test set for out-of-sample tests. We repeat the above process for 100 times to generate 100 training-and-test instances, each with 100 samples. We compare those trained using the method described above with Lasso solutions generated by the global minimizer to equation 5 and an estimator generated by solving equation 2 when Pλ is substantiated by an `2 penalty. The Lasso and `2 classifiers are solved using the scikit learn python library.
We compare the above estimators in terms of statistical performance for both `1 loss: |β∗ − βtrue| and `2 loss: ‖β∗ − βtrue‖.

heading: 4.2 NUMERICAL RESULTS
text: Table 1 contains the numbers from optimality analysis. This technique did not yield a single perturbed solution that could beat the FCP classifier obtained from the FPRAS in any of our thousands of iterations.
Table 2 shows the numerical results for the statistical performance measurements. We show the two performance measures for each of the three classifiers for tphree different problem types.
As expected, the FCP classifier generally outperformed the lasso and `2 classifiers. The margins are fairly thin between FCP and lasso, especially compared to the standard deviation. Other values of n and p were tried but the results generally followed the same pattern.
As a result we tentatively conclude that our numerical results align with our theoretical results though further testing of the global optimality probability would be valuable.

heading: 5 CONCLUSIONS
text: This paper investigates both the theoretical and empirical performance of FPRAS’s on MCP regularized GLMs. Despite such a problem being strongly NP-Hard, we have shown two FPRAS that achieve global optimality. To our knowledge this is the first probability bound for global optimization of FCP regularized GLMs using an FPRAS. Further, the same technique can be used to extend other results in order to obtain global optimization bounds for a wide variety of problems.
Though this paper focuses on GLMs, further exploration will focus on the question whether similar results can be found for more general problem classes under weaker assumptions. High-dimensional M-estimation problems could potentially be a future avenue of investigation.

heading: A APPENDIX
text: The Appendix is organized as below: Section A.1 presents the proofs for the main results, Sections A.2 and A.3 present central lemmata to be useful in Section A.1.
A.1 PROOF OF MAIN RESULTS
A useful relationship in our proofs is that, for an S3ONC solution β∗ within {β∗ : Q(β∗) ≤ Q(βtrue) + Γ} for any Γ ≥ 0, we have the following useful inequality under Assumption (A1):
bl 2n ‖Xδ∗‖2 − 1 n W ᵀXδ∗ + ∑ j∈S Pλ(|β∗j |) ≤ ∑ j∈S Pλ(|βtruej |) + Γ, (8)
where δ∗ = β∗ − βtrue. This is obtained by invoking the strong convexity of ψ, which leads to ψ(xᵀi β ∗) ≥ ψ(xᵀi βtrue) + ψ′(x ᵀ i β true)(xᵀi β ∗ − xᵀi βtrue) + 0.5 · bl(x ᵀ i β ∗ − xᵀi βtrue)2.
Proof of Theorem 1. First, given our assumption that (A1) holds, that (i) p̃∗ ≥ 2|S|, (ii) β∗ is S3ONC satisfying Q(β∗) ≤ Q(βtrue) + Γ for some Γ ≥ 0, and (iii) Pλ(aλ) > σ 2
2nbl (1 + 2
√ t′ +
2t′)+ σ2
n |S|(1+2 √ t′+2t′)+Γbl
bl(p̃∗+1−2|S|) , we can apply Lemma 5 with p̃ = p̃ ∗. This means that ‖β∗ − βtrue‖0 ≤
p̃∗ with probability at least 1− exp(−(p̃∗+ 1)(t′− ln p)) · 1−exp(−(p−p̃ ∗)(t′−ln p))
1−exp(−t′+ln p) . From this, given the additional assumption that (A3) holds, we can apply the second part of Lemma 4 with p̃ = p̃∗ to get that for any t > 0, 1n ‖X(β ∗ − βtrue)‖2 ≤ 8σ 2
b2l n
( p̃∗ + 2 √ p̃∗t+ 2t ) + 8bl min{λ
2(|S| − ‖β∗‖0)r −1 p̃∗ , Pλ(aλ)(̇|S| − ‖β∗‖0) + Γ} holds with probability at least 1 − exp(−t + p̃∗ ln( pe p̃∗ )). Given that for 2 arbitrary sets A and B
P (A ∩B) = P (B)P (A|B) = (1− P (Bc))(1− P (Ac|B)) = 1− P (Ac|B)− P (Bc) + P (Bc)P (Ac|B) = 1− P (Ac|B)− P (Bc)(1− P (Ac|B)) ≥ 1− P (Ac|B)− P (Bc)
(9)
Therefor they hold simultaneous with probability at least 1 − exp(−t + p̃ ln( pep̃∗ )) − exp(−(p̃ ∗ + 1)(t′ − ln p)) · 1−exp(−(p−p̃ ∗)(t′−ln p))
1−exp(−t′+ln p) .
The same sequence of arguments can be used to show that βopt also satisfies ‖βopt − βtrue‖ ≤ p̃∗ and 1n ‖X(β opt − βtrue)‖2 ≤ 8σ 2
b2l n
( p̃∗ + 2 √ p̃∗t+ 2t ) + 8bl min{λ 2(|S|−‖β∗‖0)r −1 p̃∗ , Pλ(aλ)(̇|S|−
‖β∗‖0)+Γ}with the same probability. Using again the union bound and DeMorgan’s law, we say β∗ and βopt satisfy the above conditions simultaneously with probability 1− 2 exp(−t+ p̃∗ ln( pep̃∗ ))− 2 exp(−(p̃∗ + 1)(t′ − ln p)) · 1−exp(−(p−p̃
∗)(t′−ln p)) 1−exp(−t′+ln p) . With this, our Γ assumption and our minimal
signal strength assumption, we can apply Lemma 6 to show that β∗ = βopt with probability at least

heading: 1− 2 exp(−t+ p̃∗ ln( pep̃∗ ))− 2 exp(−(p̃
text: ∗ + 1)(t′ − ln p)) · 1−exp(−(p−p̃ ∗)(t′−ln p)) 1−exp(−t′+ln p) .
Proof of Corollary 1. First we need to bound Γ. In order to do this we use the lasso problem Qlasso(β) = ∑ i∈N `(β, xi, yi) + ∑ j∈P λlasso|βj | as well as the concavity of MCP over positive values to get the following 2 inequalities
Qlasso(βlasso) ≤ Qlasso(βtrue)∑ i∈N `(βlassoj , xi, yi)− `(βtruej , xi, yi) ≤ ∑ j∈P λlasso(|βtruej | − |βlassoj |)
≤ ∑ j∈P λlasso|βlassoj − βtruej |
(10)
and
∑ j∈P Pλ(β true j )− ∑ j∈P Pλ(β lasso j ) ≤ ∑ j∈P P ′λ(β lasso j )(|βtruej | − |βlassoj |) ≤ ∑ j∈P λ|βlassoj − βtruej |
(11)
We also need 2 results from the proof for φ0 in Corollary 2 which shows that both |δ`Sc | ≤ 3|δ ` S | and bln ∥∥Xδ`∥∥2 ≤ 3λlasso|δ`S | conditional on A where δ` = βlasso − βtrue. Given our restricted eigenvalue assumption ‖Xδ `‖2
n‖δ`‖2 ≥ re, this can be used to show
|δ`| ≤ 4|δ`S | ≤ 4 √ s ∥∥δ`S∥∥2∥∥δ`S∥∥ ≤ 4√s ∥∥δ`∥∥2∥∥δ`S∥∥
≤ 4 √ s
ren ∥∥Xδ`∥∥2∥∥δ`S∥∥ ≤ 4 √ s re 3λlasso|δ`S | bl ∥∥δ`S∥∥ ≤ 4 √ s re 3λlasso √ s ∥∥δ`S∥∥ bl ∥∥δ`S∥∥
(12)
which means |δ`| ≤ 12λ lassos blre with conditional on A which occurs with probability at least 1 − 2p exp(−(λ
lasso)2nbua 8σ2 )
Finally we are able to bound gamma by combining the above
Γ ≤ Q(β∗)−Q(βtrue) ≤ Q(βlasso)−Q(βtrue) (13) ≤ ∑ i∈N `(βlassoj , xi, yi)− ∑ j∈P Pλ(β lasso j )− [ ∑ i∈N `(βtruej , xi, yi)− ∑ j∈P Pλ(β true j )] (14)
≤ ∑ j∈P (λlasso|βlassoj − βtruej |+ λ|βlassoj − βtruej |) (15)
≤ (λlasso + λ)|δ`| ≤ (λlasso + λ)12λ lassos
blre (16)
Next consider the conditions necessary to apply Theorem 1. We have assumptions (A1) and (A2) and (A3) per our assumption that the RE condition holds combined with 7. That leaves the 3 requirements on Γ, Pλ(aλ) and ‖βtrueS ‖min. We will convert each of these to inequalities on n
Utilizing 16 and substituting λ = Qσre √ ln p nγ/2 and λlasso = σ √ ln p n1−γ/2 , whereQ, > 0 are arbitrary constants, and setting p̃∗ = 4s, t = p̃∗nγ/2 ln p, t′ = nγ/2 ln p we get the following
Pλ(aλ) > σ2
2nbl (1 + 2
√ t′ + 2t′) + σ2
2 s(1 + 2 √ t′ + 2t′) + Γbl
bl(p̃∗ − 2s+ 1) (17)
n > 8 + 12 2 + 12 Q
blaQ2 = C1/bl (18)
Γ < Pλ(aλ)− σ2
bln
( p̃∗ + 2 √ p̃∗t+ 2t ) (19)
n >
[( 12
aQ +
√ 20 + 12 2
aQ2
) s
bl
] 2 1−γ
= [ C2 s
bl
] 2 1−γ
(20)
∥∥βtrueS ∥∥min > √ 8σ2
rp̃b2l n
( p̃∗ + 2 √ p̃∗t+ 2t ) + 8
rp̃bl min{λ
2
rp̃ |S|, Pλ(aλ)|S|+ Γ} (21)
n > [ (160 + 8Q2)
sσ2 ln p
(‖βtrueS ‖min r4sblre)2
]2/γ = [ C3
sσ2 ln p
(‖βtrueS ‖min r4sblre)2
]2/γ (22)
For some constants C1, C2 and C3
We can then apply Theorem 1 (conditional onA) substitute our values and simplify to get that β∗ is the global solution with probability at least
1− 2 exp(−t+ p̃∗ ln(pe p̃∗
))− 2 exp(−(p̃∗ + 1)(t′ − ln p)) · [
1− exp(−(p− p̃∗)(t′ − ln p)) 1− exp(−t′ + ln p)
]
≥1− 2 exp(−(nγ/2 − 1)4s ln p)− 2 p−p̃∗∑ k=1 exp(−(p̃∗ + k)(nγ/2 − 1) ln p)  ≥1− 2 exp(−(nγ/2 − 1)4s ln p)− 2 exp(−[(4s+ 1)(nγ/2 − 1)− 1] ln p) ≥1− C4 exp(−C5snγ/2 ln p)
(23)
We then use the same technique as in Theorem 1 to combine this number with the probability of A to get the final non-conditional probability that β∗ is the global solution with probability at least
≥ 1− C4 exp(−C5snγ/2 ln p)−2p exp( −(λlasso)2nbua
8σ2 )
≥ 1− C4 exp(−C5snγ/2 ln p)−2 exp( −( 2buanγ/2 − 8) ln p
8 )
≥ 1− C4 exp(−C5snγ/2 ln p)−C6 exp(−C7bunγ/2 ln(p))
(24)
For some constants C4, C5 C6 and C7.
Note that these constants, as well as C1, C2 and C3, are dependent only on the value of a Q and , as far as problem dependencies are concerned. Thus given that a Q and are chosen to be any positive constant value, as in the statement of Corollary 1, C1 through C7 are problem independent, which is the desired result.
Proof of Corollary 2. The first result is simply Corollary 2 in (Fan et al., 2014b). If we initialize the LLA algorithm with βlasso, the solution to LASSO using λlasso as the LASSO constant, then the LLA algorithm converges to the oracle solution in 2 iterations with probability 1 − φ0 − φ1 − φ2. However we still need to solve for the actual values of φ0, φ1, φ2 for GLM.
First consider φ0 = P ( ∥∥βlasso − βtrue∥∥
max > a0λ). Similar to Lemma B.1. in Bickel et al.
(2009), to bound this we will start by noticing that for the lasso penalized loss functionQlasso(β) =∑ i∈N l(β, xi, yi) + λ lasso ∑ j∈P |βj | we have that Qlasso(βlasso) ≤ Qlasso(βtrue). If we then
let δ` = βlasso − βtrue we can use the same tactic as in the derivation of 8 to get bl2n ∥∥Xδ`∥∥2 − 1 nW ᵀXδ` ≤ λlasso ∑ j∈P |βtruej | − |βlassoj |, which can then be rearranged to get.
bl 2n ∥∥Xδ`∥∥2 − 1 n ∑ j∈P |W ᵀXj ||δ′j | ≤ λlasso ∑ j∈P |βtruej | − |βlassoj |. (25)
Next letA = ⋂ j∈P{| 1 nW ᵀXj | ≤ λlasso/2}. We can combine this with 25 to get that bl2n ∥∥Xδ`∥∥2 +
λlasso/2 ∑ j∈P |βlassoj − βtruej | ≤ λlasso ∑ j∈P |βlassoj − βtruej |+ λlasso ∑ j∈P |βtruej | − |βlassoj | conditional on A. From this notice that the right term goes to zero when βtruej = 0 so we then have that bl2n
∥∥Xδ`∥∥2 + λlasso/2∑j∈P |βlassoj − βtruej | ≤ λlasso∑j∈S |βlassoj − βtruej | + |βtruej | − |βlassoj |. Using the triangle inequality and the definition of δ` we can simplify this to
bl 2n ∥∥Xδ`∥∥2 + λlasso 2 |δ`| ≤ 2λlasso|δ`S | (26)
conditional on A. By relaxing different parts of the equation, this can be further simplified to both bl n ∥∥Xδ`∥∥2 ≤ 3λlasso|δ`S | ≤ 3λlassos1/2||δ`S ||2 and |δ`Sc | ≤ 3|δ`S |. Note that the second of these shows that δ` satisfies the constraint for the RE condition 1. Therefor we have that ‖Xδ
`‖2 n‖δ`‖2 ≥ re.
If this is combined with the first of the two equations, we can get that 1 n1/2 ∥∥Xδ`∥∥ ≤ 3λlassos1/2 bl(re)1/2 conditional on A. Next, using this we can show that conditional on A we have that ∥∥δ`∥∥ max
≤ ∥∥δ`∥∥ ≤∥∥Xδ`∥∥2 2 /( ∥∥δ`∥∥nre) ≤ 3λlassos1/2blre < a0λ if λ > 3λlassos1/2bla0re . This is the inverse of the con-
dition that defines φ0. Thus, we can bound φ0 with φ0 ≤ P (Ac) = P ( ⋃ j∈P | 1 nW
ᵀXj | > λlasso/2) = P ( ⋃ j∈P |W ᵀXj |/ ‖Xj‖ > nλlasso/(2 ‖Xj‖)) ≤ pP (|〈W, v〉| > λlasson 2‖Xj‖ ) ≤ pP (|〈W, v〉| > λ lasso(nbua) 1/2
2 ) ≤ 2p exp −(λlasso)2nbua 8σ2 which uses both (A1)(ii) and (A2) as long
as λ > 3λ lassos1/2
bla0(re) per (A2). Next consider φ1 = P ( ∥∥∥OScp̃`n(βoracle∥∥∥max ≥ a1λ)
φ1 =P ( ∥∥∥OScp̃`n(βoracle∥∥∥max ≥ a1λ) (27)
=P (∃j ∈ P : |Oj`n(βoracle)| ≥ a1λ) (28)
=P (∃j ∈ P : | 1 n ∑ i∈N [ψ′(xᵀi β oracle)xi,j − yixi,j ]| ≥ a1λ) (29)
=P (∃j ∈ P : | 1 n ∑ i∈N [ψ′(xᵀi β oracle)xi,j − ψ′(xᵀi β true)xi,j +Wixi,j ]| ≥ a1λ) (30)
≤P ( 1 n |Xᵀj (ψ ′(Xβoracle)− ψ′(Xβtrue) +W )| ≥ a1λ) (31) ≤P ( 1 n |Xᵀj (ψ ′(Xβoracle)− ψ′(Xβtrue))|+ |W ᵀXj | ≥ a1λ) (32) ≤P ( 1 n ‖Xj‖
∥∥ψ′(Xβoracle)− ψ′(Xβtrue)∥∥+ |W ᵀXj | ≥ a1λ) (33) ≤P ( 1
n ∥∥ψ′(Xβoracle)− ψ′(Xβtrue)∥∥+ |W ᵀXj |/ ‖Xj‖ ≥ a1λ ‖Xj‖−1) (34) ≤P (
∥∥ψ′(Xβoracle)− ψ′(Xβtrue)∥∥+ |W ᵀXj |/ ‖Xj‖ ≥ (abun)1/2a1λ) (35) ≤P (bu
∥∥Xβoracle −Xβtrue∥∥+ |W ᵀv| ≥ (abun)1/2a1λ) (36) ≤P (bu ‖Xδo‖+ |W ᵀv| ≥ (abun)1/2a1λ) (37)
where v ∈ <n is some vector with ‖v‖ = 1 as indicated in (A2) and δo = βoracle − βtrue. From this, using Demorgan’s law and the union bound, we notice that P (A + B ≥ C) ≤ P (A ≥ C/2) + P (B ≥ C/2) which can be used to further simplify
φ1 ≤ P (bu ‖Xδo‖+ |W ᵀv| ≥ a1λ(abun)1/2) (38) ≤ P (‖Xδo‖ ≥ (1/2)a1λ(an/bu)1/2) + P (|W ᵀv| ≥ (1/2)a1λ(abun)1/2) (39)
We can then simplify both terms individually. For the first term, P (bu ‖Xδo‖ ≥ (1/2)a1λ(abun)
1/2), given the fact that the oracle solution and true solution have the same support, the oracle solution must be in the Γ = 0 level set of the true solution. Using similar arguments to Lemma 5, we have that bl2n ‖Xδ o‖2 ≤ 1nW ᵀXδo From here Lemma 2 can
be applied since we know ∥∥βoracle − βtrue∥∥
0 ≤ s. With some simplification this gives that ‖Xδo‖ ≤ 2bl (maxSp̃:|Sp̃|=s ∥∥∥ŨᵀSp̃W∥∥∥) Utilizing Lemma 3 with s in place of p̃ shows that
P [ maxSp̃:|Sp̃|=s 2 bl ∥∥∥ŨᵀSp̃W∥∥∥ ≥ 2blσ√s+ 2√st1 + 2t1] ≤ (pes )s exp(−t1). This is the first half of φ1 as long as (1/2)a1λ(an/bu)1/2 ≥ 2blσ √ s+ 2 √ st+ 2t which is equivalent to the as-
sumed condition λ ≥ 4σ √ s+2 √ st1+2t1
bla1(an/bu)1/2 Next, the second term can be easily bounded using (A2):
P (|W ᵀv| ≥ (1/2)a1λ(abun)1/2) ≤ 2 exp(−a 2 1λ 2abun 8σ2 )
Therefor φ1 ≤ (pes ) s exp(−t1) + 2 exp(−a 2 1λ 2abun 8σ2 ) Next consider φ2 = P ( ∥∥βoracleS ∥∥min ≤ aλ) First, given the assumption ∥∥βtrue]∥∥min > (a +
1)λ we can see that φ2 = P ( ∥∥βoracleS ∥∥min ≤ aλ) ≤ P (∥∥βoracleS − βtrueS ∥∥max > λ) ≤
P ( ∥∥βoracle − βtrue∥∥
2 > λ) = P (‖δo‖2 > λ). Next, since we know that the support of βoracle
and βtrue is S, we know that |δoSc | = 0 ≤ 3|δ o S | which is the constraint for the RE condition. Therefor we know that ‖Xδ o‖2
n‖δo‖2 ≥ re. With this and a similar line of argument as in φ1 we get
that φ2 ≤ P (‖δo‖ > λ) ≤ P (‖Xδo‖ > λ √ nre) ≤ P ( 2bl (maxSp̃:|Sp̃|=s ∥∥∥ŨᵀSp̃W∥∥∥ > λ√nre) =
P (maxSp̃:|Sp̃|=s ∥∥∥ŨᵀSp̃W∥∥∥ > λ bl√nre2 ≥ σ√s+ 2√st2 + 2t2) ≤ (pes )s exp(−t2) assuming that λ bl √ nre 2 ≥ σ √ s+ 2 √ st+ 2t which is equivalent to the condition λ ≥ 2σ √ s+2 √ st2+2t2 bl √ nre
This, combined with the fact that for MCP, a0 = a1 = a2 = 1 shows the first result.
The second result can be seen by first noting all assumptions of Theorem 1 part 2 are satisfied, where (A3) with r4s is implied by 7. Thus by using the same arguments as in Theorem 1 part 2 which shows that the oracle solution is unique and that the global solution is the oracle solution with some probability, since the global solution is almost surely S3ONC with Γ = 0. If we use t = t3 and t′ = t4 we get that the probability that the global solution is not the oracle solution as φ3 ≤ exp(−t3 + p̃ ln(pep̃ )) + exp(−(p̃ ∗ + 1)(t4 − ln p)) · 1−exp(−(p−p̃ ∗)(t4−ln p))
1−exp(−t4+ln p) This combined with the first result shows that the LLA algorithm converges to the global solution in 2 iterations with probability 1− φ0 − φ1 − φ2 − φ3 which is the second result.
A.2 CENTRAL LEMMAS AND THEIR PROOFS
Lemma 1. Let β∗ be a S3ONC solution to 2. If assumption (A1) holds, then P [|β∗j | /∈ (0, aλ),∀j ∈ {1, 2, ..., p}] = 1.
Proof of Lemma 1. First, define events γj and δj as
γj := { ∂2Q(β) (∂βj)2 ∣∣∣∣ β=β∗ ≥ 0 } (40)
δj := { |β∗j | ∈ (0, aλ) } . (41)
First, for any given j ∈ P , we solve for P [γj ∩ δj ] given our assumptions. We can start with ∂2Q(β) (∂βj)2 ∣∣∣ β=β∗ ≥ 0 which gives us 1/n n∑ i=1 ψ′′ (xᵀi β ∗)x2i,j + P ′′ λ (|β∗j |) ≥ 0. We can rearrange this
to get bu n∑ i=1 x2i,j ≥ n∑ i=1 ψ′′(xᵀi β ∗)x2i,j ≥ −nP ′′λ (|β∗|) = n/a where we get the leftmost inequality from assumption (A1) part (i) and the rightmost equality from the definition of MCP. More concisely we have that bu ‖Xj‖2 ≥ n/a which contradicts (A1) part (ii). Therefor we know P [γj ∩ δj ] = 0. It should also be noted that P [γcj ] = 0 since β
∗ satisfies S3ONC conditions. Thus, by applying Demorgan’s law and then the union bound, it can be obtained that
0 = P [γj ∩ δj ] = 1− P [γcj ∪ δcj ] ≥ 1− P [γcj ]− P [δcj ] = 1− P [δcj ] = P [δj ]. (42)
We can then apply this result to all indices to get that P [δj ] = 0 for all j ∈ {1, 2, ..., p}, which is the desired result.
Lemma 2. Consider an arbitrary S3ONC solution β∗ to 2 with MCP. Given the event that for some integer p̃ : ‖β∗ − βtrue‖0 ≤ p̃, then |W ᵀXδ∗| ≤ ( maxSp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥) ‖Xδ∗‖ , a.s. Where
(ŨSp̃)i,j := { USp̃ , if j ∈ Sp̃ 0, else
and USp̃ ∈ <n×p̃ is defined as in the following Thin SVD: XSp̃ = USp̃DSp̃VSp̃ .
Proof. Denote δ∗ := (δ∗j ) = β ∗ − βtrue, Sp̃ := (j : δ∗j 6= 0) ⊆ P , δ∗Sp̃ := (δ ∗ j : j ∈ Sp̃) and XSp̃ := (xij : i ∈ N , j ∈ Sp̃). By assumption, we know that ‖δ∗‖0 ≤ |Sp̃| = p̃.
First decompose XSp̃ using Thin SVD to get XSp̃ = USp̃DSp̃VSp̃ where USp̃ ∈ <n×p̃. Note that since and UᵀSp̃USp̃ = I we have that for any υ ∈ < p̃ we have ∥∥DSp̃VSp̃υ∥∥2 =
(DSp̃VSp̃υ) ᵀI(DSp̃VSp̃υ) = υ ᵀV ᵀSp̃D ᵀ Sp̃ UᵀSp̃USp̃DSp̃VSp̃υ = υ ᵀXᵀSp̃XSp̃υ = ∥∥XSp̃υ∥∥2. therefor we can obtain that
|W ᵀXδ∗| = |W ᵀXSp̃δ∗Sp̃ | ≤ ∥∥W ᵀUSp̃∥∥ ∥∥∥DSp̃VSp̃δ∗Sp̃∥∥∥
= ∥∥∥UᵀSp̃W∥∥∥ ∥∥∥XSp̃δ∗Sp̃∥∥∥ ≤ ( maxSp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥) ‖Xδ∗‖ , a.s. (43) Where
(ŨSp̃)i,j := { USp̃ , if j ∈ Sp̃ 0, else .
Lemma 3. Consider an arbitrary S3ONC solution β∗ to 2 with MCP. If (A2) holds, then for some integer p̃ ≤ p, P [ maxSp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥ ≤ σ√p̃+ 2√p̃t+ 2t] ≥ 1− (pep̃ )p̃ exp(−t). Where (ŨSp̃)i,j := { USp̃ , if j ∈ Sp̃ 0, else
and USp̃ ∈ <n×p̃ is defined as in the following Thin SVD: XSp̃ = USp̃DSp̃VSp̃ . Proof. We attempt to bound (
maxSp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥). Given that we now have W multiplied by a square matrix, we can apply Lemma 9. In the Lemma, let Σu = ŨSp̃Ũ ᵀ Sp̃
. The fact that ΣuΣu = Σu means that Σu is an idempotent matrix with ‖Σu‖ ≤ 1 and Tr(Σu) = rank(Σu) ≤ rank(ŨSp̃) ≤ rank(USp̃) ≤ p̃. Lemma 9 then states that that P [∥∥∥ŨᵀSp̃W∥∥∥ ≤ σ√p̃+ 2√p̃t+ 2t] ≥ 1−exp(−t). From this we can show that
P [ max
Sp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥ ≤ σ√p̃+ 2√p̃t+ 2t] ≥ 1− (pp̃ ) exp(−t) ≥ 1− (pe p̃ )p̃ exp(−t). (44)
Where the first inequality can seen by noting that if ηk ∈ <k is a sequence of i.i.d random variables and θ ∈ < is a scalar, by applying De Morgan’s Law and then using the union bound, it can be obtained that P [maxk∈K ηk ≤ θ] = P [ ⋂ k∈K ηk ≤ θ] = 1−P [ ⋃ k∈K ηs ≥ θ] ≥ 1− ∑ k∈K P [ηk ≥ θ] = 1− |K|(1− P [ηk ≤ θ]) which yields the same inequality as in 44. This is the desired result.
Lemma 4. Consider an arbitrary S3ONC solution β∗ to 2 with MCP. Let Assumptions (A1) and (A2) hold. Given the simultaneous occurrence of (i) the event that Q(β∗) ≤ Q(βtrue) + Γ holds for some Γ ≥ 0; (ii) the event that for some integer p̃ : ‖β∗ − βtrue‖0 ≤ p̃. Then for any t > 0, 1 n ‖X(β ∗ − βtrue)‖2 ≤ 4σ 2
b2l n (p̃+2 √ p̃t+2t)+ 8bl min{ ∑ j∈S P ′ λ(|β∗j |)|β∗j |, Pλ(aλ)(|S|−‖β∗‖0)+
Γ} holds with probability at least 1− exp(−t+ p̃ ln(pep̃ )).
If in addition (A3) holds with p̃∗ ≥ p̃, then 1n ‖X(β ∗ − βtrue)‖2 ≤ 8σ 2
b2l n
( p̃+ 2 √ p̃t+ 2t ) +
8 bl min{λ2(|S| − ‖β∗‖0)r −1 p̃ , Pλ(aλ)(̇|S| − ‖β∗‖0) + Γ} holds where rp̃ > 0 for any t > 0 with probability at least 1− exp(−t+ p̃ ln(pep̃ )).
Proof. First, denote δ∗ := (δ∗j ) = β ∗ − βtrue, Sp̃ := (j : δ∗j 6= 0) ⊆ P , δ∗Sp̃ := (δ ∗ j : j ∈ Sp̃) and XSp̃ := (xij : i ∈ N , j ∈ Sp̃). By assumption, we know that ‖δ∗‖0 ≤ |Sp̃| = p̃. Further, let us denote
T1 := min ∑ j∈S P ′λ(|β∗j |)|βtruej |, ∑ j∈S P ′λ(|β∗j |)|β∗j − βtruej |, Pλ(aλ)(|S| − ‖β∗‖0) + Γ  . (45) We now start to define the desired bound by applying the second part of Lemma 8. The result simplified using the above definitions becomes
bl 2n ‖Xδ∗‖2 ≤ 1 n W ᵀXδ∗ + T1, a.s. (46)
Next, since all assumptions for Lemma 1 are satisfied, we can apply it to get
bl 2n ‖Xδ∗‖2 ≤ 1 n
( max
Sp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥) ‖Xδ∗‖+ T1. (47) We can then complete the square, solving for 1√
n ‖Xδ∗‖ to get
1√ n ‖Xδ∗‖ ≤ 1 bl √ n max Sp̃:|Sp̃|=p̃
∥∥∥ŨᵀSp̃W∥∥∥+ √( 1
bl √ n max Sp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥)2 + 2bl T1 (48) ≤2 √( 1
bl √ n max Sp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥)2 + 2bl T1, (49) where the last inequality holds due to the value inside the square root being larger than the term outside. From here, squaring both sides gives us
1 n ‖Xδ∗‖2 ≤ 4
b2l n
( max
Sp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥)2 + 8bl T1. (50) Finally, by applying the second part of Lemma 1 we get
1 n ‖Xδ∗‖2 ≤ 4σ
2
b2l n
( p̃+ 2 √ p̃t+ 2t ) + 8
bl T1. (51)
With probability at least 1 − (pep̃ ) p̃ exp(−t). Thus by the definition of T1, the first result of the lemma has been shown.
For the second part we look to bound the central term of T1. We first notice (a) that since assumption (A1) holds, Corollary 4 indicates that if β∗j 6= 0 ⇒ |β∗j | ≥ aλ for all j ∈ P; (b) that for this range of β∗j , P ′ λ(|β∗j |) = 0; (c) that per the definition of MCP 0 ≤ P ′λ(|β∗j |) ≤ λ for any β∗j ∈ <. If we
combine these observations with 45 and the definition of δ∗, we can see that T1 ≤ ∑ j∈S P ′λ(|β∗j |)|δ∗|
≤ λ √ |S| − ‖β∗S‖0 · ‖δ
∗‖. From this, given that assumption that, for this second result (A3) holds with p̃∗ ≥ p̃, and rp̃ ≥ rp̃∗ ≥ 0 we can use (A3) part (iii) to show that T1 ≤ λ √ |S| − ‖β∗S‖0 · ‖Xδ∗‖√ nrp̃
. Since this holds almost surely, it can then be combined with 47 to get
bl 2n ‖Xδ∗‖2 ≤ 1 n
( max
Sp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥) ‖Xδ∗‖+ λ√|S| − ‖x∗S‖0 · ‖Xδ∗‖√nrp̃ . (52) We can then multiply by 2 √ n/bl ‖Xδ∗‖ to get
1√ n ‖Xδ∗‖ ≤ 2 bl √ n max Sp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥+ 2λbl√rp̃ √ |S| − ‖x∗S‖0. (53)
We then square both sides and use the rule that (A+B)2 ≤ 2A2 + 2B2 to get
1 n ‖Xδ∗‖2 ≤
[ 2
bl √ n max Sp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥+ 2λbl√rp̃ √ |S| − ‖x∗S‖0 ]2 (54)
≤ 8 b2l n max Sp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥2 + 8λ2b2l rp̃ (|S| − ‖x∗S‖0). (55) Combining this with 50 yields that
1 n ‖Xδ∗‖2 ≤ 8
b2l n max Sp̃:|Sp̃|=p̃ ∥∥∥ŨᵀSp̃W∥∥∥2 + 8bl min { λ2 rp̃ (|S| − ‖x∗S‖0), T1 } . (56)
Finally, by applying the second part of Lemma 1 and noting from (45) that T1 ≤ Pλ(aλ)(|S| − ‖β∗‖0) + Γ we see that
1 n ‖Xδ∗‖2 ≤ 8 b2l n σ2(p̃+ 2
√ p̃t+ 2t) + 8
bl min
{ λ2
rp̃ (|S| − ‖x∗S‖0), Pλ(aλ)(|S| − ‖β
∗‖0) + Γ } .
(57)
With probability at least 1− (pep̃ ) p̃ exp(−t) which is the desired result.
Lemma 5. Let Assumptions (A1) and (A2) hold. Consider a solution β∗ satisfying S3ONC of 2. Assume thatQ(β∗) ≤ Q(βtrue) + Γ holds for an arbitrary Γ > 0. For any integer p̃ : 2|S| ≤ p̃ ≤ p if the penalty parameters (a, λ) satisfy that Pλ(aλ) > σ 2
2nbl (1 + 2
√ t+ 2t) + σ2 n |S|(1+2 √ t+2t)+Γbl
bl(p̃+1−2|S|) , for an arbitrary t > 0, then ‖β∗ − βtrue‖0 ≤ p̃ with probability at least 1 − exp(−(p̃ + 1)(t − ln p)) · 1−exp(−(p−p̃)(t−ln p))1−exp(−t+ln p) .
Proof. We start from the useful inequality defined in 8
bl 2n ‖Xδ∗‖2 − 1 n W ᵀXδ∗ + ∑ j∈S Pλ(|β∗j |) ≤ ∑ j∈S Pλ(|βtruej |) + Γ, (58)
where δ∗ = β∗ − βtrue. Next, conditioning on the fact (i) that β∗ is S3ONC, (ii) that all assumptions for Corollary 4 are satisfied (which implies that Pλ(|β∗j |) ∈ {0, Pλ(aλ)}) and (iii) that Pλ(|βtruej |) ≤ Pλ(aλ) we have that
bl 2n ‖Xδ∗‖2 − 1 n W ᵀXδ∗ + ‖β∗‖0 · Pλ(aλ) ≤ |S| · Pλ(aλ) + Γ (59)
Now consider an event E1 := {‖β∗ − βtrue‖0 = p̃+ k} for an arbitrary integer k : 1 ≤ k ≤ p− p̃ Conditioning on this event, we may denote and Sp̃+k ⊆ P such that δ∗j 6= 0 for all j ∈ Sp̃+k. By assumption we can ensure that |Sp̃+k| = p̃+ k. Also denote by XSp̃+k = (xij : i ∈ N , j ∈ Sp̃+k) and let δ∗Sp̃+k := (δ ∗ j : j ∈ Sp̃+k). Note that conditional on E1, the first part of Lemma 1 (using p̃ + k in place of p̃ in Lemma 1) can be used to bound W ᵀXδ∗ in 59. Additionally, by definition ‖βtrue‖0 = |S| and conditional on E1 we can apply the substitution ‖β∗‖0 ≥ p̃ + k − |S|. This gives us
bl 2 ∥∥∥∥Xδ∗√n ∥∥∥∥2− 1√n ( max Sp̃+k:|Sp̃+k|=p̃+k ∥∥∥ŨᵀSp̃+kW∥∥∥)∥∥∥∥Xδ∗√n ∥∥∥∥ ≤ −(p̃+k− 2|S|) ·Pλ(aλ) + Γ (60)
In order for this equation to be feasible, we know that the quadratic formula must have real roots. Therefor
( max
Sp̃+k:|Sp̃+k|=p̃+k
∥∥∥∥∥ Ũ ᵀ Sp̃+k W √ n ∥∥∥∥∥ )2 − 4[bl/2][(p̃+ k − 2|S|) · Pλ(aλ)− Γ] ≥ 0 (61)
Now consier another event E2(t) := {max|Sp̃+k|=p̃+k ||U ᵀ Sp̃+k
W || ≤ σ √ p̃+ k · √ 1 + 2 √ t+ 2t}
for an arbitrary t > 0. Conditioning on the occurence of E1 ∩ E2(t) we can show, using first
E2(t) and then 61, that σ 2(p̃+k) n · (1 + 2
√ t + 2t) ≥ ( max|Sp̃+k|=p̃+k ∥∥∥ŨᵀSp̃+kW∥∥∥√ n )2 ≥ 2bl[(p̃ +
k − 2|S|) · Pλ(aλ) − Γ] almost surely, which contradicts with the assumption on the parameters (a, λ). This can be seen starting from our original assumption that Pλ(aλ) > σ 2
2nbl (1 + 2
√ t +
2t)+ σ2
n |S|(1+2 √ t+2t)+Γbl
bl(p̃+1−2|S|) ≥ σ2 2nbl (1+2
√ t+2t)+ σ2 n |S|(1+2 √ t+2t)+Γbl
bl(p̃+k−2|S|) We can then multiply both
(outer) sides by 2bl(p̃+k−2|S|) and rearrange to get σ 2 n (p̃+k) · (1 + 2 √ t+ 2t) < 2bl[(p̃−2|S|+ k) ·Pλ(aλ)−Γ]. Given this contradiction, we know that P [E1 ∩ E2(t)] = 0. Therefore, again using the union bound combined with DeMorgan’s law we get that P [E1∩E2(t)] ≥ 1−P [Ec1 ]−P [E2(t)c] which, can be simplified to
P [E2(t)c] ≥ P [E1] (62)
Since all assumptions of the Lemma 3 are satisfied, we can next use it to bound P [E2(t)c]. By taking the compliment of the result in the second half of Lemma 1, we get, for some t′, that P [ maxSp̃+k:|Sp̃+k|=(p̃+k) ∥∥∥ŨᵀSp̃+kW∥∥∥ ≥ σ√(p̃+ k) + 2√(p̃+ k)t′ + 2t′] ≤ ( pep̃+k ) p̃+k exp(−t′) ≤ pp̃+k exp(−t′). If we then let t′ = (p̃ + k)t we get
P [ maxSp̃+k:|Sp̃+k|=(p̃+k) ∥∥∥ŨᵀSp̃+kW∥∥∥ ≥ σ√p̃+ k ·√1 + 2√t+ 2t] ≤ pp̃+k exp(−(p̃ + k)t). Thus we have that pp̃+k exp(−(p̃+ k)t) ≥ P [E2(t)c] which can be combined with 62 to show
pp̃+k exp(−(p̃+ k)t ≥ P [ ∥∥β∗ − βtrue∥∥
0 = p̃+ k] ∀k ∈ Z : 1 ≤ k ≤ p− p̃. (63)
With this, we can solve for our desired value
P [∥∥β∗ − βtrue∥∥ 0 ≤ p̃ ] = 1− P [∥∥β∗ − βtrue∥∥ 0 ≥ p̃+ 1 ] = 1− p−p̃∑ k=1 P [∥∥β∗ − βtrue∥∥ 0 = p̃+ k ] ≥ 1−
p−p̃∑ k=1 exp((p̃+ k)(ln p− t))
= 1− exp(−(p̃+ 1)(t− ln p)) · 1− exp(−(p− p̃)(t− ln p)) 1− exp(−t+ ln p)
(64)
Which is the desired result.
Lemma 6. Consider an arbitrary S3ONC solution β∗ to 2 with MCP. Let Assumptions (A1) and (A3) with p̃∗ ≥ p̃ hold. Assume the satisfaction of ‖β∗ − βtrue‖0 ≤ p̃ and Event Ea(p̃) := { 1n ‖X(β ∗ − βtrue)‖2 ≤ 8σ 2
b2l n
( p̃+ 2 √ p̃t+ 2t ) + 8bl min{ λ2 rp̃ (|S|−‖β∗‖0), Pλ(aλ)·(|S|−‖β∗‖0)+
Γ}. If the sub-optimality gap satisfies Γ < Pλ(aλ) − σ 2
bln
( p̃+ 2 √ p̃t+ 2t ) . If the minimum signal
strength satisfies ‖βtrueS ‖min > √ 8σ2
rp̃b2l n
( p̃+ 2 √ p̃t+ 2t ) + 8rp̃bl min{ λ2 rp̃ |S|, Pλ(aλ)|S|+ Γ} then
β∗ is the oracle solution to 2.
If in addition we have the satisfaction of ‖βopt − βtrue‖ ≤ p̃ and the event Eb(p̃) := { 1n ‖X(β opt − βtrue)‖2 ≤ 8σ 2
b2l n
( p̃+ 2 √ p̃t+ 2t ) + 8bl min{ λ2 rp̃ (|S| − ‖β∗‖0), Pλ(aλ) · (|S| −
‖β∗‖0) + Γ} then β∗ is both the oracle solution and the global solution to 2.
Proof. First, let us denote β∗ − βtrue = δ∗. We start by combining Eα(p̃) and (A3)iii, which is possible due to our assumption ‖β∗ − βtrue‖0 ≤ p̃. This gives us
8σ2
b2l n
( p̃+ 2 √ p̃t+ 2t ) + 8
bl min{λ
2
rp̃ (|S| − ‖β∗‖0), Pλ(aλ) · (|S| − ‖β ∗‖0) + Γ}
≥ 1 n ‖Xδ∗‖2 ≥ rp̃ ‖δ∗‖2 a.s.
(65)
From here if we relax |S| − ‖β∗S‖0 to just |S|, the definition of δ ∗ and note that ‖δ∗‖ ≥ ∥∥δ∗j∥∥, we can obtain the following
√ 8σ2
rp̃b2l n
( p̃+ 2 √ p̃t+ 2t ) + 8
rp̃bl min{λ
2
rp̃ |S|, Pλ(aλ)|S|+ Γ}
≥ ∥∥β∗j − βtruej ∥∥ ≥ |βtruej | − |β∗j |,
(66)
almost surely. From this we can bound |β∗j | using the square root term and |βtruej |, so we know that if |βtruej |− √ 8σ2
rp̃b2l n
( p̃+ 2 √ p̃t+ 2t ) + 8rp̃bl min{ λ2 rp̃ |S|, Pλ(aλ)|S|+ Γ} > 0 then |β∗j | > 0. From
this we can obtain the inequality
‖β∗S‖0 ≥ ∑ j∈S I ( |βtruej | − √ 8σ2 rp̃b2l n ( p̃+ 2 √ p̃t+ 2t ) + 8 rp̃bl min{λ 2 rp̃ |S|, Pλ(aλ)|S|+ Γ} > 0 ) (67)
almost surely. We can then combine this with our minimum signal strength assumption to get that
‖β∗S‖0 = |S| a.s. (68)
We can combine this with equation 65, by focusing on the second part of the minimum term and noting the right side is always positive to get
8σ2
b2l n
( p̃+ 2 √ p̃t+ 2t ) + 8
bl (−Pλ(aλ) ‖β∗Sc‖0 + Γ) ≥ 0 a.s. (69)
which can be simplified into
σ2
bln
( p̃+ 2 √ p̃t+ 2t ) + Γ ≥ Pλ(aλ) ‖β∗Sc‖0 a.s. (70)
thus, it can be seen that if Pλ(aλ) > σ 2
bln
( p̃+ 2 √ p̃t+ 2t ) + Γ then 1 > ‖β∗Sc‖0 = 0. This is
satisfied by the assumption that Pλ(aλ)− σ 2
bln
( p̃+ 2 √ p̃t+ 2t ) > Γ.
Finally, because β∗ is an S3ONC solution, it has to satisfy FONC. Per 2, this means that β∗ ∈ arg inf{ 1n ∑ i∈N `(β, xi, yi) + ∑ j∈P P ′ λ(|β∗j |)|βj | : β ∈ <p}. Due to Corollary 4, we know that the penalty term goes to 0 since either β∗j = 0 or P ′(|β∗j |) = P ′(|aλ|) = 0. Further we know that β∗j = 0 for all j ∈ Sc. Threfore we know that
β∗ ∈ arg inf{ 1 n ∑ i∈N `(β, xi, yi) : β ∈ <p, βj = 0,∀j ∈ Sc} a.s. (71)
Given that the expression on the right is the definition of the oracle solution, we have shown the first result.
Next, Consider βopt which is the global optimal solution to 2. Given that the S3ONC conditions are necessary, βopt must be an S3ONC solution. With this fact and the assumption of Eb(p̃), we have the same set of assumptions for βopt as we had for β∗. Thus the same sequence of arguments can be used to show that
βopt ∈ arg inf{ 1 n ∑ i∈N `(β, xi, yi) : β ∈ <p, βj = 0,∀j ∈ Sc} a.s. (72)
Finally, per the strict convexity of our loss function as implied by (A1) we can see that the infimum of the above problem is unique. Therefor
β∗ = arg inf{ 1 n ∑ i∈N `(β, xi, yi) : β ∈ <p, βj = 0,∀j ∈ Sc} = βopt a.s. (73)
Which is the second result.
A.3 ADDITIONAL LEMMAS
Lemma 7. The RE condition in 1 implies (A3) with r4s ≥ re > 0 and p̃∗ ≥ 4s.
Proof. As in Lemma 1 in (Liu et al., 2017).
Lemma 8. Let β∗ be a S3ONC solution to 2 Given (A1) and that Q(β∗) ≤ Q(βtrue) + Γ holds for some Γ ≥ 0 then
bl 2n ‖Xδ∗‖2 − 1 n W ᵀXδ∗
≤ min ∑ j∈S P ′λ(|β∗j |)|βtruej |, ∑ j∈S P ′λ(|β∗j |)|β∗j − βtruej |, Pλ(aλ)(|S| − ‖β∗‖0) + Γ  , a.s. (74)
Proof. First, we know that β∗ ∈ arg min β { n∑ i=1 `(β, xi, yi) + p∑ j=1 P ′λ(|β∗|)|βj |} because the KKT
conditions are the same as FONC which β∗ satisfies. This gives us that n∑ i=1 `(β∗, xi, yi) +
p∑ j=1 P ′λ(|β∗|)|β∗j | ≤ n∑ i=1 `(βtrue, xi, yi) + p∑ j=1 P ′λ(|β∗|)|βtruej |. This can be used along the same lines as the level set inequality in the derivation for 8 to get bl2n ‖Xδ ∗‖2 − 1nW
ᵀXδ∗ ≤ p∑ j=1 P ′λ(β ∗ j )(|βtruej | − |β∗j |)
The first two terms of the min function are easily obtained from this. The last term can be obtained from 8 by noting that due to Corollary 4, β∗ /∈ (0, aλ) and that Pλ(aλ) = Pλ(β) ∀β ≥ aλ. This gives us that bl2n ‖Xδ ∗‖2− 1nW ᵀXδ∗ ≤ Pλ(aλ)(S−‖β∗‖0)+Γ Which is the final term to complete the desired result.
Lemma 9. Consider a subgaussian ñ-dimensional random vector W̃ ∈ <ñ as defined in (A2). Then for any V ∈ <ñ×ñ and Σv = V ᵀV then P [ ∥∥∥V W̃∥∥∥2 ≤ σ2 · (Tr(Σv) + 2√Tr(Σ2v)t+ 2 ‖Σv‖ t)] ≥ 1− exp(−t) for any t > 0 where Tr(·) denotes the trace of a matrix.
Proof. We apply Theorem 2.1 in (Hsu et al., 2012) where our W̃ ,V and Σv are equivalent to their x, A and Σ. Note their expectation condition is equivalent to our (A2) with µ = E[W ] = 0. This gives us that for all t > 0
P [ ‖VW‖2 > σ2 · (tr(Σv) + 2 √ tr(Σv)t+ 2 ‖Σv‖ t)
+ tr(Σvµµ ᵀ) · (1 + 2( ‖Σv‖
2
tr(Σ2v) t)1/2)
] ≤ exp(−t).
(75)
Given that µ = 0, the term involving tr(Σvµµᵀ) goes to zero and therefor the statement in 9 can be obtained by taking the complement of the probability bound.

Review summaries:
The paper studies the problem of global optimization of high-dimensional sparse estimators regularized by PCP (folded concave penalty). The main result is showing that under certain conditions, with high probability, the desired global solution is an oracle stationary point satisfying the so-called S^3ONC conditions. In light of this result and an existing polynomial-time algorithm for finding an S^3ONC solution, the global solution can be recovered with polynomial computational complexity. Numerical evidence is provided to show the theoretical predictions. 

The authors study the high-dimensional sparse estimation problems, which is one of the fundamental topics in both machine learning and optimization communities. In the literature, the folded concave penalty (FCP) methods have been shown to enjoy the strong oracle property for high-dimensional sparse estimation.  While LASSO solutions can be easily computed, such solutions do not admit unbiasedness and oracle property or require strong conditions. Therefore,  the problems with FCP have received much attentions and been studied in terms of hardness and approximability of the problems recently .