Robust group and simultaneous inferences for high-dimensional single index model

Authors: weichao  yang
hongwei  shi
xu  guo
changliang  zou
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 03:35:49

Abstract: The high-dimensional single index model (SIM), which assumes that the response is independent of the predictors given a linear combination of predictors, has drawn attention due to its flexibility and interpretability, but its efficiency is adversely affected by outlying observations and heavy-tailed distributions. This paper introduces a robust procedure by recasting the SIM into a pseudo-linear model with transformed responses. It relaxes the distributional conditions on random errors from sub-Gaussian to more general distributions and thus it is robust with substantial efficiency gain for heavy-tailed random errors. Under this paradigm, we provide asymptotically honest group inference procedures based on the idea of orthogonalization, which enjoys the feature that it does not require the zero and nonzero coefficients to be well-separated. Asymptotic null distribution and bootstrap implementation are both established. Moreover, we develop a multiple testing procedure for determining if the individual coefficients are relevant simultaneously, and show that it is able to control the false discovery rate asymptotically. Numerical results indicate that the new procedures can be highly competitive among existing methods, especially for heavy-tailed errors.

Main paper: 
 introduction: 
 With fast development of information technologies, high-dimensional data are frequently collected in various scientific disciplines. To efficiently analyze high-dimensional data, various regularization methods were developed and thoroughly investigated in both theory and methodology. For a comprehensive review, see [19] . Though significant progress has been made about the theory of estimation and variable selection, statistical inference for a high-dimensional model remains largely unexplored until the seminal works of [58] , [47] and [29] . They developed debiased estimators that correct the bias introduced by LASSO and established the estimators' asymptotic normality. [40] proposed a decorrelated score method which is a general procedure and is applicable to a large family of penalized M-estimators. [45] introduced a recursive online-score estimation for high-dimensional generalized linear model. [46] proposed the nuisance penalized regression which does not penalize the parameters of interest. [21] and [20] studied inference problems for the high-dimensional Cox model and longitudinal data, respectively. Logistic regression model was investigated in [36] . Other recent developments include [44] , [50] , [5] , and [19] .Most of current inference procedures focus on parametric regression models. To address the potential mis-specification issue of parametric models, the single index model (SIM) which is one of the most popular semiparametric modelling techniques, has received extensive attention and in-depth research in the past decade. Let Y ∈ R be the response variable along with predictor vector X = (X 1 , • • • , X p ) ⊤ ∈ R p . We consider a general SIM,EQUATIONwhere β = (β 1 , • • • , β p ) ⊤ , the link function g(•) is unknown and ⊥ means independence. Equivalently, we haveY ⊥ X | X ⊤ β.The above statement is that, given X ⊤ β which is called index, the response Y is independent of the predictor vector X. Clearly, model (1.1) covers linear models, generalized linear models and also the classical SIM [26] , that is, Y = g 1 (X ⊤ β) + ϵ with g 1 (•) being an unknown function. It relaxes restrictive assumptions on parametric models and is flexible enough to capture complex relationship between the response and the predictors.For high-dimensional SIMs, variable selection has been considered by many authors. Examples include [30] , [60] , [57] , [1] , [51] , [42] , [52] , [41] and [43] . Our primary interest is to detect whether a set of predictors contributes to the response Y or not given the other predictors, that is testing the following group inference problem:H 0,G : β j = 0 for all j ∈ G versus H 1,G : β j ̸ = 0 for some j ∈ G,where G is a prespecified subset of {1, 2, . . . , p} with p 0 = |G|. Such a hypothesis naturally arises in the high-dimensional setting. For example, researchers may want to test whether a gene pathway, consisting of multiple genes for the same biological functions, is important for certain clinical outcome The above group inference problem also includes global significance testing as a special case when we set G = {1, 2, . . . , p} and p 0 = p. The group inference problem is more difficult than purely global significance testing. Actually, for group inference problem, we are dealing with a high-dimensional interested parametric vector with a high-dimensional nuisance parameter. Such a group inference problem was considered by [59] for high-dimensional linear models. See also [16] , and [23] . All of [59] , [16] , and [23] did not establish asymptotical honesty of their procedures, which is critical for reliable inference in high-dimensional setting.For SIM, due to the existence of the unknown link function g(•), the group inference problem is more difficult. Actually, the parameter of interest and the nonparametric function are bundled together, which makes the inference being challenging. For the model (1.1), [17] developed a debiased LASSO procedure for individual coefficient, while [53] considered simultaneous confidence interval for optimal treatment regimes under the classical SIM [26] . Their procedures, as well as many other existing works, rely heavily on the sub-Gaussian assumption about the error term. In practice, data with heavy-tailed distribution or in the presence of outliers are very common [18] , and the efficiency of those procedures developed under sub-Gaussian assumption would be largely deteriorated. Recently [25] investigated the robust inference problem of high-dimensional SIM by adopting the Huber loss with original response. However, their procedure still requires bounded first moment condition of the error term and needs to select a suitable robustification parameter. In this paper, we aim to develop robust group inference methods for high-dimensional SIM without any moment condition on the error term and without choosing robustification parameter.Interestingly, under a mild condition on the predictors, which is called linearity condition in the literature of sufficient dimension reduction [31] , we find that the SIM can be recast into a pseudo-linear model with transformed response, allowing us to make robust inference in a simple fashion. With a specific "response-distribution" transformation, the sub-Gaussian response assumption required in existing methods is totally avoided, and accordingly our proposed procedures are robust to outliers or heavy-tailed error distributions.Besides group inference of a prespecified subset of predictors, we are also interested in identifying relevant predictors. To this aim, large-scale simultaneous hypotheses are considered H 0j : β j = 0 versus H 1j : β j ̸ = 0, 1 ≤ j ≤ p.Apart from identifying as many nonzero β j as possible, to obtain results with uncertainty quantification, we would like to control the false discovery rate (FDR) which is an extremely popular tool to maintain the ability to reliably detect true alternatives without excessive false positive results when large-scale hypotheses are simultaneously tested [4] . [28] , [20] , [36] and [9] considered FDR control in high-dimensional regression models. We follow this line and develop a FDR control procedure for high-dimensional SIM.Our major contributions are listed from the following three aspects.• We extend the rank-LASSO procedure in [43] to include both convex and non-convex penalties and establish error bound of any local optimum of the empirical objective. These theoretical results are summarized in subsection 3.1. • In Section 2, we provide asymptotically honest group inference procedures based on the idea of orthogonalization for testing the joint effect of many predictors, which enjoys the feature that it does not require the zero and nonzero coefficients to be well-separated. We demonstrate the superiority of our test procedures both theoretically and empirically. Please see Section 3 for theoretical justification and Section 4 for numerical studies respectively. • We develop a multiple testing procedure for determining if the individual coefficients are relevant simultaneously, and show that it is able to control the FDR asymptotically. To this end, we develop suitable multiple testing procedures and show that the proposed methods can control the false discovery rate (FDR, [4] ) both theoretically and empirically. We refer the readers to Appendix A.3 and A.6 for details.Notation. For a d-dimension vector U, we write ∥U∥ r = ( d k=1 U r k ) 1/r and ∥U∥ ∞ = max 1≤k≤d |U k | to denote l r and l ∞ norms of U. Further we define ∥U∥ 0 = # {k : U k ̸ = 0}. A random variable X is sub-Gaussian if the moment generating function (MGF) of X 2 is bounded at some point, namely E exp(X 2 /K 2 ) ≤ 2, where K is a positive constant. A random variable Y is sub-Exponential if the MGF of |Y | is bounded at some point, namely E exp(|Y |/K ′ ) ≤ 2, where K ′ is a positive constant. For a, b ∈ R, we write a ∨ b = max{a, b}. 
group inference based on distribution transformation: 
 Without loss of generality, assume that E(X) = 0, and Σ = Var(X) > 0. Now let σ h = Cov{X, h(Y )} for a given transformation function h(•) of the response. Define β h = Σ -1 σ h . We then have the following result.Proposition 2.1. Assume that E(X | X ⊤ β) is a linear function of X ⊤ β. Then β h is proportional to β, that is, β h = κ h × β for some constant κ h .The above proposition follows directly from Theorem 2.1 in [32] . The assumption in above proposition is known as linearity condition (LC) for predictors. It is satisfied when X has an elliptical distribution and widely assumed in the sufficient dimension reduction literature [31, 14] . [24] showed that in high-dimensional setting, the LC holds to a reasonable approximation. Throughout of the paper, we assume that κ h ̸ = 0. This assumption is mild. In fact, when h(•) is monotone and g(•, •) is monotone with respect to the first argument, this assumption is satisfied.From the definition of β h , we can writeEQUATIONwhere the error term e must satisfy E(eX) = 0. Different from the existing literature which usually imposes independence between the regression error and the predictors, the regression error e in the transformed model is only uncorrelated with the predictors. This implies that under LC, we can recast the general SIM into a pseudo-linear model with transformed response h(Y ). Therefore, by Proposition 2.1, testing H 0,G is equivalent toH ′ 0,G : β hj = 0 for all j ∈ G versus H ′ 1,G : β hj ̸ = 0 for some j ∈ G.This reformulation is important, allowing us to circumvent the issue of estimating the unknown link function g(•, •) and to make inference of β hj , j ∈ G in a linear regression model with transformed response instead of β j in SIM. As we show later, this reformulation greatly facilitates the construction of our test statistic and simplifies the computation as well.In practice, we need to choose a suitable transformation function h(•). We note that [17] 's procedure also relies on the above proposition and they essentially work with h(Y ) = Y . For robustness consideration, throughout this paper, we consider the distribution function of Y , denoted by F (Y ) as the transformation function. Actually with the equation (2.3), given the widely imposed subgaussian assumption on the predictors, any bounded transformation function h(Y ) would lead the transformed error term e being subgaussian, even if the original error term ϵ in the single index model Y = g(X ⊤ β, ϵ) comes from Cauchy distribution. Further as noted by [43] , in the empirical distribution function, the termn j=1 I(Y j ≤ Y i ) is the rank of Y i .Since statistics with ranks such as Wilcoxon test and the Kruskall-Wallis ANOVA test, are well-known to be robust, this then intuitively explains why our procedures with response-distribution transformation are robust with respect to outliers in response. Moreover the distribution function is very easy to estimate and thus our approach is straightforward to implement and understand.Our test statistic for group inference of H ′ 0,G relies on individual inference of H ′ 0j : β hj = 0. Now we first consider individual hypothesis H ′ 0j . Let Z j be the subvector of X without X j , and γ j be the subvector of β h without β hj . Suppose that {X i , Y i } n i=1 is a random sample from the population (X, Y ). Similarly, we denote Z ij as the sample of Z j . Defineβ h = E(XX ⊤ ) -1 E[X{F (Y ) -1/2}] and θ j = E(Z j Z ⊤ j ) -1 E(Z j X j ), where β h is the regression coefficient of the model (2.3) with h = F .Our approach is based on the idea of orthogonalization. The main idea of orthogonalization is to construct a statistic for target parameter which is locally insensitive to the nuisance parameters. Dealing with high-dimensional models, it plays an important role to make the statistic of target parameter immune to the bias from the estimators of high-dimensional nuisance parameters, which in turn enables the statistical inference of parameter of interest. For relevant references, see for example, [3] , [40] and [2] .As pointed by [53] , the adoption of orthogonalization is nontrivial for high-dimensional semiparametric setting, particularly for index model, where the challenge of bundled parameter arises. Fortunately, with the Proposition 2.1, the SIM can be recast into a pseudo-linear model with transformed response (2.3) . Note that under the null hypothesisH 0j E F (Y ) -1/2 -Z ⊤ j γ j X j -Z ⊤ j θ j = 0.∂ ∂γ j E F (Y ) -1/2 -Z ⊤ j γ j X j -Z ⊤ j θ j = E Z j X j -Z ⊤ j θ j = 0.This then motivates us to consider the following quantityT * nj = 1 √ n n i=1 F (Y i ) -1/2 -Z ⊤ ij γ j X ij -Z ⊤ ij θ j . (2.4)In practice, we need to use suitable estimates of F , γ j and θ j in (2.4) as those quantities are unknown. Naturally, F n (y) = n -1 n i=1 I(Y i ≤ y), the empirical distribution of Y 1 , . . . , Y n , can be used to estimate F (Y ). For θ j , we estimate it by the penalized least-squares method θj = arg minEQUATIONwhere p λ X (•) is a penalty function with a tuning parameter λ X . Similarly, for γ j , we adopt the following penalized least-squares βh = arg minEQUATIONwhereL n (β h ) = (2n) -1 n i=1 {F n (Y i ) -1/2 -X ⊤ i β h } 2 and p λ Y (•)is a penalty function with a tuning parameter λ Y . The γj is then set as the subvector of βh without βhj . Accordingly, for each individual hypothesis H ′ 0j : β hj = 0 we define the standardized test statisticEQUATIONwhereEQUATION)êij = F n (Y i ) -1/2 -Z ⊤ ij γj and mj (y) = n -1 n i=1 (X ij -Z ⊤ ij θj ){I(Y i ≥ y) -F n (Y i )}. Note that F n (Y i ) = n -1 n j=1 I(Y j ≤ Y i ).Then given predictors X i 's being fixed, perturbations in the responses would not make the value of T nj change as long as the ranks of Y i 's remain unchanged.Denote T n,G = ( T nj ) j∈G . To test the null hypothesis H 0,G , we consider test statistic based on the max norm of T n,G . That is,M n,G = max j∈G T 2 nj .(2.9)Based on the limiting null distribution obtained in subsection 3.2, we can reject null hypothesis H 0,G at the significant level α if and only if M n,G ≥ c G (α), where c G (α) = 2 log p 0 -log log p 0 + q α and q α is the 1 -α quantile of the Gumbel distribution with the cumulative distribution function exp{-EQUATIONOur inference procedure is summrized in the Algothrim 1 as follows. n,G defined in (2.9) ; if M n,G ≥ c G (α), reject H 0,G ; otherwise accept H 0,G .Remark 2.1. The use of quantiles of limiting null distribution in c G (α) is attractive from a computational point of view. On the other hand, the validity of limiting null distribution requires additional assumptions regarding the dependence structure of the components of the covariates X. See Assumption A.5 in Appendix A.8 for details. Besides, it is well known that this weak convergence is typically slow. To solve these problems, we propose a multiplier bootstrap approach and show its validity in theory, please see Appendix A.1 for details. 
estimation error bound: 
 The theoretical analysis of βh requires a substantial modification of the proof technique as compared to existing works on high-dimensional inference. It is related to the fact that empirical distribution F n (Y i ) are dependent, and thus L n (β h ) is a sum of dependent random variables. Note that [43] considered a similar penalized procedure with LASSO penalty. While [60] and [51] considered nonconvex penalties such as SCAD and MCP to reduce estimation bias incurred by LASSO. However, their results only allow the dimension to be polynomial order of the sample size and thus cannot work in ultrahigh-dimensional setting. Further their results only provide guarantees for global optima. In this paper, we consider both convex and non-convex penalties and establish error bound of any local optimum of the empirical objective.To be specific, assume that βh satisfies the first-order necessary condition to be a local minimum of 2.6, that is, EQUATION∥ βh -β h ∥ 2 ≤ cλ Y √ s Y and ∥ βh -β h ∥ 1 ≤ c ′ λ Y s Y with probability at least 1 -c 1 exp(-c 2 log p).Here, (c, c ′ , c 1 , c 2 ) are universal constants and λ X , λ Y ≍ log p/n.In Theorem 3.1 we establish error bound of penalized least-squares estimators with empirical distribution function of the response. For consistency in L 2 -loss, we require the sparsity level s Y satisfy that s Y = o(n/ log p). While in terms of L 1 -loss, the limitation becomes s Y = o( n/ log p). Clearly the sparsity level is allowed to be diverging. Our results unify both convex [43] and nonconvex penalties [60, 51] . Compared with the results in [60] and [51] , the dimension can be exponential order of the sample size, and no minimal signal condition is imposed to obtain those error bounds. In the proof, we modify Hoeffding's inequality by a probability inequality for U -statistic to handle tail probability in dependent case, which may be interesting in their own rights.In this subsection, we derive the asymptotic null distribution for statistic defined in (2.9). Denote Ω = β h ∈ R p : ∥β h ∥ 0 ∨ max 1≤j≤p ∥θ j ∥ 0 ≤ s and β hG = (β hj ) j∈G . We consider the following parameter space for H 0,G Ω 0 G = {β h ∈ R p : β hG = 0} ∩ Ω.β h ∈Ω 0 G Pr M n,G -2 log p 0 + log log p 0 ≤ t -exp - 1 √ π exp(-t/2) = 0.Theorem 3.2 implies that the type I error of the proposed test statistic M n,G converges to any pre-specified significance level uniformly over β h ∈ Ω 0 G . The hypothesis test with such uniform convergence property is called honest test. The advantage of honest test is that the limiting distribution of our procedure is uniformly valid over s-sparse high-dimensional models despite the possible imperfect model selection via estimator βh . An immediate implication is that it relaxes the assumption on signal strength and does not require the zero and nonzero effects to be well-separated. In particular, this procedure does not require the initial estimator to select zero and nonzero signals perfectly, which is nearly impossible in practice. Due to these excellent statistical properties, honest test has recently drawn lots of attention. See [47] , [29] , [3] , [2] , [11] , [53] for further examples.On the basis of Theorem 3.2 , we can reject null hypothesis H 0,G at the significant level α if and only if M n,G ≥ c G (α), where c G (α) = 2 log p 0 -log log p 0 + q α and q α is defined in (2.10).We next to consider the asymptotic power analysis of the M n,G . In this section, we show that our test statistic is powerful in the sense that the separation rate is of order (2 + ϵ 0 ) log p 0 /n for some ϵ 0 > 0 when p 0 → ∞. At the beginning we define the following parameter space for H 1,GΩ 1 G (c 0 ) = β h ∈ R p : max j∈G β hj δ j σ j ≥ c 0 log p 0 n ∩ Ω,where c 0 is a positive constant andδ j = E X 2 j -E X j Z ⊤ j E Z j Z ⊤ j -1E Z j X j . Theorem 3.3. Suppose that conditions in Theorem 3.2 are satisfied, we have for some ϵ 0 > 0, lim(n,p0)→∞ inf β h ∈Ω 1 G (2+ϵ0) Pr M n,G ≥ c G (α) = 1.(3.12) Theorem 3.3 implies that our test procedure can still be powerful even when there exists few components of β G with a magnitude being larger than (2 + ϵ 0 ) log p 0 /n. Therefore our testing procedure is powerful against "sparse" alternative aforementioned. This separation rate is widely discussed in the literature, such as [8] , [59] and [36] . Theorem 3.3 shows our tests based on distribution transformation can also achieve this lower bound. Remark 3.1. In the proof of Theorem 3.3, we show thatM n,G ≍ M 0 G + max j∈G β hj δ j σ j 2 ,where the M 0 G has the type I extreme limiting distribution, which is mentioned in Theorem 3.2. Then the power is largely determined by the second termEQUATIONand it can be explained as the maximum of the signal-to-noise ratio for j ∈ G. The power of the test based on M n,G increases with the growth of SNR(M n,G ). When SNR(M n,G ) ≥ (2 + ϵ 0 ) log p 0 /n, the asymptotic power is tending to 1.In this paper, we focus on the distribution transformation of the response. In the following we make some comparisons between this choice and another natural one. Actually, same to [17] , one may consider use h(Y ) = Y to conduct the following test statisticEQUATIONHere γ+ j is the subvector of β+ h without j-th component, and β+ h is a penalized least-squares estimator of β + h = E(XX ⊤ ) -1 E(XY ) as follows:β+ h = arg min β + h ∈R p 1 2n n i=1 (Y i -X ⊤ i β + h ) 2 + p l=1 p λ Y (|β + hl |).θj is estimated by (2.5) . And σ+2 j is an appropriate estimator of σ +2 j :=E{(X j -Z ⊤ j θ j ) 2 (Y -Z ⊤ j γ + j ) 2 }. For the group inference problem (1.2), it's natural to consider the following testing statistic:M + n,G = max j∈G ( T + nj ) 2 ,Now we discuss the power properties of M + n,G . By similar arguments as Theorem 3.3, the power of the test based on M + n,G is determined bySNR(M + n,G ) = max j∈G β + hj δ j σ + j 2 . (3.15)For the convenience of illustration, we consider the classical SIM, that is, Y = g 1 (X ⊤ β) + ϵ and ϵ is independent of X. To compare the theoretical power of M n,G and M + n,G , it suffices to compare the SNRs aforementioned. In the classical SIM case, the ratio of SNR can be simplified asSNR(M n,G ) SNR(M + n,G ) ≥ c 3 Var(ϵ),where c 3 is a constant depending on (X, Y ). The detailed form of c 3 is presented in the Appendix A.2. When Var(ϵ) > c -1 3 , M n,G is more powerful than M + n,G . Actually when ϵ follows a heavy-tail distribution, Var(ϵ) can be very large even be infinite.The above results illustrate the robustness of our test statistic M n,G . In fact, with distribution transformation, no moment assumption on ϵ is required. Even if the ϵ comes from Cauchy distribution, our test procedure still works well but the M + n,G requires sub-Gaussian assumption on the error term and thus in this situation would fail.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Answer: [Yes]• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
numerical studies: 
 In this section, extensive simulation studies are carried out to evaluate the numerical performance of the proposed methods for the global inference problem described in (1.2). We consider the high-dimensional SIM in equation (1.1) and generate the data from the following models:Model 1: Linear model: Y = X ⊤ β + ϵ. Model 2: Non-linear model: Y = exp(X ⊤ β + ϵ).Here, the predictors X i , i = 1, • • • , n are generated from multivariate normal distribution N p (0 p , Σ).The covariance matrix Σ ∈ R p×p is block diagonal, that is, Σ = diag (Σ 1 , • • • , Σ 10 ), where Σ 1 is p 10 × p10 dimensional identity matrix, and each Σ k is p 10 × p 10 dimensional with an AR(1) correlation structure, that is(Σ k ) ij = (0.1k -0.1) |i-j| , k = 2, . . . , 10; i, j = 1, • • • , p/10.We consider two different error distributions for ϵ which is independent of X: (1) standard normal distribution N (0, 1); (2) Cauchy distribution or equivalently Student's t distribution with 1 degree of freedom, t(1). Without loss of generality, we set the active set be {1, 2, . . . , 6}. The regression coefficients β are generated from an arithmetic sequence from 0.1 to 2, that is β j = 0.2 + 0.038(j -1) for 1 ≤ j ≤ 6 and β j = 0 otherwise. We consider n = 200, 500 and p = 800. We use Intel(R) Xeon(R) Silver 4208 CPU @ 2.10GHz.In order to explore the robustness behaviors of our proposed statistics, we add outliers to pollute the observations: p out of the responses are picked at random and increased by m out -times maximum of original responses, shorted as p out + m out • max(responses). Here p out is the proportion of outliers and m out is a pre-set constant standing for outlier strength. Throughout the simulation study, we analyze the results with both the original data (shorted as Ori.) and the data with outliers (shorted as Out.).We assess the empirical type I error and the empirical power of group inference with the nominal level α = 0.05 based on 500 simulation runs. We consider the following five different choices for G to test the hypothesis (1.2):G 1 = {10, 11, p -4, p -3, p -2, p -1, p}, G 2 = {3, . . . , 6} ∪ G 1 , G 3 = {10, . . . , 59, p -149, . . . , p}, G 4 = {3, . . . , 6} ∪ G 3 .Note that G 1 and G 2 are small groups, G 3 and G 4 are large groups. G 1 and G 3 consist of only zero coefficients, while G 2 and G 4 includes nonzero elements. Then we report the numerical results of empirical rejection rate (ERR) for different scenarios, where ERR is the proportion of rejected hypotheses among the total 500 simulations. For G 1 and G 3 , the ERR is the empirical type I error; For G 2 and G 4 , the ERR is the empirical power.From the numerical results for all scenarios displayed in table 1, we have the following findings. Firstly, when there are outliers or heavy errors, our methods still have very high powers and control empirical sizes well. Secondly, for both linear model and nonlinear model, small groups and large groups, our procedures have similar satisfactory performance in this experiment. Thirdly, we find that the dimensionality of all predictors p does not effect the empirical performances. Our test procedures with LASSO penalty perform similarly to those with SCAD or MCP. Therefore, we only present the results of test statistics with LASSO penalty. We use the R-package ncvreg [6] .Moreover, we compare the proposed methods with the three-step testing procedure based on the studentized statistics in [59] , denoted as ST. Here, we use R package SILM [59] to implement the three-step testing procedure, set the number of bootstrap replications as 500 and choose the splitting proportion of 30% for screening.Table 2 summarizes the results of empirical type I errors and powers with α = 0.05 for different methods and models based on 500 simulation runs (due to the computation limit, we only run 200 times for ST). Here we consider (n, p) = (200, 400), and other settings are same as before. It is obvious that our methods outperform the ST. Firstly, when there are no outliers and heavy errors, all methods have empirical power 1 for linear model. But for non-linear model, our test procedures T (1) α and T(2) α are more powerful than ST. For instance, under non-linear model with original data and ϵ ∼ N (0, 1), the empirical powers of ST for G 2 and G 4 are 0.785 and 0.430, while the corresponding Table 1 : Simulation results for the group inference problemsϵ ∼ N (0, 1) ϵ ∼ t(1) G1 G2 G3 G4 G1 G2 G3 G4 n Ori. Out. Ori. Out. Ori. Out. Ori. Out. Ori. Ori.Ori.Ori.Y = X ⊤ β + ϵ 200 0.028 0.034 1.000 1.000 0.026 0.046 1.000 1.000 0.026 1.000 0.058 1.000 500 0.018 0.028 1.000 1.000 0.024 0.040 1.000 1.000 0.032 1.000 0.054 1.000 Y = exp(X ⊤ β + ϵ) 200 0.018 0.024 1.000 1.000 0.030 0.052 1.000 1.000 0.044 1.000 0.084 1.000 500 0.034 0.040 1.000 1.000 0.026 0.040 1.000 1.000 0.030 1.000 0.028 1.000 empirical powers of our methods are all 1. Secondly, the outliers and heavy errors have dramatic effect on the power performances of ST, which has very low powers even as low as the nominal level. For example, under linear model with outliers, the empirical power of ST for G 4 is 0.050 and that of our proposal is 1; under non-liner model with Cauchy error, the corresponding empirical powers of ST and our proposal are 0.040 and 1, respectively. Thirdly, all methods can control type I errors reasonably. Overall our methods have better powers and are robust with respect to outliers and heavy errors, whereas ST in [59] can completely fail when outliers exist.Table 2 : Simulation results for our method and ST.ϵ ∼ N (0, 1) ϵ ∼ t(1) G1 G2 G3 G4 G1 G2 G3 G4Method Ori. Out. Ori. Out. Ori. Out. Ori. Out. Ori. Ori.Ori.Ori.Y = X ⊤ β + ϵ Our 0.028 0.030 1.000 1.000 0.038 0.044 1.000 1.000 0.012 1.000 0.050 1.000 ST 0.045 0.040 1.000 0.115 0.015 0.040 1.000 0.050 0.055 0.370 0.055 0.250 Y = exp(X ⊤ β + ϵ) Our 0.032 0.054 1.000 1.000 0.034 0.070 1.000 1.000 0.042 1.000 0.054 1.000 ST 0.025 0.040 0.785 0.045 0.060 0.075 0.430 0.080 0.045 0.085 0.035 0.040 
discussions: 
 In this paper, we investigate robust inference for high-dimensional single index model (SIM). Under the linearity condition about the predictors, we recast the SIM into pseudo-linear model with transformed response. A response-distribution transformation is considered. This transformation choice avoids the sub-Gaussian assumption for the response. Our introduced procedures are thus robust with respect to outliers or heavy-tailed errors. We develop asymptotically honest group inference procedures. Asymptotic distribution and bootstrap implementation are both established. For testing the individual coefficients simultaneously, multiple testing procedures are proposed and shown to control the false discovery rate asymptotically. Our numerical studies illustrate the robustness of our procedures against outliers or heavy errors.In this paper, the linearity condition about the predictors plays a very important role in the methodology development. How to make inference for high-dimensional SIM without this condition is of great importance and interest. The second-order Stein's method with score function-based corrections investigated in [56] and [55] would be a powerful alternative. The condition κ h ̸ = 0 excludes even link functions, and in particular the problem of sparse phase retrieval. [39] introduced a novel procedure to deal with the problem of sparse phase retrieval when κ h = 0. It would be of interest to extend their approach to make robust group inference. Further our paper doesn't provide quantitative bounds on robustness. It is also of interest to make inference for partially linear single-index regression model [15] . We will pursue these challenging problems in near future.A Appendix / supplemental material Notation. For functions f (n) and g(n), we writef (n) ≲ g(n) to mean that f (n) ≤ cg(n) for some universal constant c ∈ (0, ∞), and similarly, f (n) ≳ g(n) when f (n) ≥ c ′ g(n) for some universal constant c ′ ∈ (0, ∞). We write f (n) ≍ g(n) when f (n) ≲ g(n) and f (n) ≳ g(n) hold simultaneously.For a function h : R p → R, we write ∇h and ∇ k h to denote a gradient or subgradient and its k-th component, if they exist. The sub-Gaussian norm of a sub-Gaussian random variable X is defined as∥X∥ ψ2 := inf{t > 0 : E exp(X 2 /t 2 ) ≤ 2}. Similarly, the sub-Exponential norm of a sub-Exponential random variable Y is defined as ∥Y ∥ ψ1 := inf{t > 0 : E exp(|Y |/t) ≤ 2}.As discussed in Remark 2.1, the critical value obtained from the asymptotic distribution may not work well in practice since this weak convergence is typically slow and assumptions for dependence structure of the covariates. We next describe a simple gaussian multiplier bootstrap method to obtain an accurate critical value. Let w 1 , . . . , w n be i.i.d. N (0, 1) random variables that are independent of {X i , Y i } n i=1 and define the gaussian multiplier bootstrap statistic:EQUATIONĥj (Y i , X i ; Y ĩ, X ĩ) = I(Y ĩ ≤ Y i ) - 1 2 -Z ⊤ ij γj X ij -Z ⊤ ij θjis the kernel funtion in T nj . The bootstrap critical value is given by c(2) G (α) = inf t ∈ R : Pr w M ♯ n,G ≤ t ≥ 1 -α ,EQUATIONM ♯ n,G -L ♯ n,G ≥ ζ 1 → 0 in probability, where L ♯ n,G is the ideal version of M ♯ n,Gand it is defined in (1.27) in the Appendix A.9. Theorem A.2. Suppose that conditions in Theorem A.1 are satisfied, we have for some ϵ 0 > 0, lim(n,p0)→∞ inf β h ∈Ω 1 G (2+ϵ0) Pr M n,G ≥ c G (α) = 1.(1.17)WhereΩ 1 G (•) is defined in subsection 3.3.The theoretical guarantee of the bootstrap procedure is based on the gaussian approximation theory for nondegenerate U -statistic. Symmetrization of the kernel function ĥj (•) plays an important role in theory. Thus ĥsym j (•) is necessary in our gaussian multiplier bootstrap procedure. Meanwhile, Theorem A.1 holds whether p 0 is divergent or not. While the limiting null distribution established in Theorem 3.2 asks p 0 be divergent.By the definition of SNR and Lemma A.8 in the Supplementary Material, for a given β ∈ R p we haveSNR(M n,G ) = Cov 2 (F (Y ), X ⊤ β) (β ⊤ Σβ) 2 max j∈G β j δ j σ j 2 .Similarly, we haveSNR(M + n,G ) = Cov 2 (Y, X ⊤ β) (β ⊤ Σβ) 2 max j∈G β j δ j σ + j 2 .Thus the ratio of SNR(M n,G ) and SNR(M + n,G ) can be simplified asSNR(M n,G ) SNR(M + n,G ) ≥ Cov 2 (F (Y ), X ⊤ β) Cov 2 (Y, X ⊤ β) • (σ + j ′ ) 2 σ 2 j ′, where j ′ = arg j∈G max βj δj σ + j 2 . In the case of classical SIM, that is, Y = g 1 (X ⊤ β) + ϵ and ϵ is independent of X. We can derive(σ + j ′ ) 2 = E{(X j ′ -Z ⊤ j ′ θ j ′ ) 2 (Y -Z ⊤ j ′ γ j ′ ) 2 } ≥ E(X j ′ -Z ⊤ j ′ θ j ′ ) 2 Var(ϵ).Thus we haveSNR(M n,G ) SNR(M + n,G ) ≥ Cov 2 (F (Y ), β ⊤ X)E(X j ′ -Z ⊤ j ′ θ j ′ ) 2 Cov 2 (Y, β ⊤ X)σ 2 j ′ Var(ϵ),andc 3 = Cov 2 (F (Y ), β ⊤ X)E(X j ′ -Z ⊤ j ′ θ j ′ ) 2 Cov 2 (Y, β ⊤ X)σ 2 j ′ .In the previous section, we consider joint significance testing of a prespecified subset of predictors. In applications, however, the parameter of interest may not be specified in advance. Denote H 0 = {j : β j = 0, j = 1, . . . , p}, and H 1 = {j : β j ̸ = 0, j = 1, . . . , p}. It is of interest to identify the elements in H 1 . To this aim, we consider simultaneous testing of the following hypothesesH 0j : β j = 0 versus H 1j : β j ̸ = 0, 1 ≤ j ≤ p.In this section, we aim to develop data-driven procedures to control the FDR.Recall that in Section 2, for each individual hypothesis H 0j : β j = 0 we define the standardized statisticT nj = σ-1 j T nj defined in (2.7). At a given threshold level t > 0, H 0j is rejected if | T nj | ≥ t. For each t, let R 0 (t) = j∈H0 I | T nj | ≥ t , and R(t) = p j=1 I | T nj | ≥ tbe the total number of false discoveries and the total number of discoveries, respectively. Accordingly, the false discovery proportion (FDP) and FDR areFDP(t) = j∈H0 I | T nj | ≥ t max p j=1 I | T nj | ≥ t , 1 , FDR(t) = E{FDP(t)}.Let G 0 (t) be the proportion of the nulls falsely rejected by the procedure among all the true nulls at the threshold level t, namely, G 0 (t) = q -1 0 j∈H0 I | T nj | ≥ t , where q 0 = |H 0 |. In practice, it is reasonable to assume that the true alternatives are sparse, that is, p -q 0 =: q 1 = o(p). If the sample size is large, we can use the tails of normal distribution G(t) = 2 -2Φ(t) to approximate G 0 (t). In fact, it will be shown that, for b p = √ 2 log p -log log p,sup 0≤t≤bp G 0 (t) G(t) -1 → 0 in probability as (n, p) → ∞.To summarize, we have the following multiple testing procedure controlling the FDR and FDP at a pre-specified level 0 < α < 1, which is summarized in Algorithm A.1.Algorithm 2: multiple testing procedure EQUATIONPr(| T nj | ≥ t) 2 -2Φ(t) -1 → 0. (1.19)If in addition Assumption A.9 holds, thensup 0≤t≤bp G 0 (t) G(t) -1 → 0 (1.20)in probability.In the results given in Section 2, we could allow p to grow exponentially fast in n. However, in order to control FDR, we only allow p = O(n c ) to apply the moderate derivation result in [34] . Further we assume s = o( √ n/(log p) 2 ) which is slightly stricter than s = o( √ n/(log p log p 0 )) to make the difference between T nj and its ideal version S nj uniformly bounded with the rate o(1/ √ log p). Note that [28] imposed a similar condition. The following theorem provides the asymptotic FDR control of our procedure. Theorem A.4. Under the assumptions in Theorem A.3, for t defined in our multiple testing procedure, we have lim sup(n,p)→∞ FDR( t) ≤ α.(1.21)In this subsection, we carefully discuss the robustness of our inference procedure based on the efficient influence function. We find that our statistics are robust with respect to perturbations in the responses. The details are given below. Recall that our test procedure is inspired by the quantityI := E[{F (Y ) -1/2 -Z ⊤ j γ j }(X j -Z ⊤ j θ j )].Next we derive the efficient influence function (EIF) of I. We rewrite I asI = Ψ(P) = E P [{F P (Y ) -1/2 -Z ⊤ j γ j }(X j -Z ⊤ j θ j )],where P is distribution of (X j , Z ⊤ j , Y ) ⊤ . Consider the following parametric submodel indexed by t, i.e.P t = t P + (1 -t)P,where t ∈ [0, 1], and P is a point mass at a single observation õ = (x j , z⊤ j , ỹ) ⊤ . As mentioned in [27] , the EIF for I at observation õ isϕ(õ, P) = dΨ(P t ) dt t=0 ,whereΨ(P t ) = E Pt [{F Pt (Y ) -1/2 -Z ⊤ j γ j }(X j -Z ⊤ j θ j )]. Denote o = (x j , z ⊤ j , y) ⊤ and O = (X j , Z ⊤ j , Y ) ⊤ . Let m t (o) = {F Pt (y) -1/2 -z ⊤ j γ j }(x j -z ⊤ j θ j ) and Ψ(P t ) = E Pt [m t (O)].Further, the operator, ∂ t , applied to an arbitrary function g(t), is defined as∂ t g(t) = dg(t) dt t=0 .Some calculation entails that∂ t Ψ(P t ) = ∂ t m t (O)dP t | t=0 + m t (O)| t=0 d∂ t P t = E P [∂ t m t (O)] + m 0 (õ) -Ψ(P).Note thatEQUATIONE P [∂ t m t (O)] = E P [{I(Y ≥ ỹ) -F P (Y )}(X j -Z ⊤ j θ j )].As m 0 (õ) = {F P (ỹ) -1/2 -z⊤ j γ j }(x j -z⊤ j θ j ), by some calculations we deriveϕ(õ, P) = {F P (ỹ) -1/2 -z⊤ j γ j }(x j -z⊤ j θ j ) + E P [{I(Y ≥ ỹ) -F P (Y )}(X j -Z ⊤ j θ j )] -Ψ(P). Recall that T nj = 1 √ n n i=1 F n (Y i ) - 1 2 -Z ⊤ ij γj X ij -Z ⊤ ij θj , S nj = 1 √ n n i=1 F (Y i ) - 1 2 -Z ⊤ ij γ j X ij -Z ⊤ ij θ j + m j (Y i ) ,wherem j (y) = E[(X j -Z ⊤ j θ j ){I(Y ≥ y) -F (Y )}].Here T nj is an individual test statistic corresponding to the parameter I. While S nj is an individual quantity corresponding to the EIF of I. From our proof, we know that T nj = S nj + o p (1) .By the formula of ϕ(õ, P), given (x j , z ⊤ j ) ⊤ , it can shown that ϕ(õ, P) is bounded for any ỹ ∈ R. Thus our test statistics are robust with respect to the perturbations in the responses.In this subsection, we conduct simulation studies to compare our procedure with other methods based on different transformation functions. The simulation settings are summarized as follows. We consider the following two models:• Model 1: Linear model: Y = X ⊤ β + ϵ. • Model 2: Non-linear model: Y = exp(X ⊤ β + ϵ).The regression coefficients β = (β 1 , β 2 , . . . , β p ) ⊤ are generated as β j = δ for j = 1, . . . , 6 and β j = 0 otherwise, where δ can be regarded as a signal strength parameter. We generate the error term from the standard normal distribution. We add outliers to pollute the observations: p out of the responses are picked at random and increased by m out -times maximum of original responses. Specifically, the detailed settings for the above parameters are as follows. Firstly, we consider the sample size n = 500 and the dimension p = 800. Secondly, we set the signal strength parameter δ to vary from {0.1, 0.3, 0.5}. Thirdly, we fix m out to 10. Lastly, we vary p out from 0 to 0.5 in increments of 0.1. For more complete comparisons, we consider three transformation procedures: ( 1)h(Y ) = F (Y ) (Our method); (2) h(Y ) = Y ; (3) h(Y ) = sigmoid(Y ) = 1/{1 + exp(-Y )}.Other simulation settings are the same as described in section 4 of the main text.Figure 1 summarizes the results of empirical type I error and empirical power with the significant level of α = 0.05 for different methods when the error term follows the standard normal distribution. Firstly, under the null hypothesis (G 1 ), h(Y ) = sigmoid(Y ) cannot control the type I error for Model 2, while other procedures perform well. Secondly, under the alternative hypothesis (G 2 ), the empirical powers of other procedures decrease rapidly with the increase of p out , while the powers of our procedure remain stable, which is particularly noticeable when δ = 0.5. This finding indicates that our method has strong robustness when the responses are polluted. Thirdly, our method performs well for both Model 1 and Model 2, indicating that our method is robust across different single-index models. Figure 1 : Under different settings of the generated model, the testing group and the signal strength, simulated results of the proposed method for the proportion of outliers p out from 0 to 0.5 in increments of 0.1 when the error term follows the standard normal distribution.In this subsection, we examine the empirical FDR results of our proposed procedure. We also make comparisons with two alternatives. Let BH-V n denote the BH procedure [4] with p-values obtained from V n . BH-W n is similarly defined. Here V n and W n denote the debiased lasso procedure introduced by [17] and the decorrelated score method proposed by [40] , respectively. Here, we use the R function p.adjust(...,method="BH") to implement the BH procedure. We simultaneously test H 0j : j ∈ H 0 for all j = 1, • • • , p. We set the dimension p = 1000 and sample size n = 1000, 1500.The sparsity level |H 1 | = q 1 is set as 4 or 8. The data is generated from the following linear and nonlinear models:Y = X ⊤ β + ϵ, Y = sin(0.5X ⊤ β) • exp(X ⊤ β + 1) + ϵ, where β = β 1 with δ 0 = 1, that is β = ( s0 1, • • • , 1, 0, • • • , 0)⊤ and ϵ comes from N (0, 1) or t(1).We evaluate the empirical FDRs and the empirical powers of our test statistic T n and two existing high-dimensional inference procedures BH-V n and BH-W n . Here the power is defined as the number of correctly discovered variables divided by the number of truly active variables. We set the desired FDR level as α = 0.2 and the results are calculated based on 500 replications.We illustrate our proposed procedure under different models, error distributions, sample sizes, and sparsity levels. Table 3 shows the simulation results. We have the following important observations. Firstly our introduced procedure controls the FDR well, even when the random error comes from Cauchy distribution. Secondly, our method has power 1, indicating that the proposed method can detect all important predictors. Thirdly, the other two methods lose their powers when random error follows Cauchy distribution and sometimes fail to control FDR. In this subsection, we illustrate the practical applications of our procedures by a real data analysis. The group inference is helpful to decide whether a group of predictors are important or not for the response. If we find a group of predictors are important, we would like to know which specific predictors in the group are significant. For this aim, our developed multiple testing procedure is useful. For instance, researchers may aim to test whether a gene pathway, consisting of high-dimensional genes for the same biological function, is important for a certain clinical outcome, given the other high-dimensional genes. When determining that a certain gene pathway is important, researchers need to further identify specific genes within the pathway which are important for a certain clinical outcome.We apply our methods on a dataset about riboflavin (vitamin B2) production rate with Bacillus Subtilis. This dataset is made publicly by [7] and has been analyzed by many authors, for instance [37] , [47] , [29] , and [22] . The dataset riboflavin can be obtained from the R package hdi. It consists of n = 71 observations of strains of Bacillus Subtilis and p = 4088 covariates, measuring the log-expression levels of 4088 genes. The response variable is the logarithm of the riboflavin production rate.Our goal is to detect which genes are associated with riboflavin production rate. Like most existing studies, we first reduce ultrahigh-dimension to a moderate high-dimension. Here we pick out first 300 genes by distance correlation based screening [33] . We first conduct global testing on these 300 genes. The p-value of our group inference procedure is 1.286e -04, indicating that the null hypothesis is rejected and the selected 300 genes are influential for riboflavin production rate. Next, we further use FDR control procedure to select the important genes in these 300 genes.By implementing our proposed FDR control procedure with the FDR level of 0.1, we identify 10 genes that are significantly associated with the response. That is G I = {YTGB_at, YCKE_at, YXLE_at, YXLD_at, YJCJ_at, XHLA_at, xepA_at, YCGO_at, RPLP_at, XKDS_at}. If the FDR level is set as 0.2, 5 more genes will be selected. That is G II = G I ∪{SPOIISA_at, YHCB_at, XKDI_at, YJCF_at, XHLB_at}.We further conduct group inference on the selected subsets G I , G II and their complement sets G c I , G c II . As expected, our group inference procedure finds again that G I , G II are significant while G c I , G c II are not. The corresponding p-values of G I , G c I , G II and G c II are 6.748e-06, 7.257e-01, 9.334e-06 and 9.370e-01, respectively. These results suggest that the genes selected by the FDR control procedure are really influential.We compare these selected genes with other methods. For example, the multi-sample-splitting method proposed in [37] identified YXLD_at; [47] did not select any gene using the de-sparsified Lasso; [29] only selected two genes: YXLD_at and YXLE_at and [22] claimed YCKE_at, XHLA_at, YXLD_at, YDAR_at and YCGN_at as significant. From G I and G II , we can see clearly that the gene YXLD_at is detected by not only [37] , [29] , [22] , but also by our procedure. Besides, the genes YCKE_at, YXLE_at and XHLA_at which are detected by [29] and [22] , are also found by our method. Further, our procedure detects some additional important genes.Assumption A.1. We assume that (i) X j , X ⊤ β h and Z ⊤ j θ j , j = 1, . . . , p are all sub-Gaussian with uniformly bounded sub-Gaussian norms.(ii) c log p n ≤ λ X , λ Y ≤ C log p n for some constants 0 < c ≤ C. (iii) Σ = Cov(X) > 0.The above conditions are mild and widely assumed in the high-dimensional literature.Condition (i) is used to obtain l ∞ norm of 1 √ n n i=1 (X ij -Z ⊤ ij θ j )Z ij in the proof.Sub-Gaussian assumption is a standard condition while working with random predictors in high-dimensional models [49] . However, we should emphasize here that we do not make sub-Gaussian assumption for the error term (or equivalently for the response), which is widely assumed in the literature of high-dimensional inference [40, 17] . Actually we even do not assume the existence of E(ϵ). Assumption A.2. We assume that (i) The function p λ satisfies p λ (0) = 0 and is symmetric around zero (i.e., p λ (t) = p λ (-t), for all t ∈ R).(ii) On the nonnegative real line, the function p λ is nondecreasing.(iii) For t > 0, the function t → p λ (t) t is nonincreasing in t.(iv) The function p λ is differentiable for all t ̸ = 0 and subdifferentiable at t = 0, with lim t→0 + p ′ λ (t) = λL for some L > 0. (v) There exists µ > 0 such that p µ λ (t) := p λ (t) + µ 2 t 2 is convex. Many commonly used penalty functions p λ satisfy all the conditions in Assumption A.2, such as LASSO, SCAD and MCP penalties. Proofs can be found in [35] . By the result in [35] , we have the following proposition for θj . Assumption A.3. The RSC condition in [35] isEQUATIONEQUATIONwhere the α 1 , α 2 are strictly positive constants, the τ 1 , τ 2 are nonnegative constants and the ∆ ∈ R p is a p-dimension column vector.Here we denoteS nj = 1 √ n n i=1 F (Y i ) - 1 2 -Z ⊤ ij γ j X ij -Z ⊤ ij θ j + m j (Y i ) .andS nj = σ -1 j S nj , where m j (y) = E[(X j -Z ⊤ j θ j ){I(Y ≥ y) -F (Y )}].It could be shown that S nj = T nj + o p (1) and Var( S nj ) = 1. Denote S n = ( S nj ) p j=1 and Λ = Var( S n ), whereΛ jk = Cov( S nj , S nk ) is (j, k) component of Λ. And denote S n,G = ( S nj ) j∈G and Λ G = Var( S n,G ) correspondingly. Further denote Ω = β h ∈ R p : ∥β h ∥ 0 ∨ max 1≤j≤p ∥θ j ∥ 0 ≤ s and β hG = (β hj ) j∈G . Assumption A.4. {log(np 0 )} 7 /n = o(1). Assumption A.5. C -1 0 ≤ λ min (Λ) ≤ λ max (Λ) ≤ C 0 for some constant C 0 > 0. Assumption A.6. max 1≤j̸ =k≤p |Λ jk | ≤ c 0 < 1 for some constant c 0 > 0.Assumptions A.4-A.6 are mild and commonly used in the high-dimensional settings. Assumption A.4 is the technical condition to bound the difference between ∥ S n,G ∥ ∞ and ∥N G ∥ ∞ , where N G ∼ N p0 (0, Λ G ). Specially, Assumption A.4 imposes suitable restrictions on the growth rate of p 0 and it is commonly used in the group inference literature, such as [59] and [16] . Assumption A.5 and A.6 are used to establish the limiting distribution of ∥N G ∥ ∞ . Assumptions A.5 requires that the eigenvalues of Λ G are bounded. Assumption A.6 is widely used to establish the limiting distribution of the max type statistic in the high-dimensional settings, see [8] , [54] , and [36] for further examples.The validity of the above bootstrap method requires the following assumptions. Assumption A.9. Suppose that for some b > 0 and q > 0,(j,k)∈A(b) p 2 | Λ jk | 1+ | Λ jk | +q = O p 2 /(log p) 2 .Remark A.2. Recall Λ jk is the covariance between S nj and S nk . Thus, A(b) contains all the strongly correlated pairs and Assumption A.9 requires that the number of these pairs cannot be too large. If |Λ jk | ≤ a 0 for some constant 0 < a 0 ≤ 1, then this assumption holds under |A(b)| = O p 2 1+a 0 -q /(log p) 2 for some b > 0 and q > 0. Similar assumptions were considered in [10] , [20] and [36] .The proof of Theorem 3.1: We denote βh = arg minβ h ∈R p EL n (β h )where L n (β h ) is defined in equation (2.6) . By triangle inequality of l 1 and l 2 norms, it suffices to show that the l 1 and l 2 norms of βhβh and βhβ h satisfy the requirements in this theorem respectively.Firstly we calculate l 1 and l 2 norms of difference between βh and βh . The results are presented in Lemma A.2.Secondly, by Lemma A.8 in Supplementary Material and Theorem 1 in [43] , we can show that βh = (n -1)β h /n, and thus∥ βh -β h ∥ 2 = 1 n ∥β h ∥ 2 , and ∥ βh -β h ∥ 1 = 1 n ∥β h ∥ 1 .By Pr M n,G -2 log p 0 + log log p 0 ≤ t -exp - 1 √ π exp(-t/2) = 0.Note that all the universal constants do not depend on n, p, β h and θ j , j = 1, . . . , p. Thus this theorem is proved.The proof of Theorem A.1: For the convenience of the proof, we introduce some notations. Denote U nj asEQUATIONwhereEQUATIONandh j (Y i , X i ; Y ĩ, X ĩ) = I(Y ĩ ≤ Y i ) - 1 2 -Z ⊤ ij γ j X ij -Z ⊤ ij θ j . (1.26)Denote U nj =: σ -1 j U nj . U nj plays an important role in the proof of this theorem. Let w 1 , . . . , w n be i.i.d. N (0, 1) random variables that are independent of {X i , Y i } n i=1 and denote L ♯ n,G as :L ♯ n,G = max j∈G 2 σ j √ n n i=1 1 n -1 ĩ̸ =i h sym j (Y i , X i ; Y ĩ, X ĩ) -U nj w i 2 ,(1.27)where U nj and h sym j (•) are defined in (1.24) and (1.25) . The proof follows from the general results for the multiplier bootstrap of U -statistic in [12] . Here, we extend these results for approximate high-dimensional U -statistic by results in [13] and rewrite it in a suitable form for our analysis. It is now presented in Lemma A.22 in the Supplementary Material. Suppose ζ 1 is a sequence satisfying ζ 1 1 ∨ log(p 0 /ζ 1 ) → 0 as n → ∞. According to Lemma A.22, to prove this theorem, it suffices to prove ∥ T n,G ∥ ∞ -∥ S n,G ∥ ∞ = o p (ζ 1 ) and Pr w M ♯ n,G -L ♯ n,G ≥ ζ 1 → 0. We denote ζ 1 = C(log(np 0 )) -1/2 ,P M n,G > c (2) G (α) -α = o(1)when H 0,G holds. Note that all the universal constants in the proof do not depend on n, p, β h and θ j , j = 1, . . . , p. We thus have supβ h ∈Ω 0 G sup α∈(0,1) P M n,G > c (2) G (α) -α = o(1).The proof of Theorem 3.3: Using the inequality 2a 1 a 2 ≤ ϵ -1 a 2 1 + ϵa 2 2 for any ϵ > 0, we havemax j∈G σ-1 j √ nβ hj δ j 2 ≤ (1 + ϵ)M n,G + (1 + ϵ -1 ) max j∈G σ-2 j (T nj - √ nβ hj δ j ) 2 . (1.28) Denote η G = | max j∈G σ-2 j (T nj - √ nβ hj δ j ) 2 -max j∈G σ -2 j (S nj - √ nβ hj δ j ) 2 |, and η G = o p ((log p 0 ) -1/2 ) by Lemma A.16 in the Supplementary Material. We can derivelim (n,p0)→∞ inf β h ∈Ω 1 G (2+ϵ0) Pr max j∈G σ-2 j (T nj - √ nβ hj δ j ) 2 -2 log p 0 + log log p 0 ≤ t ≥ lim (n,p0)→∞ inf β h ∈Ω 1 G (2+ϵ0) Pr max j∈G σ -2 j (S nj - √ nβ hj δ j ) 2 -2 log p 0 + log log p 0 ≤ t -η G ≥ lim (n,p0)→∞ inf β h ∈Ω 1 G (2+ϵ0) Pr max j∈G (σ 2 j -β 2 hj δ 2 j ) -1 (S nj - √ nβ hj δ j ) 2 -2 log p 0 + log log p 0 ≤ t -η G = exp - 1 √ π exp(-t/2) .The last equality holds by Lemma A.19, A.20 in the Supplementary Material and Lemma 6 of [8] . Denote t = log log p 0 /2, it implies thatlim (n,p0)→∞ inf β h ∈Ω 1 G (2+ϵ0) Pr max j∈G σ-2 j (T nj - √ nβ hj δ j ) 2 ≤ 2 log p 0 -log log p 0 /2 = 1. (1.29)For any β h ∈ Ω 1 G (2 + ϵ 0 ), applying Lemma A.14 in the Supplementary Material, we derivemax j∈G σ-1 j √ nβ hj δ j 2 = (1 + o(1)) max j∈G σ -1 j √ nβ hj δ j 2 > (2 + ϵ 0 ) log p 0 (1.30)holds with probability 1 as (n, p 0 ) → 1.Substitute (1.29) and (1.30) into (1.28), we derive thatM n,G ≥ 1 1 + ϵ max j∈G σ-1 j √ nβ hj δ j 2 - 1 ϵ max j∈G σ-2 j (T nj - √ nβ hj δ j ) 2 > (2 + ϵ 0 ) log p 0 1 + ϵ - 2 log p 0 ϵ + log log p 0 2ϵ .(i)EQUATIONG (α) ≤ 2 log p 0 by the definition of c (2)EQUATIONG (α) = 1 by (1.29) . By choosing (ϵ 0 , ϵ)satisfying 2+ϵ0 1+ϵ -2 ϵ ≥ 2, lim (n,p0)→∞ inf β h ∈Ω 1 G (2+ϵ0) Pr M n,G > c (i) G (α) = 1, i = 1, 2can be proved.The proof of Theorem A.3: At first we calculate the uniform bound of the difference between constructed statistic T nj and its ideal version S nj . The result is presented in Lemma A.4.For (1.19) , by Lemma 6.1 in [34] , we havemax 1≤j≤p sup 0≤t≤2 √ log p Pr | S nj | ≥ t G(t) -1 ≤ C(log p) -2-γ1 (1.31)for some constant 0 < γ 1 < 1/2. So (1.19) follows from Lemma A.4, and the fact that G(t+ o((log p)-1/2 ))/G(t) = 1 + o(1) uniformly in 0 ≤ t ≤ √ 2 log p.For (1.20) , it suffices to show that and 0 < ν < 1, which will be specified later. We havesup 0≤t≤bp j∈H0 I | S nj | ≥ t q 0 G(t) -1 → 0 in probability. (1.32) Let z 0 < z 1 < . . . < z dp ≤ 1 and t i = G -1 (z i ), where z 0 = G (b p ) , z i = c p /p + cG (t i ) /G (t i+1 ) = 1 + o(1)uniformly in i, and t 0 / 2 log (p/c p ) = 1 + o(1). Note that uniformly for 1 ≤ j ≤ m, G (t i ) /G (t i-1 ) → 1 as p → ∞. The proof of (1.32) reduces to show thatmax 0≤i≤dp j∈H0 I | S nj | ≥ t i q 0 G(t i ) -1 → 0 in probability. (1.33)In fact, for each ϵ > 0, we have Pr max0≤i≤dp j∈H0 I | S nj | ≥ t i -G(t i ) q 0 G(t i ) ≥ ϵ ≤ dp j=0 Pr j∈H0 I | S nj | ≥ t i -G (t i ) q 0 G(t i ) ≥ ϵ/2 . Set I(t) = j∈H 0 {I(| Snj |≥t)-Pr(| Snj |≥t)} q0G(t). By Markov's inequality Pr (|I (t i )| ≥ ϵ/2) ≤ E{I(ti)} 2 ϵ 2 /4 , and it suffices to show dp j=0 E{I(t i )} 2 = o(1). To see this, by (1.31) ,EI 2 (t) = j∈H0 Pr | S nj | ≥ t -Pr 2 | S nj | ≥ t q 2 0 G 2 (t) + j,k∈H0,k̸ =j Pr | S nj | ≥ t, | S nk | ≥ t -Pr | S nj | ≥ t Pr | S nk | ≥ t q 2 0 G 2 (t) ≤ C q 0 G(t) + 1 q 2 0 (j,k)∈A(b):j,k∈H0 Pr | S nj | ≥ t, | S nk | ≥ t G 2 (t) + 1 q 2 0 (j,k)∈A(b) c :j,k∈H0 Pr | S nj | ≥ t, | S nk | ≥ t G 2 (t) -1 = C q 0 G(t) + I 11 (t) + I 12 (t).For (j, k) ∈ A(b) c with j, k ∈ H 0 , applying Lemma 6.1 in [34] , we have I 12 (t) ≤ C(log p) -1-ξ for some ξ > 0 uniformly in 0 < t < √ 2 log p. By Lemma 6.2 in [34] , for (j, k) ∈ A(b) with j, k ∈ H 0 , we havePr | S nj | ≥ t, | S nk | ≥ t ≤ C(t + 1) -2 exp - t 2 1 + |Λ jk | .So thatI 11 (t) ≤ C 1 q 2 0 (j,k)∈A(b):j,k∈H0 (t + 1) -2 exp - t 2 1 + |Λ jk | G -2 (t) ≤ C 1 q 2 0 (j,k)∈A(b):j,k∈H0 G(t) - 2|Λ jk | 1+|Λ jk | .Note that for 0 ≤ t ≤ b p , we have G(t) ≥ G (b p ) = c p /p, so that by assumption A.9 it follows that for some b, q > 0,I 11 (t) ≤ C (j,k)∈A(b):j,k∈H0 p 2|Λ jk | 1+|Λ jk | +q-2 = O 1/(log p) 2 .By the above inequalities, we can prove (1.33) by choosing 0 < ν < 1 so thatdp i=0 E{I (t i )} 2 ≤ C dp i=0 {pG(t i )} -1 + Cd p (log p) -1-ν + (log p) -2 ≤ C dp i=0 1 c p + c 2/3 p e iν + o(1) = o(1).The proof of Theorem A.4: We first consider the case when t, given by (1.18), doesn't exist. In this case, t = √ 2 log p and we consider the eventF 0 = { j∈H0 I(| T nj | ≥ √ 2 log p) ≥ 1}, which mean at least one false positive. In order to show the FDR/FDP can be controlled in this case, we show thatPr(F 0 ) → 0, as (n, p) → ∞. (1.34)And we haveEQUATIONFor any ϵ > 0, we can bound the first term byPr(F 1 ) = Pr j∈H0 I( S nj + T nj -S nj ≥ 2 log p) ≥ 1 ≤ Pr j∈H0 I( S nj ≥ 2 log p -ϵ) ≥ 1 + Pr max j∈H0 | T nj -S nj | ≥ ϵ ≤ p max j∈H0 Pr S nj ≥ 2 log p -ϵ + Pr max j∈H0 | T nj -S nj | ≥ ϵ .By Lemma A.4, we know that Pr max j∈H0 | T nj -S nj | ≥ ϵ → 0. For simplify, we rewriteS nj = n i=1 ξ ij / √ n. Since E(ξ ij ) = 0, Var(ξ ij ) = 1 and {ξ ij } n i=1 is a i.i.d. sequence, by Lemma 6.1 of [34] , we have sup 0≤t≤2√ log p Pr(| Snj |≥t) G(t) -1 ≤ C(log p) -1 . Now let t = √ 2 log p -ϵ, we have Pr S nj ≥ 2 log p -ϵ ≤ G( 2 log p -ϵ) + C G( √ 2 log p -ϵ) log pfor j ∈ H 0 uniformly. Hence p max j∈H0 Pr S nj ≥ √ 2 log p -ϵ ≤ (C + 1)pG( √ 2 log p -ϵ), which goes to zero as (n, p) → ∞. By symmetry, we know that Pr(F 2 ) in (1.35) also goes to 0. Therefore (1.34) is proved. Now consider the case when 0 ≤ t ≤ b p holds. We haveFDP( t) = j∈H0 I(| T nj | ≥ t) max p j=1 I(| T nj | ≥ t), 1 ≤ q 0 G( t) max p j=1 I(| T nj | ≥ t), 1 (1 + A p ),whereA p = sup 0≤t≤bp j∈H 0 I(| Tnj |≥t) q0G(t) -1 . Note that by definition q0G( t) max{ p j=1 I(| Tnj |≥ t),1} ≤ q0αp . The proof is complete if A p → 0 in probability, which has been shown by Theorem A.3.For the convenience of the proof, we introduce some definitions. Denote U nj asU nj = 1 n(n -1) 1≤i̸ = ĩ≤n h j (Y i , X i ; Y ĩ, X ĩ), where h j (Y i , X i ; Y ĩ, X ĩ) = 1 2 I(Y ĩ ≤ Y i ) - 1 2 -Z ⊤ ij γ j X ij -Z ⊤ ij θ j (1.36) + I(Y i ≤ Y ĩ) - 1 2 -Z ⊤ ĩj γ j X ĩj -Z ⊤ ĩj θ j .And denote U nj =: σ -1 j U nj . U nj plays an important role in the proof of main results. Let s X be the sparsity level for parameter θ j as s X = max 1≤j≤p ∥θ j ∥ 0 .Lemma A.1. Under Assumptions A.1 and A.2, the θj defined in (2.5) satisfies∥ θj -θ j ∥ 1 ≤ cλ X s X with probability at least 1 -c 1 exp(-c 2 log p). Here, (c, c 1 , c 2 ) are universal constants. Lemma A.2. Under assumptions in Theorem 3.1, βh defined in (3.11) satisfies ∥ βh -βh ∥ 2 ≤ c 0 λ Y √ s Y , and ∥ βh -βh ∥ 1 ≤ c ′ 0 λ Y s Y with probability at least 1 -c 1 exp(-c 2 log p). Here βh = arg min β h ∈R p EL n (β h ) with L n (β h ) defined in equation (2.6).Proof. The main part of the proof is to verify RSC condition in assumption A.3 and get exponential inequality of ∥∇L n (β h )∥ ∞ . Based on these results, we use Lemma A.5 to prove theorem A.2. Thus, we focus on RSC condition at the beginning. Note that for ∆ ∈ R p ,E n (∆) = {∇L n (β h + ∆) -∇L n (β h )} ⊤ ∆ = 1 n n i=1 (X ⊤ i ∆) 2 .From Proportion 2 in [38] , we then haveEQUATIONwith probability at least 1 -c 1 exp{-c 2 n} for an appropriate choice of α 1 . By the arithmetric mean-geometric mean inequality,τ 1 log p n ∥∆∥ 1 ∥∆∥ 2 ≤ α 1 2 ∥∆∥ 2 2 + τ 2 1 2α 1 log p n ∥∆∥ 2 1 ,and consequently,EQUATION1 , which establishes (1.22) in assumption A.3. As square loss function is convex, condition (1.23) then follows via Lemma 8 in [35] . So we verify RSC condition completely.Next, we construct exponential inequalities of ∇L n (β h ). By conditions of lemma A.5, it suffices to show that there exist universal constants (c, c 1 , c 2 ) such thatEQUATIONBy lemma A.7, we can choose proper (c, c 1 , c 2 ) such that both RSC condition and inequality (1.38) are satisfied. As λ Y is proportional to log p n , the claimed l 1 and l 2 bounds then follow directly from Lemma A.5 by choosing R proportional to 1 λ Y .Lemma A.3. Under the assumptions in Theorem 3.2,lim n→∞ sup β h ∈Ω 0 G sup t∈R Pr max j∈G | T nj | ≤ t -Pr ∥N G ∥ ∞ ≤ t = 0,whereN G ∼ N p0 (0, Λ G ).Proof. By Lemma A.16, we obtain|max j∈G | T nj | -max j∈G | S nj || ≤ max j∈G | T nj -S nj | ≤ Cq(n, p)for some constant C, with probability 1 as (n, p) → ∞. And q(n, p) is a sequence satisfying q(n, p) = o((log p 0 ) -1/2 ). This implies thatPr max j∈G | T nj | ≤ t -Pr ∥N G ∥ ∞ ≤ t ≤Pr max j∈G | S nj | ≤ t + Cq(n, p) -Pr ∥N G ∥ ∞ ≤ t + Cq(n, p) + Pr ∥N G ∥ ∞ ≤ t + Cq(n, p) -Pr ∥N G ∥ ∞ ≤ t =:I 1 + I 2 .By sub-Gaussian condition in Assumption A.1 and Assumption A.8, lim sup n→∞ sup t∈R I 1 ≤ 0. By the Gaussian anti-concentration inequality in Lemma A.20, we derive that sup t∈R I 2 ≤ Cq(n, p) 1 ∨ log(p 0 /q(n, p)) for some constant C. As q(n, p) = o((log p 0 ) -1/2 ), we derive thatlim sup n→∞ sup t∈R Pr max j∈G | T nj | ≤ t -Pr ∥N G ∥ ∞ ≤ t ≤ 0.Similarly,lim inf n→∞ inf t∈R Pr max j∈G | T nj | ≤ t -Pr ∥N G ∥ ∞ ≤ t ≥ 0.This completes the proof.Lemma A.4. Under the assumptions in Theorem A.3,max j∈H0 | T nj -S nj | = o p ( 1 √ log p ).Proof. The proof of this Lemma is similar to Lemma A.16. We only need replace G with H 0 and replace p 0 with p.Lemma A.5. (Theorem 1 in [35] ) Suppose the regularizer ρ λ satisfies Assumption A.2, the empirical loss L n satisfies the RSC condition with 3 4 µ < α 1 , and β * is feasible for the objective. Consider any choice of λ such that4 L • max ∇L n (β * ) ∞ , α 2 log p n ≤ λ ≤ α 2 6RL (1.39)and suppose n ≥EQUATION2 ) α 2 2 log p. Then any vector β satisfying the first-order necessary condition (3.11) satisfies the error boundsEQUATIONwhere k = ∥β * ∥ 0 .Lemma A.6. Consider a U -statisticsU = 1 n(n -1) i̸ =j h(Z i , Z j )with a kernel h based on i.i.d. random variables Z 1 , . . . , Z n . Then for ∀s ∈ R,E exp{s(U -EU )} ≤ E exp s N N i=1 h(Z i , Z i+N ) -EUwhere N = ⌊ n 2 ⌋, which represents the greatest integer less than or equal to n 2 .Proof. It can be seen in the proof of Lemma 14 in [43] .Lemma A.7. Under the Assumptions A.1 and A.2, there exist universal constants (c, c 1 , c 2 ) such thatEQUATIONwhere ∇L n ( βh ) is the gradient of L n ( βh ), which is a p-dimension column vector.Proof. Suppose ∇ k L n ( βh ) is the k-th component of ∇L n ( βh ). Note that ∇ k L n ( βh ) = 1 n n i=1 X ik X ⊤ i βh -{F n (Y i ) - 1 2 } = 1 n n i=1 X ik X ⊤ i βh -{F (Y i ) - 1 2 } -X ik F n (Y i ) -F (Y i ) = 1 n n i=1 X ik X ⊤ i βh -{F (Y i ) - 1 2 } - 1 n 2 n i=1 n ĩ=1 X ik I(Y ĩ ≤ Y i ) -F (Y i ) = 1 n n i=1 X ik X ⊤ i βh -{F (Y i ) - 1 2 } - n -1 n • 1 n(n -1) 1≤i̸ = ĩ≤n X ik I(Y ĩ ≤ Y i ) -F (Y i ) - 1 n 2 n i=1 X ik 1 -F (Y i ) ,where I(•) is an indicator function. So we decompose ∇ k L n ( βh ) into three parts:∇ k L n ( βh ) = 1 n n i=1 X ik X ⊤ i βh -{F (Y i ) - 1 2 } - n -1 n A k - 1 n 2 n i=1 X ik 1 -F (Y i ) ,whereA k is a U -statistics with kernel h k (y 1 , x 1k ; y 2 , x 2k ) = 1 2 x 1k I(y 2 ≤ y 1 ) -F (y 1 ) + x 2k I(y 1 ≤ y 2 ) -F (y 2 ) . Then Pr |∇ k L n ( βh )| ≥ c log p n ≤ Pr |A k | ≥ t (1.42) + Pr 1 n n i=1 X ik X ⊤ i βh -{F (Y i ) - 1 2 } ≥ t (1.43) + Pr 1 n 2 n i=1 X ik 1 -F (Y i ) ≥ t (1.44) with t = c 3 log p n .Firstly, for probability (1.42) . For ∀s ∈ REQUATION≤ exp(-st)E{exp(sA k )}(2)≤ exp(-st)E exp s N N i=1 h k (Y i , X ik ; Y N +i , X (N +i)k ) ≤ exp(-st)E exp s 2N N i=1 (X ik + X (N +i)k ) ≤ exp -st + C 1 K 2 2N s 2 ,where C 1 is an absolute constant, K = max 1≤k≤p ∥X 1k ∥ ψ2 and N = ⌊ n 2 ⌋. Note that (1) holds by Markov inequality, (2) holds by Lemma A.6. Denote s = N t C1K 2 , we havePr A k ≥ t ≤ exp - N t 2 2C 1 K 2 .Similarly, we havePr |A k | ≥ t ≤ 2 exp - N t 2 2C 1 K 2 .(1.45)Secondly, for probability (1.43). As bothX ⊤ i βh -{F (Y i ) -1 2 } n i=1 and {X ik } n i=1 are sub- Gaussian random variables, X ik X ⊤ i βh -{F (Y i ) -1 2 } n i=1 is a zero mean sub-Exponential random variable sequence with sub-Exponential norm X ik X ⊤ i βh -{F (Y i ) -1 2 } ψ1 ≤ KK ′ , where K ′ = ∥X ⊤ 1 βh -{F (Y 1 ) -1 2 }∥ ψ2 .Then by inequality for r.v. with sub-exponential sum (Corollary 2.8.3 in [49] ),EQUATIONwhere C 2 is an absolute constant.Thirdly, for probability (1.44) . Using general Hoeffding's inequality (Theorem 2.6.3 in [49] ) directly,EQUATIONwhere C 3 is an absolute constant. As exp(-n 3 ) = o(exp(-n)), probability (1.44) is trivial compared to other two terms. Combining inequalities (1.45) , (1.46) and (1.47), we have Pr ∥∇L n ( βh )∥ ∞ ≥ c log p n = Pr max 1≤k≤p |∇ k L n ( βh )| ≥ c log p n ≤ p max 1≤k≤p Pr |∇ k L n ( βh )| ≥ c log p n ≤ 4p exp(-C 4 log p), where C 4 = min{ c 2 36C1K 2 , C2c 2 9K 2 K ′2 }.β h = Cov(F (Y ), β ⊤ X) β ⊤ Σβ β.Proof. The proof of this lemma is similar to the proof of theorem 2.1 in [32] . Using Jensen's inequality and LC condition, we can get this result. We omit the proof.Lemma A.9. Suppose p λ X satisfies Assumption A.2. And under the Assumption A.1, we have1 n n i=1 Z ij (X ij -Z ⊤ ij θj ) ∞ ≤ Lλ X ,where θj is defined in (2.5) and L is a constant only depend on p λ X .Proof. By the property of θj , we have1 n n i=1 Z ijk (X ij -Z ⊤ ij θj ) = ∇p λ X ( θjk )with probability 1. Under assumption A.2 and by Lemma 4(a) in [35] , all subgradients and derivatives of p λ X are bounded in magnitude by Lλ X , that is,max 1≤k≤p |∇p λ X ( θjk )| ≤ Lλ X .And this lemma is proved.Lemma A. Z ij {F n (Y i ) -1/2 -Z ⊤ ij γ j } ∞ ≥ c log p n ≤ c 1 exp(-c 2 log p).Proof. At first we decompose1 n n i=1 Z ij {F n (Y i ) -1/2 -Z ⊤ ij γ j } as 1 n n i=1 Z ij {F n (Y i ) -1/2 -Z ⊤ ij γ j } = ∇L n (β h ) + 1 n n i=1 β hj X ij Z ij = ∇L n (β h ) + 1 n n i=1 β hj {X ij Z ij -E(X j Z j )} + β hj E(X j Z j ) =: I 1 + I 2 + I 3 .If ∥β h ∥ ∞ = 0, then I 2 , I 3 = 0. Thus this Lemma is proved by Lemma A.7. So in the following proof, suppose ∥β h ∥ ∞ > 0, that is, there exists 1 ≤ j ≤ p, such that|β hj | > 0. As ∥β h ∥ ∞ ≲ log p/n, there exist a constant c 3 > 0 such that ∥β h ∥ ∞ ≤ c 3 log p/n. Denote c ′ = c/3, c ′ 1 = c 1 /3, we have Pr 1 n n i=1 Z ij {F n (Y i ) -1/2 -Z ⊤ ij γ j } ∞ ≥ c log p n ≤ 3 k=1 Pr ∥I k ∥ ∞ ≥ c ′ log p n . So it suffices to prove Pr(∥I k ∥ ∞ ≥ c ′ log p/n) ≤ c ′ 1 exp(-c 2 log p), k = 1, 2, 3. For k = 1,this inequality holds by lemma A.7 and choosing proper constants.For k = 2, by sub-Gaussian Assumption in Assumption A.1 and applying Bernstein inequality(Theorem 2.8.2 in [49] ), we derive (1.48)Pr ∥I 2 ∥ ∞ ≥ c ′ log p n ≤ p max 1≤k≤p Pr β hj n n i=1 {X ij Z ijk -E(X j Z jk )} ≥ c ′ log p n ≤ 2p max 1≤k≤p exp -c 4 min c ′2 log p β 2 hj K 2 n , C √ log p |β hj |K √ n n ≤ 2p exp -c 4 min c ′ c 3 K 2 , c ′ c 3 K n ,Proof. To prove this lemma, we defineσ2 j = 1 n n i=1 X ij -Z ⊤ ij θ j e ij + m j (Y i ) 2 ,whereEQUATIONj is an unbiased estimator of σ 2 j . It follows by the triangle inequality that max j∈G σ 2 j -σ2 j ≤ max j∈G σ 2 j -σ2 j + max j∈G σ2 j -σ2 j . So it suffices to prove both terms are o p (log p 0 ) -1 . We begin with estimation of the rate of max j∈G σ 2 j -σ2 j . Note thatmax j∈G σ 2 j -σ2 j = max j∈G 1 n n i=1 X ij -Z ⊤ ij θ j e ij + m j (Y i ) 2 -σ 2 j .Since Pr β h (max j∈G σ 2 j -σ2 j ≥ t) ≤ p 0 max j∈G Pr β h ( σ 2 j -σ2 j ≥ t) and for ∀t ≥ σ 2 j ,Pr β h X ij -Z ⊤ ij θ j e ij + m j (Y i ) 2 -σ 2 j ≥ t ≤Pr β h X ij -Z ⊤ ij θ j e ij + m j (Y i ) 2 ≥ t -σ 2 j =Pr β h X ij -Z ⊤ ij θ j e ij + m j (Y i ) ≥ t -σ 2 j .By sub-Gaussian assumption, X ij -Z ⊤ ij θ j e ij + m j (Y i ) is a sub-Exponential random variable with uniform bounded sub-Exponential norm. When t is large enough, by the definition of sub-Exponential random variable, we could find some constants L 1 , L 2 such that Pr β h X ij -Z ⊤ ij θ j e ij + m j (Y i ) ≥ t -σ 2 j ≤ L 1 exp(-L 2 √ t).(X ij -Z ⊤ ij θj )ê ij + mj (Y i ) 2 -(X ij -Z ⊤ ij θ j )e ij + m j (Y i ) 2 ≤ max j∈G 1 n n i=1 (X ij -Z ⊤ ij θj )ê ij + mj (Y i ) -(X ij -Z ⊤ ij θ j )e ij -m j (Y i ) 2 + max j∈G 2 n n i=1 (X ij -Z ⊤ ij θ j )e ij + m j (Y i ) (X ij -Z ⊤ ij θj )ê ij + mj (Y i ) -(X ij -Z ⊤ ij θ j )e ij -m j (Y i ) ≤ max j∈G 2 n n i=1 (X ij -Z ⊤ ij θj )ê ij -(X ij -Z ⊤ ij θ j )e ij 2 + max j∈G 2 n n i=1 mj (Y i ) -m j (Y i ) 2 + max j∈G 2 n n i=1 (X ij -Z ⊤ ij θ j )e ij + m j (Y i ) (X ij -Z ⊤ ij θj )ê ij + mj (Y i ) -(X ij -Z ⊤ ij θ j )e ij -m j (Y i ) =:2(I 1 + I 2 + I 3 ).So it suffices to estimate the converge rate of I i , i = 1, 2, 3 respectively. Decomposing I 1 , it follows thatI 1 ≤ 6 ν=1 I 1ν ,where I 1ν is given byI 11 = max j∈G 1 n n i=1 (X ij -Z ⊤ ij θ j ) 2 (ê ij -e ij ) 2 , I 12 = max j∈G 1 n n i=1 Z ⊤ ij (θ j -θj ) 2 e 2 ij , I 13 = max j∈G 1 n n i=1 Z ⊤ ij (θ j -θj ) 2 (ê ij -e ij ) 2 , I 14 = max j∈G 2 n n i=1 (X ij -Z ⊤ ij θ j )e ij Z ⊤ ij (θ j -θj )(ê ij -e ij ) , I 15 = max j∈G 2 n n i=1 (X ij -Z ⊤ ij θ j )Z ⊤ ij (θ j -θj )(ê ij -e ij ) 2 , I 16 = max j∈G 2 n n i=1 Z ⊤ ij (θ j -θj ) 2 e ij (ê ij -e ij ) .Then we deal with I 1ν , ν = 1, . . . , 6 respectively. For simplify, in the rest of proof, we abbreviate 1 ≤ i ≤ n, j ∈ G as i, j respectively.For I 11 , we have Proof. First we decompose ( T nj -σ-1I 11 ≤ max i,j (X ij -Z ⊤ ij θ j ) 2 max j 1 n n i=1 (ê ij -e ij ) 2 = O p (log(np))O p (n -1 s log p) = O p (n -j √ nβ hj δ j ) -( S nj -σ -1 j √ nβ hj δ j ) as: ( T nj -σ-1 j √ nβ hj δ j ) -( S nj -σ -1 j √ nβ hj δ j ) =σ -1 j (T nj -S nj ) + (σ -1 j -σ -1 j )(S nj - √ nβ hj δ j ) = : A 1j + A 2j .Next we deal with max j∈G |A ij |, i = 1, 2 respectively. We handle max j∈G |A 1j | at first. We decompose A 1j asA 1j = σ-1 j 1 n 3/2 n i=1 ( 1 2 -Z ⊤ ij γ j )(X ij -Z ⊤ ij θ j ) + (θ j -θj ) ⊤ √ n n i=1 Z ij {F n (Y i ) -1/2 -Z ⊤ ij γ j } + (γ j -γj ) ⊤ √ n n i=1 Z ij (X ij -Z ⊤ ij θj ) -n -1/2 U n + ( √ nU nj -S nj ) =: σ-1 j I 1j + I 2j + I 3j + I 4j + I 5jwhere U nj is defined in (1.24) . So it suffices to prove max j∈G |I ij | = O p (r 1 (n, p, p 0 , s)), i = 1, 2, 3, 4, 5 respectively. For a sufficient small constant c ′ > 0 and a sufficient large constant C > 0, we have following results.Step 1, for i = 1. max j∈G |I ij | = O p (r 1 (n, p, p 0 , s)) holds by sub-Gaussian Assumption and Bernstein inequality. The proof is similar to step 2 in the proof of Lemma A.10, we omit it.Step 2, for i = 2. By Proposition A.1 and Lemma A.10, we havePr(|I 2j | ≥ c √ nλ X λ Y s X ) ≤ c 1 exp(-c 2 p),where (c, c 1 , c 2 ) are constants defined in these lemma and independent of the choice of j. We can derive that√ nλ X λ Y s X = O(s log p/ √ n) by the fact that λ X , λ Y ≍ log p/nStep 3, for i = 3. It's similar to step 2 by using Theorem 3.1 and Lemma A.9. We omit it.Step 4, for i = 4. Note that I 4j = n -1/2 (U nj -β hj δ j ) + n -1/2 β hj δ j and n -1/2 max j∈G |β hj δ j | = O(n -1 √ log p). It suffices to deal with n -1/2 (U nj -β hj δ j ). We begin the proof in the same way as we argue in the proof of lemma A.7. Applying lemma A.6, for ∀s ∈ R we haveEQUATIONwhere N = ⌊n/2⌋ and t∈ R + . Denote W ij = h j (Y i , X i ; Y N +i , X N +i ) -β hj δ j , i = 1. . . . , N .By the sub-Gaussian assumption made in assumption A.1, {W ij } N i=1 is an i.i.d. sub-Exponential sequence and max 1≤j≤p ∥W •j ∥ ψ1 < ∞. Thus by the same way used in the proof of Bernstein inequality( [49] ), we derivemax j∈G Pr(|U nj -β hj δ j | ≥ t) ≤ 2 exp -c min t 2 (max 1≤j≤p ∥W •j ∥ ψ1 ) 2 , t max 1≤j≤p ∥W •j ∥ ψ1 n ,where c is a constant. Denote t = Cr 1 (n, p, p 0 , s), we prove it as n ≫ t -2 .Step 5, for i = 5. By the theory of U -statistics, U nj -S nj / √ n is a two-order canonical U -statistic. Thus max j∈G |I 5j | = O p {n -1 (log p 0 ) 3/2 log(np 0 ) + n -1/2 log p 0 } by Lemma A.21. In summary, we derive thatmax j∈G |A 1j | = O p (r 1 (n, p, p 0 , s)) since σj ≥ C min -ψ n ≥ C min /2 > 0 for n large enough. Next we prove max j∈G |A 2j | = o p (r 2 (n, p, p 0 , s)(log p 0 ) 1/2 ). For ∀ c ′ > 0 we have Pr β h max j∈G |A 2j | ≥ 4 3 C ′-3 c ′ Cr 2 (n, p, p 0 , s)(log p 0 ) 1/2 ≤ Pr β h max j∈G σ 2 j -σ2 j ≥ c ′ r 2 (n, p, p 0 , s) + Pr β h max j∈G S nj - √ nβ hj δ j ≥ C(log p 0 ) 1/2 = o(1).The last equality holds by Lemma A.14 and A.15. Note that all the universal constants do not depend on n, p, β h and θ j , j = 1, . . . , p. Thus this lemma is proved.Lemma A.17. Under the Assumptions A.1 and A.2,max j∈G | T nj - √ n U nj | = O p (r(n, p, p 0 , s)).Proof. The proof of this Lemma is similar to Lemma A.16.Lemma A.18. Under the Assumptions in theorem A.1, if s = o √ n/(log p(log np 0 ) 3/2 ) , thenPr w M ♯ n,G -L ♯ n,G ≥ C(log(np 0 )) -1/2 → 0 in probability,where L ♯ n,G is defined in (1.27) in the main text and C is a sufficient large constant.Proof. By definition,M ♯ n,G -L ♯ n,G ≤ max j∈G 2 n ( T nj - √ n U nj ) n i=1 w i + max j∈G 2 √ n(n -1) n i=1 ĩ̸ =i σ-1 j ĥj (Y i , X i ; Y ĩ, X ĩ) -σ -1 j h j (Y i , X i ; Y ĩ, X ĩ) w i =: I 1 + I 2 .Thus it suffices to prove Pr w (I i ≥ C(log p 0 ) -1/2 ) → 0 in probability, i = 1, 2.For I 1 , since w i is standard Gaussian r.v., the Hoeffding inequality implies thatEQUATIONwhere C ′ is a constant. Applying Lemma A.17, Pr w (I 1 ≥ C(log(np 0 )) -1/2 ) = o p (1) is proved.For I 2 , the Hoeffding inequality also implies thatPr w (I 2 ≥ C(log(np 0 )) -1/2 ) ≤ 2 exp - C ′ n(n -1) 2 (log(np 0 )) -1 max j∈G n i=1 ĩ̸ =i σ-1 j ĥj -σ -1 j h j 2 .For simplify, we abbreviate ĥj (Y i , X i ; Y ĩ, X ĩ), h j (Y i , X i ; Y ĩ, X ĩ) as ĥj , h j respectively. In the rest of proof, we abbreviate max 1≤i≤n , max ĩ̸ =i , max j∈G as max i , max ĩ, max j respectively. To prove Pr w (I 2 ≥ C(log(np 0 )) -1/2 ) = o p (1), it suffices to prove that1 n(n -1) 2 max j∈G n i=1 ĩ̸ =i σ-1 j ĥj -σ -1 j h j 2 = o p (log(np 0 )) -1 .By the Hölder inequality,1 n(n -1) 2 max j∈G n i=1 ĩ̸ =i σ-1 j ĥj -σ -1 j h j 2 ≤ max i, ĩ,j σ-1 j ĥj -σ -1 j h j ≤ max i, ĩ,j σ-1 j ĥj -σ-1 j h j + max i, ĩ,j σ-1 j h j -σ -1 j h j = : I 21 + I 22 .Thus it suffices to prove I 2k = o p ((log(np 0 )) -1/2 ), k = 1, 2 respectively. For I 21 , it holds thatI 21 = max i,j |(X ij -Z ⊤ ij θj )η ij -(X ij -Z ⊤ ij θ j )η ij | = max i,j |(X ij -Z ⊤ ij θ j )(η ij -η ij )| + max i,j |Z ⊤ ij ( θj -θ j )(η ij -η ij )| + max i,j |Z ⊤ ij ( θj -θ j )η ij |, where ηij = K -Z ⊤ ij γj , η ij = K -Z ⊤ ij γ j and K = 1/2 or -1/2. Note that ηij -η ij = Z ⊤ ij (γ j -γj ) and max i,j |η ij -η ij | ≤ max i,j |Z ij | • ∥γ j -γj ∥ 1 = O p s log(np 0 ) log p n ,where the last step follows from the s log p/n convergence of ∥γ j -γj ∥ 1 and max i,j |Z ij | = O p ( log(np 0 )), due to the sub-Gaussian properties of Z ij . Following the similar arguments, we havemax i,j |(X ij -Z ⊤ ij θ j )(η ij -η ij )| ≤ max i,j |X ij -Z ⊤ ij θ j | max i,j |η ij -η ij | = O p s (log(np 0 )) 2 log p n , max i,j |Z ⊤ ij ( θj -θ j )(η ij -η ij )| ≤ ∥ θj -θ j ∥ 1 max i,j |Z ij | max i,j |η ij -η ij | = O p s 2 log(np 0 ) log p n , max i,j |Z ⊤ ij ( θj -θ j )η ij | ≤ ∥ θj -θ j ∥ 1 max i,j |Z ij | max i,j |η ij | = O p s (log(np 0 )) 2 log p n .In summary, we haveI 21 = O p s (log(np 0 )) 2 log p n = o p (log(np 0 )) -1/2 ,where the last equality holds by Assumption A.7. For I 22 , it holds thatI 22 ≤ max j |σ -1 j -σ -1 j | max i, ĩ,j |h j | = o p (log(np 0 )) -1/2 ,where the last inequality holds by Lemma A.14 and max i, ĩ,j |h j | = O p (log(p 0 n)). In summary, this lemma is proved.Lemma A.19. Let X 1 , . . . , X n be n i.i.d random vectors, where X i = (X i1 , . . . , X id ) ⊤ . Assume that there are some constants 0< c 1 < C 1 such that E X 2 ij ≥ c 1 and X ij is sub-Exponential with ∥X ij ∥ ψ1 ≤ C 1 for all 1 ≤ j ≤ d. If (log(dn)) 7 /n = o(1), then lim n→∞ sup t∈R P 1 √ n n i=1 X i ∞ ≤ t -Pr ∥N∥ ∞ ≤ t = 0,where N ∼ N d 0, E(XX ⊤ ) .NeurIPS Paper Checklist Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: Type of compute workers and memory are provided in Section 4.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [NA] Justification: This paper is a foundational research and not tied to particular applications.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: This paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: The authors cite the original paper that produced the code package.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. 
 
 Appendix:  


Decision: The focus of the paper is an algorithm for group inference in single-index models which is robust under heavy-tailed noise in high dimension: under some sparsity hypothesis, an equivalent linear estimation objective allows to obtain a robust estimator converging to the true parameters.

The paper received positive feedback from all reviewers because of its clarity and the proposed ideas to deal with the problem of group inference under heavy tails in the response in single-index models. The results are further supported by numerical evidence. During the rebuttal discussion, limitations in the analyzed set-up were addressed and clarified in the resubmitted version, and further numerical evidence was added to the manuscript. I share positive evaluations of the Reviewers and recommend therefore acceptance.

Reviews: review: paper_summary: paper_summary: The paper proposes an algorithm for group inference in single-index models with an unknown link function, that is robust to heavy-tailed noise in the responses. The central idea of the approach is based on the property of elliptical distributions that the linear input-output correlations remain along a fixed direction for any transformation of the labels. The algorithm introduces a linear estimation objective for transformed labels resulting in an estimator converging to the true parameters under sparsity assumptions. Crucially, the obtained estimator yields a robust test for group inference satisfying the orthogonality property, which relaxes the requirements of separation between zero and non-zero coefficients. The obtained test further satisfies honesty, namely the uniform convergence of Type 1 error. The work concludes with numerical tests of the proposed approach.
strength_weakness: strength_weakness: A few limitations of the setup are not explicitly discussed:
1) The robustness is restricted to variability in the noise of the response function, not w.r.t general perturbations to the distribution. The paper also doesn’t provide quantitative bounds on robustness.
2) The condition kappa_h \neq 0 excludes even link functions and in particular the problem of sparse phase retrieval.
3) Theorem 3.1 imposes limitations on the sparsity level that should be explicit. The results in the paper seems to be restricted to constant sparsity levels/

The choice of h=F for robustness appears to be due to F(Y) being uniformly distributed in [0,1] independent of the distribution of Y. This should be more discussed explicitly.
questions: questions: Questions:

Why are λX, λY not directly described in Theorem 3.1?

Why is the transformation H=F particularly suitable for robustness? Is it because of the distribution of F(Y) being agnostic to the distribution of Y?

For sub-exponential variables, Lemma A.11 appears to be weaker than Bernstein’s inequality. 

Could you clarify the difference w.r.t Bernstein’s inequality and related inequalities for random variables with bounded Orlicz norms?
limitations: limitations: Some missing discussion on limitations is described above under “weaknesses”. The work is primarily of a theoretical nature and has no potential negative societal impact.
paper_summary: paper_summary: This paper proposes a robust group hypothesis testing procedure for a high-dimensional single index model based on a data-driven transformation of the response variable. The key observation is that under the the linearity condition (LC) of the predictors, a wide class of single index models can be equivalently (for testing purposes) converted into a linear model with a transformed response variable. The main contribution of this paper is the introduction of the distribution transformation (a data-driven approach) of the response variable, and derive the asymptotic distribution of the resulting test statistic. Numerical performance of the proposed method is evaluated through simulation studies.
strength_weakness: strength_weakness: The simulation study can be improved.
questions: questions: 1. The introduction of the distribution function $F_n$ in (2.7) is interesting but also a little arbitrary. I suspect that any other functions that can be written as linear combinations of functions of $Y_i$'s can be used in (2.7), for example, the kernel density estimator, the kernel smoother of $Y_i$'s. The resulting test statistic should still be asymptotically equivalent to a U-statistic. Could the authors comment on this?

2. Whenever the term "robustness" is involved, it implies loss of efficiency in some cases. In my opinion, it is important to make a transparent evaluation of the limitations of the proposed procedure. In the simulation study, one can add a comparison to the method with $h(Y)=Y$ as described in section 3.3. It will be informative if one can showcase in which cases which procedure is more powerful, and thus understand which method to use in practice.

3. In the simulation study, I would also suggest a graphical comparison between different methods by gradually increasing $p_{out}$ from $0$ to, say, $0.5$ for a more complete picture.
limitations: limitations: The simulation study can be improved.
paper_summary: paper_summary: This paper studies the high-dimensional single index model (SIM), which takes the form $Y=g(X^T\beta, \epsilon)$ with $\epsilon$ and $X$ being orthogonal. Although this model has flexibility and interpretability, its efficiency is adversely affected by outlying observations and heavy-tailed distributions. The paper improves this in the following 3 aspects:

(1) they extend the rank-LASSO procedure to include both convex and non-convex penalties and establish error bound of any local optimum of the empirical objective;

(2) they provide asymptotically honest group inference procedures based on the idea of orthogonalization for testing the joint effect of many predictors;

(3) they develop a multiple testing procedure for determining if the individual coefficients are relevant simultaneously, and show that it is able to control the FDR asymptotically.
strength_weakness: strength_weakness: I think the paper is quite well-written. The only issue may be that the authors can provide more about the background of the problem,
questions: questions: N/A
limitations: limitations: N/A

Comments: title: response to rebuttal
keep the score
comment: I want to thank you for addressing my concerns. They are very helpful. Congratulations for a solid work!
Thanks for the rebuttal. I tend to keep the score.
Dear authors, dear reviewers,

the discussion period has begun as the authors have provided their rebuttals.
I encourage the reviewers to read all the reviews and the corresponding rebuttals: the current period might be an opportunity for further clarification on the paper results and in general to engage in an open and constructive exchange.

Many thanks for your work.
The AC