Title: FULLY POLYNOMIAL-TIME RANDOMIZED APPROXI-MATION SCHEMES FOR GLOBAL OPTIMIZATION OF HIGH-DIMENSIONAL MINIMAX CONCAVE PENALIZED GENERALIZED LINEAR MODELS

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 (S 3 ONC). Given that the desired S 3 ONC 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 S 3 ONC 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.

1. INTRODUCTION:

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 = (x 1 , ..., x n ) be the n×p design matrix with x i = (x i1 , ..., x ip ) , i = 1, ..., n, and Y = (y 1 , ..., y n ) 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[y i ] = ψ (x i β true ) for some known link function ψ : Θ → , where Θ ⊆ and β true = (β true 1 , ..., β true p ) 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:

EQUATION

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:

EQUATION

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 (S 3 ONC) provable admit FPRAS's. The S 3 ONC are weaker conditions than the standard second-order KKT conditions. As per this result, all S 3 ONC-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 S 3 ONC 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 S 3 ONC. 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, v Q is intended as (v j : j ∈ Q). For any set Q, we denote the complement as Q c . In particular, let S be the true support set, that is, S := {j : β true j = 0} and its complement is S c . 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.

2. SETUPS, PRELIMINARIES, AND ASSUMPTIONS:



2.1. SETUPS AND ASSUMPTIONS:

Our analysis will focus on sparse GLMs that have a fixed design matrix and satisfy the following assumptions:

(A1) Assume that (i) b u ≥ ψ (x i β) ≥ b l > 0 for all x i β ∈ Θ;

(ii) There exists K > 0 such that the design matrix satisfies 1 n X j 2 < K for all j ∈ [p].

Let the tuning parameter a in P λ satisfy K < (b u a) -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(-t 2 /2σ

2 ), for all v : v = 1 and any t > 0; (A3) There exists a sequence {r d ≥ 0 : d = 1, 2, ..., p} such that the following are satisfied:

(i) For any d 1 , d 2 : 1 ≤ d 1 ≤ d 2 ≤ p, we have r d1 ≥ r d2 ;

(ii) There exists some p * : 2|S| ≤ p * ≤ p such that r p * > 0;

(iii) For all d : 1

≤ d ≤ p, it holds that n -1 Xβ 2 ≥ r d β 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. 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 r e > 0, it holds that

1 n Xδ 2 ≥ r e δ 2 for all δ ∈ | Ŝ|=s C( Ŝ)

where

C( Ŝ) := {δ := (δ i ) ∈ p : |δ Ŝc | ≤ 3|δ Ŝ |}, δ Ŝc := (δ j : j ∈ Ŝc ), and δ Ŝ := (δ j : j ∈ Ŝ).

Furthermore, the largest possible r e 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 S 3 ONC:

Our results focus on the S 3 ONC 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 S 3 ONC 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 β * ) -y i ] x i + P λ (|β * j |)∂(|β * j |), 1 ≤ j ≤ p (3)

where

∂(| • |) denotes the subdifferential of | • |.

Definition 3 (S 3 ONC). A solution β * satisfies the significant subspace second-order necessary condition (S 3 ONC) if it satisfies FONC and for all j ∈ {j :

β * j = 0}, ∂ 2 Q(β) (∂β j ) 2 β=β * ≥ 0 (4)

if the second derivative exists.

Remark 5. The S 3 ONC can be intuited as the second order necessary condition applied only to the dimensions where β j = 0, i.e., the significant dimensions. Since the S 3 ONC is weaker than the standard second-order KKT conditions, any algorithm that guarantees the second-order KKT conditions can be used to obtain an S 3 ONC solution, by requiring a more stringent optimality condition, may be slower than necessary. One specifically S 3 ONC guaranteeing approach, presented in (Liu & Ye, 2019) , utilizes an interior point trust region algorithm in order to guarantee an S 3 ON C solution in polynomial time. This is the scheme which will be used later in Section 4.

3. MAIN RESULTS:

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 min L(β) + λ|β|.

(5)

Theorem 1. Suppose assumptions (A1), (A2), and (A3) with any p * : p * ≥ 2|S|. Let β * be an arbitrary S 3 ONC 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 b l n p * + 2 √ p * t + 2t ; (ii) choose P λ (aλ) > σ 2 2nb l (1 + 2 √ t + 2t ) + σ 2 n |S|(1+2 √ t +2t )+Γb l b l ( p * -2|S|+1)

and (iii) assume that the minimal signal strength satisfy

β true S min > 8σ 2 r pb 2 l n p * + 2 p * t + 2t + 8 r pb l min λ 2 r p |S|, P λ (aλ)|S| + Γ

then the following two statements hold (a) β * is an oracle solution with probability at least 1 -exp(-t + p * ln

( pe p * )) -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( pe p * )) -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 S 3 ONC 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 S 3 ONC 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 β true S 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, b l ≤ 1, and s ≥ 1. Let β * be an S 3 ONC 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 n 1-γ/2 where γ ∈ [0, 1] is an arbitrary scalar. Let λ = σ re ln p n γ/2

and a ∈ [.8, 1). There exist problem independent constants

C 1 > 0,C 2 > 0 and C 2 > 0 such that if n > max    C 1 b l , C 2 s b l 2 1-γ , C 3 sσ 2 ln p β true S 2 min b 2 l r 4 e 2/γ    (6)

Then β * is the global solution to 2 with probability at least 1 -C 4 exp(-C 5 sn γ/2 ln p) -C 6 exp(-C 7 b u n γ/2 ln(p)) for problem independent constants C 4 , C 5 , C 6 and C 7 Remark 9. Corollary 1 indicates that for γ > 0, the global optimal solution coincides with computable S 3 ONC 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 S 3 ONC 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, â, k = 0 and β Lasso by solving equation 5Step 2. If Case 1: |β k j | ∈ (0, aλ) for some j = 1, ..., p, then choose an arbitrary ι ∈ {j : |β k j ∈ (0, aλ)} and solve

β k+1 ι ∈ arg min β [∇L(β k )] ι • β + P λ (|β|) s.t. (β -β k ι ) 2 ≤ δ 2 and let β k+1 j = β k j ∀j = ι. Go to Step 3.

Else Case 2: If -β k j / ∈ (0, aλ) for all j = 1, ..., p then for all j = 1, ..., p:

-If β k j = 0 then β k+1 j = â • |[∇L(β k )] j | -λ + • sign(-[∇L(β k )] j ).

-If 

|β k j | ≥ aλ, then β k+1 j = β k j -â • [∇L(β k )] j . Go to Step 3. Step 3. Algorithm stops if |β k j | / ∈ (0,

Q(β Lasso ) -Q(β opt ) • max{(1/(2a) -b u /2) -1 , 2b -1 u , (1/a -b u /2) -1 } • 1/δ 2

for any γ-accuracy S 3 ONC 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 S 3 ONC (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 S 3 ONC solution, we require that Q(β * ) ≤ Q(β Lasso ). However, this can generally be obtained by initializing any S 3 ONC 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 C 1 through C 7 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 min L(β) + j∈[p] P λ (|β k j |) • |β 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

β true S min > (a + 1)λ, λ > max{ 3λ Lasso s 1/2 b l re , 4σ √ s+2 √ st1+2t1 b l (an/bu) 1/2 , 2σ √ s+2 √ st2+2t2 b l √ 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 ) 2 nbua 8σ 2

),

φ 1 := P ( S c p n (β oracle ) max ≥ λ) ≤ ( pe s ) s exp(-t 1 ) + 2 exp( -λ 2 abun 8σ 2

),

φ 2 := P ( β oracle S min ≤ aλ) ≤ ( pe s ) s exp(-t 2 ),

(b) If in addition (A1) and (A2) holds, while the parameters of (a, λ) satisfy that 

P λ (aλ) > σ 2 2nb l (1 + 2 √ t 4 + 2t 4 ) + σ 2 n |S|(1+2 √ t4+2t4)b l b l ( p * -2|S|+1)

and t 1 , t 2 , t 3 , t 4 > 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.

4. NUMERICAL EXPERIMENTS:



4.1. EXPERIMENTAL SETUP:

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, β true S 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 β true S c are set to be 0. The design matrix X ∈ n×p is constructed by generating n iterations of x i ∼ N p (0, Σ) where Σ = (.5 |j-j | ) p×p . We then generate Y using a Bernoulli distribution where P (y i = 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 S 3 ONC 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 S 3 ON C 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/p 1/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 .

4.2. NUMERICAL RESULTS:



Table 1: Percent of time FCP beat all perturbations:

n = 100 n = 100 n = 100 n = 100 p = 500 p = 1000 p = 1500 p = 2000 % Best FCP 100% 100% 100% 100%

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.

5. CONCLUSIONS:

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.

A APPENDIX:

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 S 3 ONC solution β * within {β * : Q(β * ) ≤ Q(β true ) + Γ} for any Γ ≥ 0, we have the following useful inequality under Assumption (A1):

EQUATION

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 • b l (x i β * -x i β true ) 2 .

Proof of Theorem 1. First, given our assumption that (A1) holds, that (i

) p * ≥ 2|S|, (ii) β * is S 3 ON C satisfying Q(β * ) ≤ Q(β true ) + Γ for some Γ ≥ 0, and (iii) P λ (aλ) > σ 2 2nb l (1 + 2 √ t + 2t )+ σ 2 n |S|(1+2 √ t +2t )+Γb l b l ( 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,

1 n X(β * -β true ) 2 ≤ 8σ 2 b 2 l n p * + 2 √ p * t + 2t + 8 b l 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 (B c ))(1 -P (A c |B)) = 1 -P (A c |B) -P (B c ) + P (B c )P (A c |B) = 1 -P (A c |B) -P (B c )(1 -P (A c |B)) ≥ 1 -P (A c |B) -P (B c ) (9)

Therefor they hold simultaneous with probability at least 1 -exp(-t + p ln( pe p * )) -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 1 n X(β opt -β true ) 2 ≤ 8σ 2 b 2 l n p * + 2 √ p * t + 2t + 8 b l 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

( pe p * )) - 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

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)

.

Proof of Corollary 1. First we need to bound Γ. In order to do this we use the lasso problem

Q lasso (β) = i∈N (β, x i , y i ) + j∈P

λ lasso |β j | as well as the concavity of MCP over positive values to get the following 2 inequalities 

Q lasso (β lasso ) ≤ Q lasso (β true ) i∈N (β lasso j , x i , y i ) -(β true j , x i , y i ) ≤

|δ | ≤ 4|δ S | ≤ 4 √ s δ S 2 δ S ≤ 4 √ s δ 2 δ S ≤ 4 √ s r e n Xδ 2 δ S ≤ 4 √ s r e 3λ lasso |δ S | b l δ S ≤ 4 √ s r e 3λ lasso √ s δ S b l δ S (12)

which means |δ | ≤ 12λ lasso s b l re with conditional on A which occurs with probability at least 1 -

2p exp( -(λ lasso ) 2 nbua 8σ 2 )

Finally we are able to bound gamma by combining the above

EQUATION

EQUATION

EQUATION

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 β true S 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 n 1-γ/2 , where Q, > 0 are arbitrary constants, and setting p * = 4s, t = p * n γ/2 ln p, t = n γ/2 ln p we get the following

EQUATION

EQUATION

> 12 aQ + 20 + 12 2 aQ 2 s b l 2 1-γ = C 2 s b l 2 1-γ (20) β true S min > 8σ 2 r pb 2 l n p * + 2 p * t + 2t + 8 r pb l min{ λ 2 r p |S|, P λ (aλ)|S| + Γ} (21) n > (160 + 8Q 2 ) sσ 2 ln p ( β true S min r 4s b l r e ) 2 2/γ = C 3 sσ 2 ln p ( β true S min r 4s b l r e ) 2 2/γ (22)

For some constants C 1 , C 2 and C 3

We can then apply Theorem 1 (conditional on A) substitute our values and simplify to get that β * is the global solution with probability at least

EQUATION

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

EQUATION

For some constants C 4 , C 5 C 6 and C 7 .

Note that these constants, as well as C 1 , C 2 and C 3 , 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, C 1 through C 7 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 > a 0 λ). Similar to Lemma B.1. in Bickel et al. (2009) , to bound this we will start by noticing that for the lasso penalized loss function Q lasso (β) = i∈N l(β, x i , y i ) + λ lasso j∈P |β j | we have that Q lasso (β lasso ) ≤ Q lasso (β true ). If we then let δ = β lasso -β true we can use the same tactic as in the derivation of 8 to get b l (re) 1/2 conditional on A.

Next, using this we can show that conditional on A we have that

δ max ≤ δ ≤ Xδ 2 2 /( δ nr e ) ≤ 3λ lasso s 1/2 b l re < a 0 λ if λ > 3λ lasso s 1/2 b l a0re

. This is the inverse of the condition that defines φ 0 . Thus, we can bound φ 0 with φ 0 ≤ P (

A c ) = P ( j∈P | 1 n W X j | > λ lasso /2) = P ( j∈P |W X j |/ X j > nλ lasso /(2 X j )) ≤ pP (| W, v | > λ lasso n 2 Xj ) ≤ pP (| W, v | > λ lasso (nbua) 1/2 2 ) ≤ 2p exp -(λ lasso ) 2 nbua 8σ 2

which uses both (A1)(ii) and (A2) as long as λ > 3λ lasso s 1/2 b l a0(re) per (A2).

Next consider

EQUATION

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

EQUATION

We can then simplify both terms individually. For the first term, P (b u Xδ o ≥ (1/2)a 1 λ(ab u n) 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 

2 b l Ũ Sp W ≥ 2 b l σ s + 2 √ st 1 + 2t 1 ≤ ( pe s ) s exp(-t 1 ). This is the first half of φ 1 as long as (1/2)a 1 λ(an/b u ) 1/2 ≥ 2 b l σ s + 2 √ st + 2t which is equivalent to the as- sumed condition λ ≥ 4σ √ s+2 √ st1+2t1 b l a1(an/bu) 1/2

Next, the second term can be easily bounded using (A2): 

P (|W v| ≥ (1/2)a 1 λ(ab u n) 1/2 ) ≤ 2 exp(- a 2 1 λ 2 abun 8σ 2 ) Therefor φ 1 ≤ ( pe s ) s exp(-t 1 ) + 2 exp( -a 2 1 λ 2 abun 8σ 2 ) Next consider φ 2 = P ( β oracle S min ≤ aλ) First, given the assumption β true] min > (a + 1)λ we can see that φ 2 = P ( β oracle S min ≤ aλ) ≤ P ( β oracle S -β true S max > λ) ≤ P ( β oracle -β true 2 > λ) = P ( δ o 2 > λ).

Ũ Sp W > λ b l √ nre 2 ≥ σ s + 2 √ st 2 + 2t 2 ) ≤ ( pe s ) s exp(-t 2 ) assuming that λ b l √ nre 2 ≥ σ s + 2 √ st + 2t which is equivalent to the condition λ ≥ 2σ √ s+2 √ st2+2t2 b l √ nre

This, combined with the fact that for MCP, a 0 = a 1 = a 2 = 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 r 4s 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 S 3 ON C with Γ = 0. If we use t = t 3 and t = t 4 we get that the probability that the global solution is not the oracle solution as

φ 3 ≤ exp(-t 3 + p ln( pe p )) + exp(-(p * + 1)(t 4 -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 S 3 ON C 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

EQUATION

EQUATION

First, for any given j ∈ P, we solve for P [γ j ∩ δ j ] given our assumptions. We can start with

∂ 2 Q(β) (∂βj ) 2 β=β * ≥ 0 which gives us 1/n n i=1 ψ (x i β * ) x 2 i,j + P λ (|β * j |) ≥ 0. We can rearrange this to get b u n i=1 x 2 i,j ≥ n i=1 ψ (x i β * )x 2 i,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 b u X j 2 ≥ n/a which contradicts (A1) part (ii). Therefor we know P [γ j ∩ δ j ] = 0.

It should also be noted that P [γ c j ] = 0 since β * satisfies S 3 ON C conditions. Thus, by applying Demorgan's law and then the union bound, it can be obtained that

EQUATION

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. 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 V W 2 > σ 2 • (tr(Σ v ) + 2 tr(Σ v )t + 2 Σ v t)

+ tr(Σ v µµ ) • (1 + 2( Σ v 2 tr(Σ 2 v )

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.

:



Figures and Tables:

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.

aλ) and β k j -β k+1 j < δ. Otherwise, let k := k + 1 and go to Step 2. Remark 11. Here, the above algorithm has iteration complexity of O

2 results from the proof for φ 0 in Corollary 2 which shows that both |δ Sc | ≤ 3|δ S | and b l n Xδ 2 ≤ 3λ lasso |δ S | conditional on A where δ = β lasso -β true . Given our restricted eigenvalue assumption Xδ 2 n δ 2 ≥ r e , this can be used to show

let A = j∈P {| 1 n W X j | ≤ λ lasso /2}. We can combine this with 25 to get that b l 2n Xδ 2 + λ lasso /2 j∈P |β lasso j conditional on A. From this notice that the right term goes to zero when β true j = 0 so we then have that b l 2n Xδ 2 + λ lasso /2 j∈P |β lasso j Using the triangle inequality and the definition of δ we can simplify this tob l 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 bothb l n Xδ 2 ≤ 3λ lasso |δ S | ≤ 3λ lasso s 1/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 ≥ r e . If this is combined with the first of the two equations, we can get that 1 n 1/2 Xδ ≤ 3λ lasso s 1/2

Consider an arbitrary S 3 ON C solution β * to 2 with MCP. Given the event that for some integer p : β * -β true 0 ≤ p, then |W Xδ * | ≤ max Sp:|Sp|= p Ũ Sp W Xδ * , a.s.Where ( ŨSp ) i,j := U Sp , if j ∈ S p 0, else and U Sp ∈ n× p is defined as in the following Thin SVD:X Sp = U Sp D Sp V Sp .Proof. Denote δ * := (δ * j ) = β * -β true , S p := (j : δ * j = 0) ⊆ P , δ * Sp := (δ * j : j ∈ S p) and X Sp := (x ij : i ∈ N , j ∈ S p). By assumption, we know that δ * 0 ≤ |S p| = p. First decompose X Sp using Thin SVD to getX Sp = U Sp D Sp V Sp where U Sp ∈ n× p.Note that since and U Sp U Sp = I we have that for any υ ∈ p we haveD Sp V Sp υ 2 =If in addition we have the satisfaction of β opt -β true ≤ p and the eventE b (p) := { 1 n X(β opt -β true ) + 8 b l min{ λ 2 rp (|S| -β * 0 ), P λ (aλ) • (|S|β * 0 ) + Γ} then β * isboth 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 p (|S| -β * 0 ), P λ (aλ) • (|S| -β * 0 ) + Γ} ≥ 1 n Xδ * 2 ≥ r p δ * 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 r pb 2 l n p + 2 pt + 2t + 8 r pb l min{ λ 2 r p |S|, P λ (aλ)|S| + Γ} ≥ β * j -β true j surely. From this we can bound |β * j | using the square root term and |β true j + 8 rpb l min{ λ 2 rp |S|, P λ (aλ)|S| + Γ} > 0 then |β * j | > 0. From this we can obtain the inequality p |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 λ (aλ)β * S c 0 + Γ) ≥ 0 a.s. + 2t + Γ ≥ P λ (aλ) β * S c 0 a.s.(70)thus, it can be seen that if P λ (aλ) > σ 2 b l n p + 2 √ pt + 2t + Γ then 1 > β * S c 0 = 0. This is satisfied by the assumption that P λ (aλ) -σ 2 b l n p + 2√ pt + 2t > Γ.Finally, because β * is an S 3 ON C solution, it has to satisfy FONC. Per 2, this means that β * ∈ arg inf{ 1 n i∈N (β, x i , y i ) + 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 ∈ S c . Threfore we know that β * ∈ arg inf{ 1 n i∈N (β, x i , y i ) : β ∈ p , β j = 0, ∀j ∈ S c } a.s. (71)

Lasso converges to the global solution in two iterations with probability at least 1 -φ 0 -φ 1 -φ 2 -φ 3 where φ 3 := P (B oracle = B opt ) ≤ exp(-t 3 + p * ln(

Table content: <table><tr><td>pe p *  ))</td><td/></tr><tr><td>+ exp(-(p</td><td/></tr><tr><td>σ 2 b l n rpb 2 p *  + 2 8σ 2 l n p *  + 2 √ p *  t 3 + 2t 3 and let the minimal signal strength satisfy β true S √ p *  t 3 + 2t 3 + 8 rpb l min{ λ 2</td><td>min &gt;</td></tr></table>

Statistical performance of the four classifiers.

Table content: <table><tr><td/><td/><td colspan="2">n = 100, p = 1000</td><td colspan="2">n = 100, p = 1500</td><td colspan="2">n = 100, p = 2000</td></tr><tr><td colspan="2">Classifier Measure</td><td>Mean</td><td>Std. dev</td><td>Mean</td><td>Std. dev</td><td>Mean</td><td>Std. dev</td></tr><tr><td>MCP</td><td>1 loss</td><td colspan="6">13.909907 1.471911 14.818059 1.698191 14.506226 1.480686</td></tr><tr><td/><td>2 loss</td><td colspan="6">4.108019 0.320061 4.304993 0.374453 4.489184 0.399441</td></tr><tr><td>Lasso</td><td>1 loss</td><td colspan="6">15.015975 1.039529 15.882654 1.29422 17.079414 1.545309</td></tr><tr><td/><td>2 loss</td><td>4.3255</td><td>0.25996</td><td colspan="4">4.397969 0.326336 4.433467 0.362707</td></tr><tr><td>2 penalty</td><td>1 loss</td><td colspan="6">22.211963 0.791955 26.026067 0.966091 28.485075 0.993699</td></tr><tr><td/><td>2 loss</td><td colspan="6">4.734209 0.241683 4.738025 0.296726 4.755959 0.296746</td></tr></table>

b l 2n Xδ o 2 ≤ 1 n W Xδ o From here Lemma 2 can be applied since we know β oracle -β true 0 ≤ s. With some simplification this gives that Xδ o ≤ 2 b l (max Sp:|Sp|=s Ũ Sp W ) Utilizing Lemma 3 with s in place of p shows that P max Sp:|Sp|=s

Table content: <table/>

Next, since we know that the support of β oracle and β true is S, we know that |δ o Sc | = 0 ≤ 3|δ o S | which is the constraint for the RE condition. Therefor we know that Xδ o 2 n δ o 2 ≥ r e . With this and a similar line of argument as in φ 1 we get that φ 2 ≤ P ( δ

Table content: <table/>

Peer review (iclr2020):
=Summary= 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 . In the paper, the authors focus on the minimax concave penalty (MCP) that is a special class of FCP. Inspire of the problem being NP-hard, they proposed pseudo-polynomial time algorithms and showed the global optimality; using the characterization of a significant subspace second-order necessary condition (S^3ONC), they showed that all local solutions within an efficiently achievable sub-level set are globally optimal. They also investigates the empirical evaluation of proposed methods. =Significance= In the technical view, the part of (a) in Theorem 1 is very similar to Theorem 4 in [Hongcheng Liu, Tao Yao, Runze Li, and Yinyu Ye (Mathematical programming, 2017)]. The main result of this paper may be the part of (b) on Theorem 1. However, the motivation why the authors wish to compute global optimization is not clear to me. Indeed, for FCP, [Liu et al, 2017] showed that 1. Any local solution is a sparse estimator. 2. Any local solution satisfying a significant subspace second-order necessary condition yields a bounded error in approximating the true parameter with high probability 3. For MCP (with restricted eigenvalue condition), S^3ONC solution has oracle property. Therefore, the global optimality is not necessarily stipulated to ensure the recovery quality. I would like to confirm the motivation of this work. =Missing references= Although the topic addressed in the paper is related to the complexity of the sparse estimation problems, the following papers are not included in the references. Huo, X. and Chen, J. Complexity of penalized likelihood estimation. Journal of Statistical Computation and Simulation, 80(7):747–759, 2010. Bian, W. and Chen, X. Optimality conditions and complexity for non-lipschitz constrained optimization problems. Chen, X., Ge, D., Wang, Z., and Ye, Y. Complexity of unconstrained L2 ? Lp minimization. Mathematical Programming, 143(1-2):371–383, 2014. =Questions or comments= 1. Time complexity: Why and how the running time of the proposed algorithm can be polynomial time from pseudo-polynomial time? Where is the time complexity of the proposed algorithm explicitly presented in the paper? At least, the definitions of FPTAS and FPRAS and other technical terms of the complexity are missing in the paper. If the author wants to publish the paper in this conference, the definitions of them should be written in the paper, since many of readers in machine learning areas may not know the exact definitions. 2. Can the results presented in the paper extend to the case for SCAD? Since SCAD and MCP are shared with similar properties, I’d like to see more discussion on SCAD case. 3. About the problem setting: All the inputs and parameters are assumed to be rational number? Are there any assumption on the loss function? 4. After the definition 1, “some independence assumptions “ or “the same setting” are unclear. Please clarify the assumptions or setting to be self contained. =Minor comments= About the title: If the paper did not provide any results on SCAD, it is better to use the term *the minimax concave penalized* not *folded concave penalized* (Page 3: “The oracle solution is a hypothetical assumes the prior knowledge on the true support set…” is grammatically wrong?) Remark 1: remove one “gradient* Remark 4: S^3ONC References: Cun-Hui Zhang, Tong Zhang, et al. -> add *and*, remove *et al.* References: Cun-Hui Zhang et al. => remove et al. Page 18: delete the memo “ move to intro somehwere?”