Scale-invariant Optimal Sampling for Rare-events Data and Sparse Models
Authors: Jing Wang
HaiYing Wang
Hao Zhang
Venue: NeurIPS 2024 poster (2024)
Keywords: adaptive lasso
imbalance data
optimality criterion
subsampling
Primary area: learning_theory

Abstract: Subsampling is effective in tackling computational challenges for massive data with rare events.  Overly aggressive subsampling may adversely affect estimation efficiency, and optimal subsampling is essential to mitigate the information loss. However, existing optimal subsampling probabilities depends on data scales, and some scaling transformations may result in inefficient subsamples. This problem is more significant when there are inactive features, because their influence on the subsampling probabilities can be arbitrarily magnified by inappropriate scaling transformations. We tackle this challenge and introduce a scale-invariant optimal subsampling function in the context of sparse models, where inactive features are commonly assumed. Instead of focusing on estimating model parameters, we define an optimal subsampling function to minimize the prediction error,  using adaptive lasso as an example to outline the estimation procedure and study its theoretical guarantee.  We first introduce the adaptive lasso estimator for rare-events data and establish its oracle properties, thereby validating the use of subsampling. Then we derive a scale-invariant optimal subsampling function that minimizes the prediction error of the inverse probability weighted (IPW) adaptive lasso.  Finally, we present an estimator based on the maximum sampled conditional likelihood (MSCL) to further improve the estimation efficiency. We conduct numerical experiments using both simulated and real-world data sets to demonstrate the performance of the proposed methods.

Scale-invariant Optimal Sampling for Rare-events Data with Sparse Models
Jing Wang
HaiYing Wang Department of Statistics University of Connecticut Storrs, CT 06269
Department of Statistics University of Connecticut Storrs, CT 06269 j jing. wwwgwuconn. edu
haiying- wang@uconn edu
Hao Helen Zhang Department of Mathematics
University of Arizona hzhang@math. arizona. edu
Abstract
Subsampling is effective in tackling computational challenges for massive data with rare events. Overly aggressive subsampling may adversely affect estimation efficiency, and optimal subsampling is essential mitigate the information loss. However, existing optimal subsampling probabilities depends on data scales, and some scaling transformations may result in inefficient subsamples. This problem is more significant when there are inactive features, because their influence on the subsampling probabilities can be arbitrarily magnified by inappropriate scaling transformations. We tackle this challenge and introduce scale-invariant optimal subsampling function in the context of sparse models, where inactive features are commonly assumed. Instead of focusing on estimating model parameters, we define an optimal subsampling function minimize the prediction error, using adaptive lasso example to outline the estimation procedure and study its theoretical guarantee. first introduce the adaptive lasso estimator for rare-events data and establish its oracle properties, thereby validating the use of subsampling. Then we derive a scale-invariant optimal subsampling function that minimizes the prediction eиrоr of the inverse probability weighted (IPW) adaptive lasso. Finally, we present an estimator based on the maximum sampled conditional likelihood (MSCL) to further improve the estimation efficiency. We conduct numerical experiments using both simulated and real-world data sets to demonstrate the performance of the proposed methods.
1 Introduction
Rare-events data refer to binary-response data that are highly imbalanced, i.e., the number of zeros (a.k.a "controls" or "negative instances") are possibly hundreds or thousands times as large as the number of ones (a.k.a. "cases" or "positive instances"). This type of data is common in various fields, such as medicine, natural science, political science, and social science, where examples of rare events can be rare diseases, natural disasters, wars, and financial crises, respectively. Modern technologies also prompt us to pay more attention to rare-events data. For example, in modern online recommendation systems, clicks are usually rare events compared with nonclicks. Statistical analyses, including parameter estimation and inferences, pose unique challenges for rare-events data because of high imbalance. In addition, rare-events data often involve sparse models. For instance,
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

rare diseases might be linked to limited number of key genes. Therefore, researchers frequently adopt sparse models in genome-wide association studies for analyzing rare diseases. A different yet related example is the use of deep neural networks to predict click-through rates in modern online recommendation systems. These networks are typically overparameterízed, necessitating methods that balance rare-events data with the sparsity of the underlying models. Data balancing is a popular approach to overcome challenges caused by imbalanced data and is usually accomplished through subsampling the zeros [5, 15] or oversampling the ones [3. 12, 16, 4]. In addition, rare-events data are often massive in order to obtain an adequate number of ones, and computation is demanding. Therefore, focus the subsampling approach since it addresses the imbalance issue and reduce the computational burden simultaneously.
It is shown in [20] that the efficiency of parameter estimation is essentially determined by the number of ones for rare-events logistic regression, and subsampling does not reduce the estimation efficiency as long as sufficient zeros are kept. In case of excessive removal of zeros, [22] developed an optimal sampling approach to minimize information loss. However, the optimal sampling probabilities in [22] are scale-dependent, which may lead to inefficient results. Figure illustrates the issue using simulated example, with details in Section D.I of the appendix. We generate the data from the same logistic regression model and transform one of the covariates with different scales s - 0.01, 0.1, 1, 10, and 100. Then we apply two optimal subsampling methods in [22], labeled with "A-OS" and "L-OS" in Figure 1. It is observed that the prediction errors of A-OS and L-OS significantly impacted by the data scaling. The A-OS may perform similarly to the Uni (simple random sampling or uniform sampling) in Figure la when dependent issue is not specific to logistic regression and rare-events data in [22]; it is a wide concern in literature for various data types and models, including but not limited to 29, 21, 14, 26, 25, 24). In this paper, we propose a scale-invariant optimal subsampling method to overcome the issue. is labeled "P-OS" in Figure 1.
0.01; so is the L-OS in Figure 1b when
100. This scale-
10-08
A-
10-12.4
prediction mox
L-O5
P-OS
predition
1-4
10044.
10
10
10a
10-
10
10
10
10'
10
10a
(a) Non-sparse parameter
(b) Sparse parameter
Figure 1: Prediction errors with different scale transformation of the same model. (a): with non-sparse parameter (-1, 11, -0.01, -0.01, -0.01, -0.01) (b): with sparse parameter (- 0, 0, 0, 0,0
The scale-dependence issue can seriously impact variable selection results for sparse models, where true parameters are zero for inactive covariates. In this case, inactive variables may be arbitrarily transformed without changing the underlying model, but the A-OS or L-OS would be highly influ- enced and may lead to misleading results. To resolve this issue, we investigate scale-invariant optimal subsampling in the context of variable selection, for which one main goal is to distinguish active and inactive features.
Penalty-based feature selection methods are widely used. Specifically, the adaptive lasso is a popular choice due to its oracle properties, convexity, and practical ease of implementation [see 30, 28]. While penalization methods have been used for bias reduction in rare-events analysis [7], variable selection for rare-events data has not been investigated. Conducting effective variable selection is difficult in the context of rare-events data analysis, mainly due to the scarcity of information available for ones. An inaccurate variable selection result can subsequently impact both the effectiveness of optimal subsampling and the efficiency of parameter estimation. In this paper, we address the challenge of variable selection the context of rare-events data. First, we propose the full data adaptive lasso and study its theoretical properties. Next, we introduce novel subsampling estimator that seamlessly combines penalty-based variable selection and optimal sampling into one unified framework for rare-events data. The implementation of the adaptive lasso requires pilot estimator to construct
2

---

data-dependee deteadeeeedeet weights for covariates. Given that optimal sampling also relies on pilot estimates [see 23, 1]. the adaptive lasso emerges as a natural choice for conducting variable selection method in the context of subsampled rare-events data. We validate the new estimators by proving their oracle properties and also develop an efficient algorithm to facilitate their practical implementation when handling massive real-world data sets. In summary, our main contributions are listed as follows:
We propose scale-invariant optimal subsampling to enhance parameter estimation and variable selection. Existing optimal subsampling methods are scale-dependent, which may lead to unreliable or misleading results.
We define adaptive lasso and establish its oracle properties for rare-events data, which show that the asymptotic variances are determined by the number of ones in the data and the active features in the model.
We present a practical subsampling algorithm based on optimal probabilities that significantly reduces the computational burden and accelerates the optimization for penalty-based feature selection methods.
The rest of the paper is organized as follows. Section 2 introduces the model setup. Section investigates nonuniform sampling and variable selection tailored for rare-events data. We propose new methods construct scale-invariant optimal probabilities. Section 4 discusses theoretical properties of the MSCL. estimator and presents a two-step algorithm to implement the proposed methods. Section 5 conducts numerical experiments simulated and real data sets. Section 6 concludes the paper. Proofs and mathematical details are presented in the appendix.
2 Background and model setup
We use the subscript, to indicate the true parameters. For a p-dimensional vector æ, we use z(j) to represent its i-th element. For an index subset A C {i: 1.2, P?.,,.. we use ((AA) denote the subvector of æ, whose elements correspond to the indexes A. Furthermore, we use 0: to denote a..". use "~-" denote convergence in distribution, " "_"," denote convergence in probability, and use "0.44" to denote convergence almost surely. We use I denote an identity matrix of suitable dimension and use 0 to denote a vector of zeros of a suitable dimension.
Let (Æ1, 1/1, (Iz,, .2.,..., (æN, yN) denote N sample points from the joint distribution of (æ, y). where { denote the p-dimensional predictors and { the binary responses. Assume that the probability of y being one (y 1) given or is
((=+)))))))
c.(2)..)
p(x: 0.) = Pyy = =1/x) = 1 f(x 1 es(=;0,)*
where 0, = (ar, BT is the vector of true parameters and f(æ; 3z) is a smooth function of Br. For rare-events data, N1 < No. where N1 Σ is the number of ones (i.e. yi = 1) and No . N1 is the number of zeros (i.e. y/ 0). Following the model setup used in [22], we assume that a + -00 as N + Dc, which implies that, under appropriate moment conditions,
N1 E E(p(x;0,) No
+ o(1) 1Efortt(man + o(1) + 0, almost surely.
1 - E{p(2; 0.)}
Under this assumption, the asymptotic variance of the full data maximum likelihood estimator (MLE) is of order 1/N, instead of 1/N, indicating that the estimation efficiency is determined by the number of rare ones. Therefore, we can keep all the ones and sample the zeros to save computational costs. There could be variance inflation due to aggressive subsampling, and [22] developed optimal subsampling functions to reduce the variance inflation. Specifically, the authors proposed non-uniform optimal sampling functions under the A- and L-optimality criteria, respectively, as follows: pa (x) x p(æ; M j(æ; 0.) and spf.alle (x) x p(æ; 0,) [[](;; 0,) where M = ((ef( ) (x: } and g(x; 0) denotes the derivative of g( (x;; 0) with respect to 0. However, the sampling functions pacalee (x) and pi (a:) proposed in [22] depend on the scale of æ. and may not perform well for certain measurement scale of For example, if g( 0,) a + x TT, then pf (x) is proportional to + ||x||, which will be influenced by the scale of Similarly, scale
3

---

changes in z may also change - (z), although the impact may not be in the same direction, as demonstrated in Figure 1. Besides parameter estimation, variable selection another important topic. which has not been studied in the literature on rare-events data. This work aims to fill this gap.
3 Nonuniform sampling with variable selection for rare-events data
The adaptive lasso [30, 28] is popular variable selection method because it has oracle properties and is easy to implement. We define the full data adaptive lasso for rare-events data as
1B 1=1 1ppc))
gadp = arg max Vag(a;; 0) log(1 + c0(m=)0))] AN mle
(2)
where Av and are tuning parameters, and 3pI is consistent pilot estimator of Br. In practice, it is common to set y = 1. In the literature, iterative algorithms such as coordinate descent are commonly used solve the adaptive lasso [8]. However, their computational demand can become prohibitive when dealing with massive data. It is feasible to alleviate the computational burden by subsampling zeros and create smaller subset of data for adaptive lasso. To be specific, consider Algorithm 1.
Algorithm Poisson Subsampling algorithm
1: For 1, I:...
2: if y/ then
3: include (2;; Vs) in the subsample;
4: else
5: Compute ((x;)) and generate u4 ~ 00, 1]:
6: if' u4 S (2, 91) then
7: include (æ;, 3/.) and record pạ((;; in the subsample;
8: end if
9:
end if
The inclusion probability in Algorithm for the ith observation is T(.;; vi) = + (1 yı (x:)) where p is the baseline sampling rate for the zeros and (pxx) > satisfies Ef = Let the subsample from Algorithm be {æ=²u,, you Nia. responses. We introduce an Inverse Probability Weighting (IPW) adaptive lasso estimator to correct for the bias, defined
I Ti=1 which is biased since TTx, W1) ''s depend on the
@ndp == arg max
N
y
sub
g( 0) log(1 + gla
AN
(3)
Txx bub
0
To space, put the general assumptions used throughout this paper in Section B.1 of the appendix. We use A to denote the set of indexes of active variables, i.e., A = {j Pii) * 0} and A denote the set of indexes of inactive variables, i.e., A (ú BB() = 03. We first study the asymptotic properties of godp in the following theorem.
Theorem 1. Let Spt be a consistent pilot estimate such that AN/ (v (NN ( " xo for E A Under Assumptions 1-4, if AN /NN- - 0, then the IPW adaptive lasso estimator defined in (3) has the following properties:
1. Consistency in variable selection: The estimated active set A.. {j Audp # 0} satisfies that
lim P(A. A) =
2. Asymptotic normality: The estimator of the active parameter vector satisfies that
VN1Y -1/2 (gadp w(A) w(A) - Oc(AA) ~ NO,
where Vw(A) = E ((f((a)))) M Mw
E Eff((a))) {Mid) + cVoubca). (2)0) .002
M(A) = E (e (z: 0.). VeubCA) Md ((m) [[A (o; 0.) Mūdy c

---

lim oo en /. and (x))(x; 0.) consists of the elements of gradient vector g(x; O.) with in- dexes the active set A.
Remark 1. Theorem I shows that the estimation efficiency of Ouip is predominantly determined by the number of ones instead of the full data size. The term sub((A) is the variation inflation due to subsampling. The full data adaptive lasso in (2) correspond the scenario with p and ((xx) 1, for which с = -im e O. Intuitively, с can interpreted as the imbalance rate in the subsample. If we include sufficient zeros (c 0), the subsampling does not reduce the estimation efficiency of 0
inflated variance if c * O. To minimize the information loss due to sampling, we derive optimal functions as follows, where - (a:) corresponds to the A-optimality criterion [17] and - (0l-00(()) corresponds to the L-optimality criterion [17] in design of experiments. Here, the A-optimality minimizes the trace of the asymptotic variance of Gadp the L-optimality focuses on the asymptotic variance of a linearly
From Theorem 1,
that there maybe information loss reflected
w(AA)'
transformed estimator ( w(AA) which is proportional to Mw(A). The A-optimality criterion has a more direct interpretation, while an advantage of the L-optimality criterion is that the resulting optimal function is often faster to calculate.
Proposition! 1. The A A-optimal function that minimizes tr( (r(W((A) is
pl:: M (z;; 0,)
adp PA-OS (x)
E. p(a; (p)a:00))]lm 0.) M (n 0.)
The L-optimal function that minimizes tr( (r(M)(A)) is
pla: θ 19 () x; 0.)
pr (x) = E (p(x; 0,) ((( 0.)]
(5)
os)=)) relies only the active variables. This implies that a first-step pilot estimator given by the adaptive lasso algorithm can benefit from sparse estimation methods when calculating optimal probabilities. For example, employing the standard lasso can effectively eliminate a large number of inactive variables to facilitate the computation of optimal adp (a)) and adp (z)) However, in practice, pilot estimators are often obtained from small subsample size, introducing additional uncertainty. Therefore, it becomes crucial to exercise caution and be conservative by over-selecting variables during the first step to prevent the exclusion of important variables. As a consequence, although theoretically 4x FA- (a)) and -p..-os (ar) do not depend on inactive variables, they are affected by inactive variables practical implementations.
Unlike the optimal sampling function in [22], - os (x) (or adp
3.1 Scale invariant optimal function
As discussed in Section 1, scaling dependent optimal probabilities may impact the performance of variable selection in practice. To address the issue, we propose to construct a scale invariant optimal function by focusing on the prediction error of an estimator θ, defined below.
[p(=:; 0) p((:; 0.) = / pll: 0) p(:; p(=;((x)))))) -p)x:0) dlP..
MSPE(()
Ex
where P.z the probability measure of x. The probability term p(æ;; Q2) involves both the covariates and the parameter vector 6,, and it often does not depend on the scale of æ. For example, in the logistic regression model, the value p(z; 0.) only related x B. If we change the scale of z())- the value of would change accordingly under the same data-generating model and so p(æ; 0.) remains the same. Thus, re-scaling covariates would not affect this criterion. In the following, we give an optimal function that minimizes the prediction error.
Theorem 2. Under the assumptions of Theorem 1, for the IPW adaptive lasso estimator defined in (3), its prediction error satisfies
Nie 20, MSPE(0 ())))) ~ E-1 (c((m)))) M
Mr Z (6)
w(A)
5

---

(2 { ((A) (x, 0.) The optimal function that minimizes where (()) ~N, 1), and (A)==
=
the asymplotic mean of the prediction error in (6) is given
pla; 0.) (ajiiaa)(a;; 0.)1
adp PP Op-os(2)
E pla; 0,) 001
M (ca 0.)1
We refer this prediction oriented criterion as P-optimality criterion. As we expect, the optimal function in (7) is unaffected by the scale of for class of functions g. The following proposition proves that Adp (p-os(x)) (x) is invariant to rescaling of æ.
Proposition 2. If g 0) satisfies fhr every non-singular matrix A there exists a non-singular matrix B, such that
g(A; Bd 0) = g( 0),
(8)
then, adp os( is invariant to scale changes of x.
Remark The condition (8) is not restrictive and it is quite easy to satisfy. One simple example of g(æ; 0) that satisfies the condition is a linear function g(a; 0) = a + x. B, which corresponds to the logistic regression. The condition is also satisfied by more complex models. For example, consider an L-layer neural network
+ +b*), g( W W2 W....., W ",.bb'..bb))........ -")" T T
where W are the weights and b' are the biases in each layer, 1 1.2,... L. If is rescaled to Ax, we can change W to (AT)-IWI so that the value of g does not change. That is
g(Ax; (A W 2 WL b' b'......''' f" (x ( +b))
= -(x: W W W. b' .....bb)). b²))
Penalized MSCL estimator
The IPW estimator in (3) is not the most efficient estimator, because it assigns smaller weights for more informative data points with larger sampling probabilities. To improve the estimation efficiency, propose the penalized MSCL estimator for variable selection given
N
gddp mscl
arg max
Σ Lu gaxbub, 0) log 1 + cl ta 10111 AN Σ
(9)
where Jaub log (pp(a:;u)); The estimator introduced in [22] is defined as the minimizer of the objective function in (9), excluding the penalization term. this paper, we extend this approach by proposing penalized MSCL estimator ensure model sparsity. We present the oracle properties of the penalized MSCL estimator in the following theorem.
Theorem
Let Ar be a consistent pilot estimate such that Ni - Under Assumptions 1-3 and 5, if ANx /N + 0, the estimator based on MSCL,functin with adaptive lasso penalty defined (9) have the following properties:
do for E A%
1. Consistency in variable selection: The estimated active set Anset == Gadp mecl(s) # 0} satisfies that lim P(Ammce = A)==
2. Asymptotic normality: The estimator of the active parameter vector satisfies that
VNV--1/2
(s)o(AA) (SC))A)) -Ou(A)) N(O, I),
1/2 (gudp
(10)
((m: ) 32 (w:8.) Itep-1 (mreef (=m (=18())
where VissctA.)) E f (m:B)))
AmeaicA) and Amsel(A) =
The penalized MSCL estimator has the same asymptotic variance as the MSCL estimator under the true model, indicating that it is more efficient than the penalized IPW estimator [22]. We prove this by comparing the asymptotic variances and present the result in the following theorem.

---

Theorem 4. If the asymptotic variances for finite, i.e., O (()) (A)) < 0c, then V (isc((A) S Vw(A), where the inequalities hold in the sense of Loewner ordering.
m for w(A) in (2) and mmse(() for msel( in (9). are
Thus, we give a practical two-step algorithm based on the penalized MSCL estimator. Since the optimal sampling functions contain unknown values and the adaptive lasso penalty also requires consistent pilot estimator build weights, it natural combine optimal sampling and the adaptive lasso into one unified framework. We recommend to use the lasso for pilot estimation. One reason is that it does estimation and variable selection simultaneously, and excluding some inactive variables improves the estimation accuracy of optimal probabilities. This also reduces the computational burden for subsequent steps. Another reason is that the lasso estimator tends to include more variables in practice and therefore has a low risk of excluding important variables in the pilot step. We present an outline the practical implementation in Algorithm 2. More details are given in Section C.
Algorithm 2 Two-step subsampling adaptive lasso algorithm
1: First stage screening:
Take a pilot sample of expected sample size Npt using (*((:) = Po + ye(pi -oo)
and obtain a lasso penalized MSCL pilot estimator and an estimated active set Api- Calculate approximate optimal sampling probabilities (#(
yi (1
y: W:)p((;; ] based on (4), (5), or (7).
2: Second stage screening: Use Algorithm 1 with the estimated optimal sampling probabilities
to obtain a subsample of expected sample size Naub and compute the adaptive lasso penalized
MSCL estimator based on Apl-
Numerical experiments
In this section, we use numerical experiments on both simulated and real data to investigate the performances of our proposed optimal subsampling and variable selection procedures.
5.1 Simulation design
We consider a logistic regression with g(:; 0) x and the following three true parameters B of dimension 50. We set different a, so that the proportion of ones is 0.005:
(1) Case A: B. = (0.75, 0.75,,55..... 0 0.75, 0, 0.750, 0.75, 0 and a; = -5.8.
(2) Case B: B3 (3, (3.-2,00 07, 0.85, 0, 0.75, 0% T and a, 6.2.
(3) Case C: G. = (3, 07. 0.85, 0%, and a; -7.5.
Here, O₂ denotes the zero vector of dimension d. We use PA and PA to denote the number of active and inactive variables, respectively, and assume that æ a normal random vector. The active components T(AJ))
j I P. of have variances 0.25 and the inactive components ((As))) S j S of have variances 100/p3 100 (p./AP 1)3 ;... 100/33, 3 100/23, 100/13. The correlation between the i-th and j-th elements of is 0.5li-5,, i,j S p. We repeat our experiments 5 = 500 times generating N = 500000 data points in each run and use pilot sample of size Npt = 500 for obtaining pilot estimates based on the lasso. We consider uniform sampling, the full data lasso, and the full data adaptive lasso for comparison. We use the 5-fold cross validation and Bayesian information criterion (BIC) to determine the tuning parameter for the lasso and the adaptive lasso, and choose 7 1 for the adaptive lasso.
cross-validation
5.1.1 Estimation and prediction efficiency
We present the empirical median squared error (eMSE) for parameter estimation in Figure 2. All opti- mal sampling estimators outperform the uniform sampling. As the sampling rate increases, sampling estimators outperform the full data lasso estimator eventually in all of the three cases. Among the three optimal subsampling methods, 3pdp performs better than the other two subsampling methods.
7

---

CA
-30
Iogle
logickSE)
-10
4.0
48
-45
sampling raze
compiing rate
samplin mnnping r200
(a) Case
(b) Case
(c) Case C
Figure 2: eMSE for different true parameters with different sampling rates.
1300
-125
-135
LOS
-110
-32.5
-14.0
logloNSSHE)
Alasse (140
-150
-140
233
-145
-34.5
2,005
3,000
3.012
(a) Case A
(b) Case B
(c) Case C
Figure eMPSE of estimated probability with different sampling rates.
Figure shows the results of the empirical median squared prediction error (cMSPE). Similarly to the results eMSE, optimal sampling estimators perform better than the uniform sampling. meaning optimal sampling results in less information loss. It is possible that sampling estimators outperform the full data lasso estimator as the sampling rate increases, despite that the latter uses all of the data. In general, - performs the best among the three optimal subsampling algorithms.
5.1.2 Variable selection and computational complexity
In this section, we discuss the results of variable selection in terms of the first stage screening and the second stage screening. Table presents the mean numbers of selected variables Case C. where the numbers in the parentheses are the corresponding standard errors. Results for Cases A and B are similar so are put in Table the appendix.
Table 1: Mean number of selected variables in Case C
P
first-stage
Uni
A-OS
L-OS
P-OS
0.0025
13.27(0.34)
2.84(0.02)
2.97(0.02)
2.96(0.02)
2.96(0.02)
0.005
12.46(0.32)
2.94(..02)
3.04(0.03)
3.05(0.03)
3.06(0.03)
0.0075
12.76(0.33)
2.97(0.01)
3.04(0.02)
3.03(0.02)
3.03(0.02)
0.01
12.81(0.34)
2.96(0.01)
3.03(0.02)
3.02(0.01)
3.02(0.01)
While the first stage screening significantly reduces the dimension in Table it indeed includes inactive variables expected. In the second stage screeming, the mean numbers of selected variables close to the number of selected variables from uniform sampling is smaller than the true number of active variables especially when the sampling rate is low. This indicates that the second- stage screening of uniform sampling may exclude active variable. We present the rates of missing active variables in Table 2 for Case C. shows that uniform sampling has higher rates excluding active variables than optimal subsampling procedures, so optimal sampling may be preferable in practice. Results for Cases A and are similar and are put Section E.I. We also investgate the rates of selecting the true model in that section.
numbers of active variables for all subsampling methods. However, the mean

---

Table 2: Rates of excluding active variables (false negative rate) in Case C
p
Uni
A-OS
L-OS
P-OS
0.0025
0.168(0.017)
0.086(0.013)
0.088(0.013)
0.084(0.013)
0.005
0.100(0.013)
0.068(0.011)
0.066(0.011)
0.066(0.011)
0.0075
0.066(0.011)
0.046(0.009)
0.048(0.010)
0.046(0.009)
0.01
0.068(0.011)
0.052(0.010)
0.054(0.010)
0.054(0.010)
5.1.3
Computational time
We present the mean computational times of different algorithms in Table 3. Our codes written in the julia programming language [2] and implemented on Linux workstation. The lasso pathes solved with Lasso,jl [13]. As shown in Table 3, subsampling algorithms significantly reduce the computational times compared with full data estimators. Although optimal sampling requires to calculate sampling probabilities, they use only about 0.77% of the computational time that the full data adaptive lasso requires. As we discussed in Section C.2, optimal sampling algorithms reduce both sample size and the data dimension. Therefore, the computational cost of the coordinate decent algorithm, which often requires large number of iterations, is significantly reduced.
Table
Mean computational time (seconds)
Case
Uni
A-OS
L-OS
P-OS
A-lasso (full)
Lasso (full)
A
0.29
1.09
0.91
1.06
129.62
112.97
B
0.31
1.23
1.20
1.27
129.89
122.40
C
0.31
1.02
0.93
1.00
130.33
121.29
5.2
Real data
We evaluate the performances of proposed estimators on two real data sets.
Covtype data set: It is available Attps://archinn wuw.uuissi....................ii edu/n//datasstttetteeeee/// covertype, with N = 581012 observations and 54 covariates - 10 being quantitative and 44 being qualitative with dummy coding. We drop the 14th and 54th columns avoid exact colinearity of the dummy variables. Our goal is to classify whether the forest cover type is Cottonwood/Wiloon (labeled 1) not (labeled as 0). The proportion of Cottonwood/WWila is 0.473%, which is highly imbalanced.
(ii) Font data set: It is available at https: staps://archive ics. ics.uci. www.eeeeeeaaaaaeeeeeeeeee CharacterFFonts Images, with 0.50% of the N 832670 responses being the GADUGI font. The first 10 covariates are about the value, size, and style of the characters and there are additional 400 pixel values of the 20 x 20 images. We remove the 4th, 9th, and 10th covariates because they are constants.
For both data sets, we apply Algorithm 2 on the Oggaiihmic-transformee data. We use pilot samples of size Nyp 1000 for the covtype data and Np = 1500 for the font data due to its higher dimension. Since we do not know the true parameter for real data, use area under the curve (AUC) to measure the performances of subsampling algorithms. We repeat the experiment for S = 500 and compute the empirical median AUC using the full data. The results are summarized in Figure 4. As shown in Figure 4, nonuniform sampling outperforms uniform sampling in general. There is one case for font data set that 1 is worse than the uniform sampling when the sampling rate is high. For the covtype data set, among the three estimators based on optimal sampling, A performs the best and - is worst. For the font data set, - Os and grdp os are similar, and pp-os based on the scale invariant optimal sampling function is significantly better.
6 Conclusion and limitations
In this paper, we investigated the problem of scale-invariant optimal subsampling the context of variable selection for rare-events data. We derived optimal probabilities based on the A- and L-

---

0.830
A-OS L-OS
0.9964
P-OS
Uni
AUC
0.9963
AUC
0.9962
0.826
0.996
0.824
0.0050 0.0075 0.0100 0,0125 0.0150 sampling rale
0,0050 0.0075 0.0100 0.0125 00130 sampling rate
(a) AUCs of the covtype data
(b) AUCs of the font data
Figure 4: Empirical median AUCs for two real data sets
optimality criteria, and discussed their limitations. Furthermore, proposed scale-invariant optimal probabilities based on prediction errors to overcome the limitations. Both analytical and numerical results show the desirable properties of the proposed methods.
Our investigation has the following limitations.
Our proposed criterion optimizes the probabilities by minimizing the asymptotic mean squared error in estimating rare-event probabilities. While this prioritizes the accuracy of estimation, puts less emphasis on the quality of variable selection. Further research is needed to devise optimal probabilities that focus variable selection performance metrics.
Our theoretical analysis based on asymptotic properties, with optimal probabilities defined through the asymptotic normality. Although our results may hold for sufficiently sparse models, they may not generalize to cases where the model dense or over-parameterized. because asymptotic normality may longer be applicable. Therefore, an important direction for future research is to study the non-asymptotic properties of our estimators, such as prediction error bounds. Non-asymptotic behaviors particularly of interest in high- dimensional regimes.
We employ Lasso as the pilot estimator. However, other variable selection methodologies, such sure independence screening, can also be considered. Exploring the impact of different pilot estimators on our method's 's performance represents another avenue for future investigations.
We assume that the underlying full model correctly specified and possesses sparse structure. Our analysis does not account for model misspecification. Further research is required to address scenarios where the model possibly misspecified or where the number of features vastly exceeds the number of observations.
Acknowledgments and Disclosure of Funding
The authors are grateful to Professor Kun Chen for the insightful comments and suggestions on the development of the manuscript. Funding in direct support of this work: NEI grant R21EY035710 10, NSF grant 2105571, UConn CLAS Research Funding Academic Themes, GPUs donated by NVIDIA.
References
[1] M. J. J. Yu, H. Zhang, and H. Wang. Optimal subsampling algorithms for big data regressions. Statistica Sinicu, (1)277799777, 2021.
[2] Bezanson, A. Edelman, S. Karpinski, and B. Shah. Julia: A fresh approach to numerical computing. SIAM review, 59(1):65-98, 2017.
10

---

[3] N. V. Chawla, K. W. Bowyer, L. O. Hall, and W. P. Kegelmeyer, Smote: synthetic minority over-sampling technique. Journal of artificial intelligence research, 16:321-357, 2002.
[4] G. Douzas and F. Bacao. Self-organizing map oversampling (somo) for imbalanced data set learning. Expert systems with Applications, 82:40-52, 2017.
[5] C. Drummond, R. C. Holte, et al. C4. 5, class imbalance, and cost sensitivity: why under- sampling beats over-sampling, In Workshop on learning imbalanced datasets II, volume 11, 2003.
[6] Fan and J. Lv. Sure independence screening for ultrahigh dimensional feature space. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70(5):849- 911, 2008.
[7] D. Firth. Bias reduction of maximum likelihood estimates. Biometrika, 80(1):27-38, 03 1993.
[8] J. Friedman, T. Hastie, H. Höfling, and R. Tibshirani. Pathwise coordinate optimization. The Annals of Applied Statistics, 1(2):302 332, 2007.
[9] J. H. Friedman, T. Hastie, and R. Tibshirani. Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software, 33(1):1-22, 2010.
[10] W. Fu and K. Knight. Asymptotics for lasso-type estimators. The Annals of statistics, 28(5):1356- 1378, 2000.
[11] C. J. Geyer. On the asymptotics constrained m-estimation. The Annals of statistics, pages 1993-2010, 1994.
[12] H. Han, W.-Y. Wang, and B.-H. Mao. Borderline-smote: new over-sampling method in imbalanced data sets learning. In International conference intelligent computing, pages 878-887. Springer, 2005.
[13] JuliaStats. Lasso.jl.
2022.
[14] N. Keret and Gorfine. Analyzing big ehr data -ottimal cox regression subsampling proce- dure with rare events. Journal of the American Statistical Association, 118(544):2262-25) 22752225.. 2023.
[15] X.-Y. Liu, J. Wu, and Z.-H. Zhou. Exploratory undersampling for class-imbalance learning. IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cyhermetics), 39(2):539-555) 2008.
[16] Mathew, C. K. Pang, M. Luo, and W. H. Leong. Classification of imbalanced data by oversampling in kernel space of support vector machines. IEEE transactions on neural networks and learning systems, 4969):4065 4076, 2017.
[17] F. Pukelsheim. Optimal design of experiments. SIAM, 2006.
[18] Shao. Mathematical Statistics, 2nd. Springer- Verlag, New York, 2003.
[19] G. Tripathi. A matrix extension of the cauchy-schwarz inequality. Economics Letters, 63(1):1-3, 1999.
[20] H. Wang. Logistic regression for massive data with rare events. In International Conference on Machine Learning, pages 9829-9836. PMLR, 2020.
[21] H. Wang and Y. Ma. Optimal subsampling for quantile regression in big data. Biometrika, 108(1):99-112, 2021.
[22] H. Wang, A. Zhang, and Wang. Nonuniform negative sampling and log odds correction with rare events data. Advances in Neural Information Processing Systems, 34, 2021.
[23] H. Wang, Zhu, and P. Ma. Optimal subsampling large sample logistic regression, Journal of the American Statistical Association, 113(522):829 844, 2018.

---

[24] J. Wang, J. Zou, and H. Wang. Sampling with replacement vs poisson sampling: comparative study in optimal subsampling. IEEE Transactions on Information Theory, 68(10):6605-660) 6630, 2022.
[25] Y. Yao and H. Wang. Optimal subsampling for softmax regression. Statistical Papers, pages 585-599, 12 2018.
[26] Yu. H. Wang. M. Ai, and Zhang. Optimal distributed subsampling for maximum quasi- likelihood estimators with massive data. Journal of the American Statistical Association, 0(0):1-12, 2020. D0::10. 00000001621499.220.1777777
[27] G.-X. Yuan, C.-H. Ho, and C.-J. Lin. An improved glmnet for 11-regularized logistic regression. J. Mach. Learn. Res., 13:1999-2034 2012.
[28] H. H. Zhang and W. Lu. Adaptive Lasso for Cox's proportional hazards model. Biometrika, 94(3):6117703, 05 2007.
[29] T. Zhang, Y. Ning, and D. Ruppert. Optimal sampling for generalized linear models under measurement constraints. Journal of Compulational and Graphical Statistics, 30(1): 101:))006:1110 14, 2021.
[30] H. Zou. The adaptive lasso and its oracle properties. Journal of the American Statistical Association, 101(476):1418-1429 2006.
A Appendix supplemental material
In this appendix, we present the details of the proof, the practical algorithm and simuation settings in the paper. Details of mathematical proofs are provided in Section B. Details of the practical algorithm are provided in Section C. We present the details of simulation settings in Section and in Section E, we give some additional simulation results.
B Details of mathematical proofs
In this section, we provide details of mathematical proofs.
B.1 General assumptions in the main paper
We begin with some general assumptions used throughout this paper. Assumption 1.
The first, second and third derivatives of f f(x;; 0) and (a(9)) f(x; 3) with respect to are bouned by square intergrable random variable (xx).
Assumption 2. The matrix (x; 0)} is finite and positive definite.
Assumption 3. The subsampling rate satisfies that eN = e e0'pp - e, where 0 S c < 0o is constant.
Assumption 4. The integral E [[p(()) + u (x=) B² (x)] is finite, where B(x)) square- integrable function that dominates the first, second, and third derivatives of f(x; 0) and e f(x; (3) with respect to B.
Assumption 5. The integral E (p is finite.
These assumptions are the same assumptions used in [22]. Here, we remind some notations used in the main paper:
ca+((=;;)
p(x:0) =
1 1 ++
6(x; 0) = p(z: 0) 0)
M MBBfellassae 0.)
12

---

and
( g (a:; 3()
msel
1 1+ op o (x))/(a:3))
To ease the presentation in the following sections, we denote
K 1+ cef(e;
Mw
= E
e Ma. 9 ((1) (2; ((xx)
and aN Ve in the appendix. Note that
}
N1 Σ
y NE
{1 [oo
[1+aat((;()))
1=1
Ne E ellla)))
( (1+ axB (fllss)) (p())
B.2 Proof of Theorem 1
Proof of Theurem I. We consider target of IPW adaptive lasso estimator:
8,
(w) π (x, yi)
(agla 0) logg1 1 e0(2) + AN
Σ ااار
j=1
= -0) + AN Σ
v0 1
j=l
where ún = 1/1 $ Then, we have that üN = a, (6w 0.) is the minimizer of
N (u) Qw (0. + a u) Qw.
Asymptotic normality:
We prove the asymptotic normality part in this paragraph. By Taylor's
expansion.
6z
7w (u)
uT (m)) +
(b (l))::0)(()) 0c) (u)((; 0,)3 Aw + R
ON
20²N
T (x yi)
AN wjaN (BBw
"()) 11(1sss
AN
ON
We first consider the limit behavior the IPW target function by prove the asymptotic normality. In [22], the authors established that under Assumptions to 3,
aN in (0.) ~ M22
N
8
0a2;; (2) 0,) P, M,
(z, 1h)
and
Aw = oP((1,
Rw op())
Thus,
- Lv (0.) ~ -u -uT
Mu.
2
Next, we consider limit behavior the adaptive lasso penalty. Since we assume 3n to be consistent estimator, we know that when j E A, i.e., Bi(r) * O.
ij 1/p1G)) P. (((G)~ > 0.
and
aN (
(110)
"(
- ||(6s) + sgn(Bec))ucp
13

---

Therefore, for j E A. we have that
AN
WjaN (BcS)) +
"Gi) oN
- -1Bcco) oP(1)
aN
since AN = AN/ Ne - O. On the other hand, when j E A,, i.e, Bug) = O, we have that for u()) # 0,
An
WjaN
Be(s) +
u())
AN dr |()) = aN (3p()))
An
MmGil P
aN
ON
aN
since AN/ Aw//VVer Ne (oogs11) _n co. Then, we have that n (u) Tw (u), where
Tw (w) = (A) M Ma) - utAy M 2/ ((A) if f = 0, E A
0
otherwise.
Note that the unique minimizer of N (1)) is (M 1 M ZTA
(A)) 0)T if we put all the indexes of active variables in front. Thus, following the results of [11] and [10], we have the minimizer N (u). iec. v. satisfies that
(A)
ÜN(A) MújM.. / (AA) and (N(A9) O.
Thus,
ÚN(A) = aN aN (( O((a)) - N(O. M MM
Since
Vii ANDI1/2 fereveeeeii fellai (p)))
applying Slusky's theorem, have
V ((A) (0 O((A)) ~ N)..
Consistency in variable selection We prove the consistency in variable selection in this paragraph. From the result of asymptotic normality, we know that Aw(J) - Beo) for every and therefore F(j E Aw) + 1. Thus, we only consider j' E As. When jr E Aw. we know that by K-K-T optimality conditions, have
(gn())))3)))) = (m
which means
X (0w))
aN
ON aN (ww(0w) - (o.)
(0.) + aN
I1 Iz.
a2y
We have known that 41 = (Or)//a Zw. We now prove that proof that 1/ Op(1). We apply Taylor expansion the k-th element of (0w) and have that
an ((s) - (())
1
5,
0,) 0.) o0 (1)9())( (x;;
aN
T Bi
i=1
where,
ŭN aN
(0w -0,) =Op(1),
8,
Â Áw(k)
(y) p(x,; 0:)3 Σ
J(aj (x.;; ())
N
(ii 1/i)
14

---

and
Rw((k) =
N
8;
(2 0k) {1 - 2p(2;; 6.)) D100) (an; 6s o )
20
IT
(Z; Ub)
2
N
6;
6 (121; On)
{ {
iuN
To(k) (o,; 6s)
} ( 02)
2
(2;;; Vi)
00
1
8,
o10, O1) (2:;; 0z) { ( ( (o dw)a)) }
2
7 (2 Vii
N
8;
O (as) Oz)
{ui - p(zr;; aN T (zi, wi)
ÜN.
20
002
where 6₂ is between mlc and 0.. First, we prove that R(n is (r) We have that
[Rwcks] 5
| 11
N
6,
0(0;; 0w) |1 - 2p(x;; 0.)) ( (zr;; 6b) |
203
T (Ti y1i
2 jiv1 1
N
8;
Di(k) (c:;; 0iik) 0n)
| 6w)
+
p 0x)
2a N
x (xx; vi)
00
| 1|
6;
o(I;; 6x)
(ea0) ( (2
de))
| 6n)
2a Nv
I (z M1)
| 112
6,
0 gi)) (cx;; 0x)
p(z,; 6%)
2 N
π (x Vi)
002
j 11?
N
6,
0 g(N) (z);; 6w)
2 N
(2;, V/)
Ui
802
1=1
| 112
85
|| N
62
p pla;; Cx) 6x) +
v Baxs)
2
π (Fi yii
20
1=1 x
(x yi)
|
1P -a
co
5;
(o C(z;: + | 1?
Σ
N
8,
y B(z))
20
T y1)
2axx
(2,, y1)
i=1
| 12 -c a
N
(iorib))C( 0x) + || 1|2
N
5i
Σ
y B(B()))
2NaN
I(x. Vi)
20²
(x
1=1
(1),
where
|am) 6)| {
K| | || + -|(aco::)) 00|}}
C(a;;
+ || (8j () (xi; 0x)
|l00e: Gu)11 +
82 BicRs (20 0z)
00
002
Therefore, we proved that Rw((k) (1) Next, we prove that Aw(a) = (1)... We know that Eim 8e/T(2; v.) (v+- p(=,; 0.) 0(2;; 0.)] O. We also have that for the every element Σm 81/ T (2: yi) {u p( 0,)) (ax;; O.), we have
E aN 2 N of N
22
N
8,
a
{yr- p(a,; 0.)} Begty (zo.; 0.)
n Di yi)
i=1
15

---

Σ
i p(a;; 0c 0:)0 (o:; 0.) - (0 13(a;; 0.)]( - 0.
aN
Thus, due to Chebyshev's inequality, we know that Áw or(1). Since we know that Ei 0:/=T(1 ()) 0,))8(((;; 0+) Hw = (1().. Hence, we have that
(fw = Or(1).
aN
Note that we also have
AN
P 00.
ON
aN 18
Therefore,
P(j' E Aw) S P & ANtr = lm (0m)
Ann (B6)) (0w) 0. =P
AN
aN
Thus, we prove the part of consistency of variable selection.
B.3 Proof of Proposition
We first give a lemma for general optimal functions.
Lemma 1. Assume that h(x)² and (px) are integrable function with E[p(x)} = 1. The optimal
function y²* (z) that minimize the value E (²)))) is given as (²*(()) B[h(a)]'
h(m)
sp(x)
Proof. Appying Cauchy-Schwand inequality, we have that
2
h(x:)
(x)) <<xx)
(()(())
B{h(=)}² E
/o(x))
E
B(p(#))
(x())
(((x)
Therefore, we have that (m))))
2 E(h(x))² and the equality holds and only if /p(xx) Kh(=)/ / p(x)) where K is constant. Therefore, (²"((x) Khh()), and since E(p** = we know that (**(()) = h(z)/
p(w)
Now, we prove Proposition 1.
Proof. We first calculate the optimal function that minimizes tr( (r(VV(A)) We have that
tr(V tr(Ve(()) = { M Mu(a) (A)
=ir - Mij {1 |++
caf()a)))
e (mmt)) 983 (2; 0.) Mid].
p(x)
We focus on the values that related to x(x)) We know that
B ( (1,4)
.22 (a;; 0.))
-
20s (1 + or(1)}E (x- (x:))((x;:00) (2)0 (x:: 0.)
Therefore, we need minimize
(x; O, M (1) 82
0. Mz) (4)
tr
E
(xx)
82 2² (x: O. M 9 (A) (x: O M
(x; 0,) ( ( (x:; 12
E
tr
E
(2))
(0(x))
16

---

Appying Lemma 1. We know that the minimizer is given as
plæ; 0.) M)) (x; 0:)1
(a-)s( E f 0.) || (a) 0.l1
Next, we calculate the optimal function that minimize ur M w(a)). We have that
{ B { col(mna)))
tr (r))) -
e ellan (z: 0.)
(2)
fpxx)
Therefore, we need to minimize
p² xx; Og (A) (x;; 0.)
}
(2; (1)
())((A)( (z; 02) (pxx)
E
EE
((x)
Appying Lemma We know that the minimizer is given as
p(æ: (1) 0.)
E 0 (() 0,)l
B.4 Proof of Theorem 2
Proof. In the proof of Thereom 1, we know that
an (O w(A) - Oa(A)) ~y M M 1/2
w(A) ((A)
To simplify the representation, We define function
h(0) e -2a MSPE(0) 20.E [p(x;;) 0) p(a:;
We have that
[p(x; 0) p(æ: 0
2 (p(x;; 0) - p(x;00) 00(:
00
and
8 {p(x; 0) -p(:; 0.)1
05(x; 0)
200 (2;; 0)
(1: 0) + (p(x; 0) p(:; 0.)}
((:; 0)
0000T
00
+ 2 (p(x; 0) p(=; 0.)) () ()ū(2;
Note that
|p(:; 0) p(a; 0.)l S 2 and thus,
0 (p(æ; 0) p( 0,)3
< p(=;; 0) 0.) 0))(;; 0) s 4B(x),
00
and
8 {p(x; 0) p( 0
T0a0T
(x: 0) (x: () 0) p(æ; 0.) 0ø(2; 0)
$ 2 00 (x:
ōt j( 0)
+ 2 2pp((
p(:; lo(z; 0)ā(2; 0) < 10
Hence, due to donimating convergence theorem, we know that the expectation and derivitive are exchangable. Thus, we have that
Oh(0)
-2a; E (p(x; 0) p(:; 0,) d(x; 0))(;; 0)]
00
17

---

and
Oh(O) 0000T = = 20 v; 0 0) 02 (x; 0) + 2 (2; 0)
00(x: 8)
p(:; 0,))
5(x; 0)
00
+ {p( 0) p(:: 0.) (::; ()ā(2; 0)
Due to donimating convergence theorem, we also know that the first and second derivitive of (())
are continous. We have that
ah(0)
(h((t)
0,
2e-22
2e
(0; Og ) (x; 0.)
a0
0=0.
0000T
0=0.
Note that 20g o²2 ( 0 (x; 0.)]
(x)
+ S2 due to donimating covergence theorem. Now applying
Theorem 1.12(ii) [18], we have that
(p(z;0m(a)) -ppz: p(2; ppl 21
N [ m(A) h(Or(a))) =aze 201 20uB E
2 2 M M M (A)
1/2 w(A) Z(n)
A/ M
((A) M 1/2 w(d) (())
Considering Nı az { (a))))) (1 (1()) applying Slutsky's 's theorem, have that
Nic 20 {ple: ((())) p(2:
- 0
E {ellosah) M1/3 w(A)
M M Z(A)
((A)-
Since Zw(a) = N(O, 1), we have
E
Z Z( M 1/2 M (LA) M (o) ((A) ((A) = AM
w(1)
M 1/2
= [ Mit
{1
11+
cef(o)))))
e (iees 11 (a:; 02)
(fxx)
We focus the values that related to (p(x)) We know that
E B ( (AA)
022 (x; 0.)
2ac {1 or(1)}B (20
Therefore, we need minimize
(xx; O, 0,)9 (x: 0,)
tr
M M Raa
x(xx)
12; 0z) (x; 0t)
E
tr
M M()) M
(xx)
(x; 0, 0, 1/2
82 (x: O M 1/2
(1)
E
tr
x(x)
(x: 0, 0.) 1/2 M J( (z; 1
(A)
(xx)
Appying Lemma 1. We know that the minimizer is given
p(:; 0.) 1/2 (ij 0.)
(p-
E
B(( 1222 Mi 0x))1)
18

---

B.5 Proof of Proposition
Proof. First, know that g( 0) (Ax; B) Since the equation holds for æ and if take derivitive with respect to 0 on both sides, the equation still holds. Thus, we have that
g(x; 0) Bi(Ax; B O)
If we scale the whole covariate variable ã Ax, need to reparameterize O, to 0. O to remain the problem invariant. We have that for ã and O,
g(x; أ( Aa: B² B ((;;
Now, we know that
(- M (e(e 0.) =B=BB---
B- 222
and
S =(E((E((EAA) (a (x; == "DBB")))
Thus, we have that
|| "( - 0.) 12.1 ((;;
=
-g"(x 0.) (B") BB BBBBBBBBBBBBBBBBB BT M BB j(:;; 0.)
(x;
j(: 0,) M
Therefore, the leveraging term invariant. For the probability term, we know that is only related to value g(;; 0,) ( 0,), we know that it does not change after scaling inactive variables. This complete the proof.
B.6 Proof of Theorem 3
Proof of Theorem 3, Consider maximum sampled conditional likeihood function with adaptive lasso penalty:
P [ygg 0) logg1 + (0,,0)))411 + AN Σ أاارا
Ppi Q mscl (0)
- mscl
Oyt
+ Ax
where w 1/l 12, S S Then, have that ún @ mse! 0,) is the minimizer of
8 Umsell (w) = aue
mseel (O, + ax ax²u) u) Q a (0.).
Asymptotic normality:
We prove the asymptotic normality part in this paragraph. By Taylor's
expansion,
N
Opt Nimsel (2) = (0.) +
u T on
Σ
G, or (xx; 0.) (,)(,""(x,; 0,)) A 6, ROp
msee
20
meel +
muell
aN
i=1
An wjaN (1800) +
First, we consider the limit behavior of the MSCL function. In [22], the authors proved that under Assumptions and 3,
O musel (0.) ~n (Apl Alee)) 1/2 Z.
N
[ O 02 ( 0.) P, mscl
19

---

and
- or((),
R op(1).
Thus,
Pul (0.) ~ -u (A msc! 1/2 z m
1 =
A pl
mse! u + 0p(1).
2
Next, we consider the limit behavior of the adaptive lasso penalty. Since we assume 3pl be consistent estimator, know that when j E A. i.e., P((1) # 0,
وj = 18ppco)
P, 18(c))1 0,
and
aN (166) +
"G aN
- |1865s) + sgn(Sec))uug
Therefore, for j E A, we have that
An
(Bus)
"(})
-18tss) P, 0,
WaaN
aN
aN
since AN/aN= = AN/VNNer + O. On the other hand, when j E i.., BBG) = O, we have that for (())= # 0,
An
wjaN (BCS)
u(s)
1Becsll)
AN tiy luunnn
aN 1Bpuc)) |uG)l P, oo,
An
aN
aN
since AN/ AN/(VNeee a (p(r))7)
r3 co. Then, we have that Tissel (u) Ymsol (u), where
6,
} (1) Apl
meel ((A) "[A)(
(A msel)
1/2 ((A) if ( = = Vj A
Tmscl
00
otherwise.
Note that the unique minimizer of (u) is ((A ) (
(A)) 0)* if we put all the indexes of active variables in front. Thus, following the results of [11] and [10], we have the minimizer of Dn
T msee (u). üN satisfies that
@N(A) ~ (Apl 1/2 ((A) and (NN(A)) ~ 0. musel)
Thus,
(N(A)) =a (0 mmsc((A)) -c((a)) ~ N(O, (Apl moac)
We know that
V 11/2 (A ) ) aN 1/2 (( } E-1/2 ( { {1+ aN fi (1)))
VNVV-1/2
Hence, applying Slusky's theorem, we have
1/2 NiV VNV-11/2
ms) ( (museic)) - -e( N(O,1).
Therefore, prove the part of aymptotic normality.
Consistency in variable selection We prove the consistency in variable selection in this paragraph. From the result of asymptotic normality, we know that Blssckg) P, Buy) for every j A and therefore Pij E Amscl) + 1. Thus, we only consider j' E AS. When j' E Amsol, we know that by K-K-T optimality conditions, we have
ppp (mmsc)),
which means
A ue (gngggnn(3m (Bmsses)))))))
p iusel
(mmsct)
aN
NN
20

---

Aa msscl (0.)
aN
{ [in (mscll mscl) "usel ( (0:) 11 + 12.
(
aN
ai
We have known that 1 =
a ~n maze We now prove that proof that 17 = Op()). We
msc1) apply Taylor expansion to the ke-th element of (0 and have that
(0 mscc)
- Legj [A)
(0.(
aN
1
pu 0pt (a,; 01)01))))) 0,)g" (x; (a;;; + 5(0) +
Qer Ra.1
(k))
Nx
ON
i=1
where,
üN aN (6mu msel
0.) =Or(1;
A
p0 (k)
1
ppt {us-
p (ar;; {1) Σ ((ss) 0.)
a'
and
N
(2 6x) {1 - 2pm 0u (x;; 0b) (os) on) ( GN)
mscc((k)
2 N
1=1
2
N
0pl 0pl ØT
(x,; 6x)
T
(z:n) 0x)
} { bw)
20
00
1=1
N
1
ppt 0p1
(xs: (1)) (20;) 6x) { (N 6x juN)
20
N
0 (a)) 6n)
0p1
{
pa 6p (x;; 61) uN
{ii
uN.
20
002
i=1
where 6% is between @mscl and 0,. First, we prove that Rmmak() is or(1). We have that
R
FRBSCcRs) S
||
B Onl Or
(ar;; 0k) 1 12pppt (xs;; 0s)
(0 61)
|o(o:: bh) ||
20
+
2 j 11?
N
0, pr
0ul 0pl
(ar;; 0b)
0gin) (w)
| bb)
2 N
00
| 212
pp pn or
(xr 6x)
Jacks (c 6x))
( [ 02)l
y
| 2|2
Σ
0 Px
ppi 0pl
(zx; 0x)
θ2 g(N (000 (2 Ok)
+
|
112
N
0, yi
g2 ((n))(r: bx)
2
ag2
2a
d02
|| 1
N
|uNll 112
N
20
Σ
6/
pa
(I,; (()s;
20
Σ
0, Yi B(x:)
1=1
|| || 1/
aa
N
C(ar;; 6k) + | 112
N
5.e y; B (ax))
20²
2a
1=1
1=1
|| (e ar
N
0p1 / (1 C(,,) 6x))
||
Σ
N
0 yı B(a:)
NNa
2n
|| (edv
N
- (xx)e (lo + | 11
N
00 x
Σ
0,
y Bxx))
2NaNP
20
1=1
21

---

1uN|12E0 ai
N
[vvl 11 N
(x.
Σ ==1
y ((.))
2NaNP
2a²
1=1
or(1);
where
C(
( (2,; 0)
|( 1/)) +
+ Dp(x) (z;; 0x)
| B.s1|
8 (() 0k)
00
O0
Therefore, proved that R (m (p()). Next, we prove that 2 mssc(k)) Op We know that
1= 0, Uv - p (2,; 0.) g( 0.)
0-1 = O. We also have that for the every element
a
of
2
1= {Us p( 0.)) (( 0.), we have
N 1y: pa - (o (z.;i))) 0.) (60) (x,; 0.) 0pl
V
N
1=1
Σ
pl Gut ;
ps (a,; 0 0.00 (o 0.) 621 = (me
(g(;; 0.)11 0.
NN
ai
Thus, due to Chebyshev's inequality, we know that Á 6pt
mueel op(1). Since we know that
G (a O 0,)g (0;; 0.) = O1().. Hence, we have that EN 0h
( = Op(1).
ON
Note that we also have
Awty AN aN
P oc.
aN
Therefore,
F(J' E AN) S ANisy (gn(iisgnnnna (gn(3mscclc)) = misel (0
ANt0y sgn(s(n((mms)
19 ( (0msel)
}
0.
aN
aN
Thus, prove the part of consistency of variable selection.
B.7
Proof of Theorem 4
Proof. Letting h -( -e (el (( (1x) 0), f hēv, and g h 21b20 we have that
)) =R(fg" E(wn") = E. e ( (x; 0.) M(a),
E(ff") = B(huvT) = E
(m: (x; (1)
MM(AA):
((x))
and
(a,
22 (x: 0.)
(1)
E(gg²) = E(h -1o
Amscl(A).
+ op-
(x)e
Now, applying the matrix form of Sauchy-Sehwwrtz''s inequality (see [19]). we have that
A Amsc((A) = EgggT) Elgg"T " "
M, (M M(A) E ((mm))) { (Vm(A)}-1.
Therefore, simple algebra shows that
VIscK()) = (m)a))))) (A (Amsc((A)) 5 Vo(A))
which complete the proof
22

---

C Details about the practical algorithm and its complexity
C.1 Two-step algorithm
We take a pilot sample by uniform sampling with the sampling rate P Npt/2N1 for the ones and
Po / for the zeros. Denote pilot sample of actual sample size Ni as {(=)' "."")))))") ))
the pilot estimate of f as ptt, and the pilot estimate of the active set as Ap = (i PPios # 0}. We propose the following moment estimators of (A) and (A))
Ni (x
"T 3p) ((A,)) (x/p 0p1)
MP (Api)
(11)
Npl
Po + y (p1 Po)
Noi 2a t B 9 (x:7 0pp)) (Ant)
221 (Aai)
(12)
Np Po u (p1 po)
respectively. We also use the following moment estimator to estimate the denominator of (4):
W A-OS
(13)
Npl
+ up (pl -
Po)
Po
where wA- 0A-08 = p(of" ; MM
pl (AmI) (wl (a, Opi)ll. using (5) or (7), we use
w ww-os = p(=:'' p1) |1 (ل)) (2.²1 : B2)1l.. or
wp -p-os p(a:f'; 0l)) (Ani) ) (Am))) ) ((A)) (aef") Opr)lls
respectively, instead of w in (13). Now, present the proposed two-step procedure in Algorithm 3 with more details than Algorithm in Section 4.
Algorithm 3 Subsampling adaptive lasso algorithm
1 No Take a pilot sample {(=)' yy )= of expected sample size Npt using {T(यE)
Po ys(p1 Po)) and obtain pilot estimator
Na Diu g(x""' 0) - log(1 + (a(ap))0))))
p Ap1 j=1
0 == arg max
Σ
(14)
where No is the actual pilot sample size and I No)(NN/Ni). We call this first stage screeming.
Calculate approximate optimal sampling probabilities (+ v1) = ys + (1 ())p((=)) by replacing (fx;)) with aadp ( 0p)) sudp ( 0pi), or Pp- os (z;; d,), based on (4), (5), or (7), respectively. The denominator of (4) is estimated using (13), and replace w,--ss with -/ -.--s p-os for the denomi- nator (5) or (7), respectively. If using adp A or - estimate M(A) and S(,A4) using the moment estimators in (11) and (12), respectively.
2: Use Algorithm 1 with the estimated optimal sampling probabilities to obtain subsample ((=fub, yubb]
) [=
Naub and compute the adaptive lasso estimator:
No sub
1801
gudp arg max
Σ
Σ
(x O) logg(1 cola (0(07",,)+1111 - AN
Σ
1 1B
JEAp.
where Nob, is the actual subsample size, based on the smaller model obtain from the first stage screening. We call this step the second stage screening.
23

---

Remark 3. Our algorithm naturally integrates the MSCL fimction with the adaptive lasso penalty. can also be implemented when p > N as long as the dimension of selected variables is smaller than N in the first-stage screening. If the model is sparse and the data are massive, this is usually possible in practice. Screening algorithms such as sure independence screening [6] can also be used for the first stage screening to guarantee that the dimension of second-stage screening is smaller than the subsample size. Furthermore, the first stage screening can help to speed up the computation, as shown by the analysis of computational complexity in the next section.
We consider a coordinate desent method to calculate the estimators defined in Algorithm 3. (see [8], [9] and [27]). In each cycle, we need to find an optimal direction at a starting point 0. We consider the quardratic approximation of Q ( + d) Quei meeel (0), which is
ppt (0 d)
Og (0)
P ا5ارر + ddy
811 (xg (x;; 0+ + d) + logg1 cdle, AN Σ
j=1
8,1 yı 0) logs1 +0
+1+ H + AN Σ وو أ
Put mscl (0) d + d p msel
1
d
(0)d + AN Σ
(051i00) + duحl - + 1enl
2
j=1
where
Gal macl
0)
y1
Sot Pa
(xo, 0 0( 0) and
1 meel (0)
())
Øx (æ;, 00 (x:; 0)- Thus, using coordinate desent to obtain the optimal direc- tion, the quadratic approximation for the j-th element is given as
N
ou 00 0g
ppi dai (d + ze;) - meel (d) (mscc))))) = + (mman(0)d]
(6] ()) 2 " (0)=²
2 + 1
+ Av (B0 + d()) + +| + ANuu + dzyي-
Then, we have the value of a that minimize Que (d + ze,) - QQQN (d) is
AN m)at(t)
(0)
amal (0)d) +hm;
if Bin + de) + 2
(5nel(s))) (0)
El () (0)+
ul
(0)
if B1) + d() + 00
mmeeliss)) (0)
B() d
otherwise,
which is the same as
21, - BBn deps } max { -2, + (co) } + B,0 + do),
-" max
where
/0m 21 m))
(0) + mssel
ب (8)a + x100
(is) (0)
and
p
m)) (0) + iusel (0)d) (1)
>x00
22
pOt (js)( (0)
For the special form of g(x; 0) a + f(xT (3,, we know that
N
pJ
(0) = Σ
6, 0x (xi, 01g 82 (x;; 0) &.
0,1 Om
mseel
i=1
24

---

where
fax a (a
G
and 4- - diaglo; 0x (x:;; 0)}. Thus, we have
( (0)d) (5)
= -G" e; = (Gd)" ((e)) (G
Therefore, we can store Gd and keep updating Gd with
G + ze,) = Gd + zGe; = Gd + G
Thus, we do not need to obtain the full matrix Emsel (0) = G We only need to caleu- late the diagnoal elements: (msel(s) (0) GTty G) 1, .... p 1 and { { (0)
(5)
(G)) = 1, ....PP + 1. From the analysis above,
notice that the computaional complexity of one cycle calculating optimal direction d is Op), where Çin denotes the number of inner iteration.
C.2 Computational complexity
We analyze the computational complexity of the two-step algorithm. To facilitate the presentation, we consider a special case for our model when g(x; 0) a + f(x 3), and assume that the number of variables selected at the first-stage screening q. Coordinate descent a widely used optimization algorithm for sol solving lasso and adaptive lasso [see 9]. We consider the improved coordinate descent algorithm proposed in [27], which requires inner iterations to determine an optimal direction and outer iterations to update the estimator. Considering the form of g 0) + f(x" (a:" ()), the computational complexity for coordinate descent with data of size NN and dimension p is O(inn per inner- cycle where Sin represents the number of inner iterations (detailed derivations of this complexity presented Section C). Thus, the computational complexity of full data lasso is (Soutxiia Np), where Soutxin Sout Çin and Sout is the number of outer iterations. The computational complexity of the full data adaptive lasso is O(pple Np² + Soutxin with the MLE as the pilot estimator, and OSShi OSp xinloottiinii Np + Soutxin .N with lasso as the pilot estimator, where Spi and Spl.out xin the iteration numbers in the two pilot estimators, respectively. The coordinate descent algorithm often requires large Soutxin or Spl,outxin while Newton's 's algorithm requires small omle so it is often the case that Gille < Soutxin Therefore, the time complexity of the adaptive lasso is O(CottxinNNp) which is the same as the full data lasso estimator.
clas
Now, we analyze the time complexity of Algorithm 3. We start with the computational complexity of the optimal probabilities, for which the main computational cost to approx- imate || (z;; 0)ll or ||52 || 1/2 M (1)) Il. respectively, for i 1...,N. Since ((A,)) (z,; 0) (1, f B)xT ((Ap))) (T, the computational complexity of calculating (x;; (z;;) 0)'s
O(Ng), and the computational complexity of Mel or Ω2 (Ap) (Api) O(Naig)). Taking the inverse (Ma)) and finding the square root ($2pl 1/2 both take O(q²) time. Thus, the computational complexity of calculating II(NTp) (Api) ) ((Am)) (z,; 0) 's or [K (A.i) )13(( (A,)) ) (z;; Is ONNg+NNpu' + + N O(Ng²). Therefore, the complexity of approximating the optimal probabilities in (4) or (7) is O(Nq²). The computational complexity of approximating the optimal probabilities in (5) is only O(Ng), because there is no
is O O(Npi( + 1)²} =
(Ant))
need compute M 52
(Api) (An))' Next, analyze the complexity of parameter estimation. The
average subsample size with the sampling rate p is on average
E(N1) p(N E(N1)) = NE[s(; 0.) (1' + + o(1)) = ((N +p}.
25

---

Using the coordinate descent algorithm, the computational complexity of the two-step algorithm is pll.ouuxiii NNq² + Sout in + p)g} using optimal probabilities in (4) or (7), and is Sptl Npp + Ng + Xn (e - p) using the optimal probabilities in (5). For optimal probabilities in (4) or (7), when Soutxin > q/(e" + p), the dominating term of the complexity is CoutxinN((e"" + p)g. Remember that Soutxin = Sint<in and Sout is usually large for the coordinate descent algorithm. Therefore, Sout xin > 4/(ee + p) is often satisfied in practice. The dominating term for the time complexity of optimal probabilities in (5) is also Coutxin (C + p)q. Compared with full data estimators, both the sample size and the dimension are reduced. If we set the subsample size to be the same order N,, which is often the case practice for balancing the ones and zeros, the time complexity of Algorithm 2 is of order (Soutxiin Ne q), which is significantly faster than that the full data estimator.
D Detalis of simulation settings
In this section, we present more details of the simulation settings in the main paper. In Section D.I, we provide detailed simulation settings of the example Section 1. In Section D.2, we present detailed settings in Section 5.
D.1 Simulation in Section 1
We first present the detailed settings of the simulations in Section 1, where we illustrate the scale- dependent issues of optimal subsampling probabilities. Our simulation based on logistic regression models with the true parameter Br to be 6-dimentional vectors and covariates æ ~ lognormal E) with the (i, j-ttthhhhee element of E is given as Σg = 0.5 0.5 -1, < i.j < 6. We consider two cases of parameters:
(a) Non-sparse parameter: B, = (1,,-11,-00 00011 -0001 0.01, -0.01-0.001. 0.01) and ot -4.
(b) Sparse parameter: Br (,,,0,0,0,0,000 and a =
We generate full data of size N 500000 according to the above logistic models. To investigate the effects of scale transformation, we multiply the ((6 with (s = 0.01, 0, 1, 100) and divide 3e(6) with the same to remain x =T33, 3. to be the same and thus the logistics regression model does not change. We obtain subsamples with optimal subsampling probabilities described in [22] with transformed æ under each and calculate the resultant subsampling estimators. We set the nominal pilot sample size Npi 800 and nominal subsample size Naub 1000 (see details in [22]). We repeat the experiment for 500 times under each scale and compute the mean prediction error.
D.2 Simulations in Section 5
For the estimation procedures in the second step of our two-step algorithm, we choose 7 1. which means that the weights in the adaptive lasso penalty are ش 1/1/ S j S q, with being the number of selected variables in the first stage screening. Furthermore, we consider uniform sampling, the full data lasso and the full data adaptive lasso baselines for comparison. For the uniform sampling method, we use a similar two-step algorithm as presented in Algorithm but set the sampling function in the second step as spxx) 1, which means the sampling probabilities a constant p. We use lasso to implement the first stage screening and adaptive lasso with y = 1 to implement the second stage screening for fair comparison. For full data lasso, directly apply the lasso algorithm to the full data. For the full data adaptive lasso, we use the full data MLE estimator as the pilot estimator construct the weights and then apply the adaptive lasso algorithm to the full data set.
E
Additional simuations
In this Section, we give some addtional simulation results. More simulation results of variable selection is provided in Section E. We also provide addtional simulation results to compare our approach with standardization to resolve scale dependent issues in Section E.2.
26

---

E.1 Addtional variable selection results
Table 4: Mean number of selected variables in Case and Case B
Case A (five active variables)
P
first-stage
Uni
A-OS
L-OS
P-OS
0.0025
14.87(0.30)
4.99(..02)
5.02(0.02)
5.03(0.02)
5.01(0.02)
0.005
14.68(0.28)
5.05(0.02)
5.07(0.02)
5.06(0.02)
5.08(0.02)
0.0075
14.20(0.27)
5.09(0.03)
5.08(0.02)
5.08(0.02)
5.09(0.03)
0.01
14.46(0.30)
5.13(0.02)
5.13(0.02)
5.13(0.02)
5.12(0.03)
Case B (four active variables)
P
first-stage
Uni
A-OS
L-OS
P-OS
0.0025
16.81(0.33)
3.91(0.02)
4.01(0.03)
4.00(0.02)
4.01(0.03)
0.005
17.88(0.34)
4.04(0.03)
4.10(0.03)
4.08(0.04)
4.07(0.03)
0.0075
17.19(0.33)
4.06(0.02)
4.13(0.03)
4.10(0.03)
4.12(0.03)
0.01
17.51(0.34)
4.07(0.03)
4.08(0.03)
4.08(0.03)
4.08(0.03)
Table 5: Rates excluding active variables (false negative rate) in Case and Case B
Case A
P
Uni
A-OS
L-OS
P-OS
0.0025
0.094(0.011)
0.076(0.012)
0.068(0.011)
0.070(0.012)
0.005
0.048(0.010) 0.052(0.010)
0.050(0.010)
0.044(0.010)
0.0075
0.054(0.010)
0.052(0.010)
0.052(0.010)
0.052(0.010)
0.01
0.040(0.009)
0.036(0.008) 0.036(0.008)
0.036(0.008)
Case B
P
Uni
A-OS
L-OS
P-OS
0.0025
0.108(0.014)
0.074(0.012)
0.074(0.012)
0.070(0.012)
0.005
0.034(0.008) 0.062(0.011)
0.046(0.009)
0.034(0.008)
0.036(0.008)
0.0075
0.046(0.009)
0.046(0.009)
0.044(0.009)
0.01
0.056(0.010)
0.050(0.010)
0.052(0.010)
0.050(0.010)
Table Rates of selecting the true model
Case A
P
Uni
A-OS
L-OS
P-OS
0.0025
0.856(0.016)
0.870(0.015)
0.868(0.015)
0.878(0.015)
0.005
0.872(0.015) 0.884(0.014)
0.880(0.015)
0.872(0.015)
0.0075
0.858(0.016)
0.868(0.015)
0.868(0.015)
0.868(0.015)
0.01
0,854(0.016)
0.866(0.015)
0.862(0.015)
0.868(0.015)
Case B
P
Uni
A-OS
L-OS
P-OS
0.0025
0.862(0.015)
0.872(0.015)
0.864(0.015)
0.876(0.015)
0.005
0.898(0.014)
0.888(0.014)
0.898(0.014)
0.902(0.013)
0.0075
0.848(0.016)
0.850(0.016)
0.854(0.016)
0.848(0.016)
0.01
0.874(0.015)
0.874(0.015)
0.876(0.015)
0.872(0.015)
Case C
P
Uni
A-OS
L-OS
P-OS
0.0025
0.824(0.017)
0.880(0.015) 0.874(0.015)
0.884(0.014)
0.005
0.880(0.015)
0.884(0.014)
0.884(0.014)
0.882(0.014)
0.0075
0.908(0.013)
0.916(0.012)
0.910(0.013)
0.914(0.013)
0.01
0.908(0.013)
0.910(0.013)
0.908(0.013)
0.910(0.013)
is seen in Table that, no subsampling method dominates others. Table and Table 5 shows that uniform sampling has higher rates of excluding active variables than optimal subsampling procedures.
27

---

Although uniform sampling may have a higher rate of selecting the true model in some cases, given that it more likely exclude important variables, optimal sampling may be preferable in practice.
E.2
Comparison with standardization
Another approach to avoid scale-dependency is to standardize the data. We compare the proposed scale-independent optimal probabilities with the approach of data standardization here. For the data standardization approach, we standardize the data, calcualte the optimal probabilities, and then implement subsampled adpative lasso algorithm. We used the same pilot estimation methods for fair comparisons.
We first compare the eMSE and eMSPE in Figure 5 and Figure 6, respectively. We use sP-OS to denote the approach with data standardization and use P-OS to denote the approach without data standardization.
2.8
2.9
aP-O5
3.0
33.6
-3.0
(og(o4SE)
log(oMSE)
lograMSE)
3..8
-32-
-32
3.3-
3.4
-3.4
4,2
-35
0.004
0.005
0,008
0.010
0.004
0.006
0.008
0.004
0.008
0.000
sampling rate
sampling rate
Sampling rate
Case A
(b) Case
(c) Case C
Figure 5: Empirical median squred eтror of estimated probability for different parameters with different sampling rates. The same pilot sample size Npt 500.
In Figure 5, we notice that the performances of P-OS and sP-OS are similar, this is also true for eMSPE. However, standardization may decrease the rate of selecting the true model. We present results of variable selection in Table and Table
-13.8
I9:05
-13.3
-13.9
-13.4
-14.0
(ogicMSPE)
logteMASEE)
-13.5
(o0/eMSPE)
-12
-34.1
-142
-334
137
13.8
-14.4
0.004
1.008
B.008
0.010
9.904
0.006
2.008
6.100
0.004
0.006
0.008
0.010
sampling rate
sampling
sampling rate
(a) Case A
(b) Case
(c) Case C
Figure 6: Empirical median squred error of estimated probability for different parameters with different sampling rates. The same pilot sample size Npt 500.
We notice in Table 7 and Table that the rates of selecting true models by Apdp is higher than (gdd fer-os without much increase on the rates of excluding active variables. Therefore, although standardization an approach to solve the scale-dependency issues, it may decrease the rates of selecting true models in practice.
28

---

Table Rates of selecting true models
A
B
C
sUni
P-OS
sP-OS
sUni
P-OS
sP-OS
sUni
P-OS
sP-OS
0.0025
0.848
0.878
0.862
0.850
0.876
0.862
0.828
0.884
0.880
0.005
0.878
0.872
0.850
0.890
0.902
0.890
0.882
0.882
0.880
0.0075
0.850
0.868
0.846
0.850
0.848
0.842
0.906
0.914
0.910
0.01
0.840
0.868
0.844
0.868
0.872
0.862
0.900
0.910
0.904
Table 8: Rates of excluding active variables (false negtive rate)
C
sUni
P-OS
sP-OS
sUni
P-OS
sP-OS
sUni
P-OS
sP-OS
0.0025
0.088
0.070
0.068
0.106
0.070
0.070
0.164
0.084
0.084
0.005
0.056
0.044
0.044
0.046
0.036
0.036
0.100
0.066
0.066
0.0075
0.052
0.052
0.052
0.062
0.044
0.044
0.062
0.046
0.046
0.01
0.040
0.036
0.036
0.056
0.050
0.050
0.068
0.054
0.054
29

---

NeurIPS Paper Checklist
[Yes] [No]
[NA]
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results
be expected to generalize to other settings. It is fine to include aspirational goals as motivation long it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
In Section 6.
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 encouraged to create separate "Limitations" section in their paper.
The paper should point any strong assumptions and how robust the results 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 few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably provide closed captions for online lectures because it fails 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 faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, 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 impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
30

---

Answer: [Yes]
Justification: Detailed proofs and required assumptions are provided in the appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: All details for numerical experiments provided in Section 5 and in the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, descríbing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is new model (e-g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path reproducing verifying the results. 5. Open access to data and code
31

---

Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: Codes are submitted supplement for anonymity. They will be released in public github repository after the review period.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.cm
conolicyy for more details.
While encourage the release of code and data, understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https: / /ttpss.ccpppliiiiim wwtppccpcccc////cccccccccccc..lt http//uuiieeeeoeeeeee
for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
In Section
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment Statistical Significance
Does the report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: We provide standard errors in Tables and 2. We perform large number of repetitions of the simulation experiments to calculate the empirical median squared error, so error bars are not relevant in Figures 1, 2, 3, and 4.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emor bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
32

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
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 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 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.
8. Experiments Compute Resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: In Section 5.1.3.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU 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 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 failed experiments that didn't 't make into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification:
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 deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws regulations in their jurisdiction).
10. 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 theoretical research.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
33

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess 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 direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that 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 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 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 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 monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper deseribe 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)?
Question:
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have 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 best faith effort.
12. Licenses for existing assets
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?
Question:
Answer: [Yes]
Justification: URLs for real data provided in Section 5.2.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include 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 use in package should be provided. For popular datasets, paperswithcode. thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
34

---

For existing datasets that are re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/mooel as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included the main paper.
According the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: [NA]
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
35


Reviews:
[{"review_id": "mv3kS0eDcL", "reviewer": "Reviewer_RXcc", "summary": "This paper introduces an optimal subsampling algorithm designed for imbalanced data with inactive features. The primary algorithm can be summarized in three steps:\n\n1. Train a pilot model using a lasso-penalized objective.\n2. Subsample the data based on the pilot model.\n3. Train the subsampled data using an adaptive lasso-penalized algorithm.\n\nAdditionally, the paper proposes two options for the model estimator: the IPW estimator and the more accurate MSCL estimator. Empirical results demonstrate the algorithm's effectiveness, showing that it outperforms the baseline estimator which does not apply variable selection.", "questions": "In the first stage screening of Algorithm 2, we use the Lasso estimator to manage computational costs and noise in the sampled data. If more computational resources are available, is Lasso always the best practical choice for pilot estimation?", "limitations": "No significant limitations were observed.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "1. The paper addresses a key drawback of a previous \"optimal estimator\" for rare-events data by incorporating a lasso-based estimator.\n2. The proposed algorithm is theoretically sound with solid asymptotic guarantees.\n3. Empirical results on both synthetic and real data highlight the algorithm's effectiveness.", "weaknesses": "No technical weaknesses were observed in this paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your positive feedback on our work along with the\ninsightful question. Please see below for our response.\n\n**Q1. If more computational resources are available, is Lasso always the best\npractical choice for pilot estimation?**\n\n**R1.** In general, better pilot estimation leads to a better final\nestimator. While we recommend the Lasso method, it may not be the best choice in\nall cases. More computationally intensive methods can achieve higher accuracy\nthan the Lasso. However, if more computational resources are available, it is\nmore beneficial to allocate additional resources to subsequent steps rather than\nfocusing solely on pilot estimation. Although we cannot assert that the Lasso is\nalways the best choice, it is a robust and practical option for the problem at\nhand.\n\nWe recommend using the Lasso method for pilot estimation due to its simultaneous\nparameter estimation and active set selection. In addition, it tends to be\nconservative, meaning it has a low risk of excluding important variables. Other\nmethods that combine parameter estimation and variable selection, such as sure\nindependence screening, may also be considered. We do not recommend using the\nmaximum likelihood estimator (MLE) under the full model because including many\ninactive variables reduces the efficiency of the approximated optimal\nsubsampling function.\n\nFor further illustration, we conducted a preliminary simulation to compare the performance of the Lasso\npilot (Lasso), Sure Independence Screening (SIS) pilots, and the\nLasso-followed-by-MLE pilot (LassoMLE). The SIS requires to pre-assign the\nnumber of variables to select, for which we considered 3, 6, 13, and 20, and they are labeled as SIS3, SIS6, SIS13, and SIS20, respectively. The LassoMLE calculates the MLE\nbased on the variables selected by the Lasso.  Below are the eMSEs of the final\nP-OS estimator from 200 iterations of the simulation for Section 5.1 Case C of the main paper:\n\n| $\\rho$ | Lasso |  SIS3 |  SIS6 | SIS13 | SIS20 | LassoMLE |\n|-------:|------:|------:|------:|------:|------:|---------:|\n| 0.0025 | 0.031 | 0.735 | 0.026 | 0.039 | 0.049 |    0.053 |\n| 0.0050 | 0.026 |  0.73 | 0.017 | 0.029 | 0.033 |    0.039 |\n| 0.0075 | 0.015 | 0.729 | 0.015 | 0.021 | 0.026 |    0.030 |\n\nThe table above reveals that the SIS pilot may yield better results than the Lasso pilot when the pre-assigned number of variables is 6 (SIS6). However, determining the optimal pre-assigned number of variables in practice requires further investigation. On\naverage, the Lasso pilot selects 13 variables and shows better performance than the SIS13\npilot. Unexpectedly, we found that the LassoMLE pilot performed worse than\nthe Lasso pilot. Upon further examination, we discovered that the LassoMLE produces larger absolute estimates for some inactive variables. Consequently, these inactive variables receive a smaller penalty in the final adaptive Lasso\nestimator, as the penalty is proportional to $1/|\\hat{\\beta}|$.\n\nWe greatly appreciate your insightful question, which has led to these additional interesting findings. The choice of the optimal pilot estimation method is indeed critical in practice, especially for sparse models. We plan to further investigate the effect of different pilot estimation procedures on the performance of the proposed algorithm.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "w40N0TAmXy", "reviewer": "Reviewer_2g5r", "summary": "The scale-invariant optimal subsampling function proposed in the paper addresses the challenge of inactive features in rare-events data by overcoming the issue of scale-dependence in exitsing optimal subsampling methods. In the context of variable selection for rare-events data, where distinguishing active and inactive features is very important, the scale-invariant optimal subsampling function minimizes prediction error regardless of scaling transformations applied to inactive features.\n\nThe key is to provide reliable and unbiased subsampling probabilities even when inactive variables are transformed without altering the underlying model. This is important because inappropriate scaling transformations of inactive features can significantly impact traditional optimal subsampling methods, leading to unreliable or misleading results. This resluts in improving the accuracy of parameter estimation and variable selection.", "questions": "Why is the tr(variance) sufficient ot be minimized in the presented results ? In general most of the discussion revolves around assymptotic results, is it interesting to ask questions about non-assymptotic behaviors, such as how would performance of variable selection impact non-assymptotic prediction error? Further, many practical applications could be overparametrized/high-dimensional, how would the results carry over? If not, what would be complications that arise while generalizing them to such settings?", "limitations": "Yes", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Good novelty:The paper introduces a novel scale-invariant optimal subsampling function to address the challenge of inactive features in rare-events data, offering a new approach to improving parameter estimation and variable selection in sparse moedls. This novel contribution could advance the field of rare-events data analysis.\n\nStrong theory:The paper establishes a theoretical foundation for the proposed scale-invariant optimal subsampling function, including discussions on assumptions, proofs, and the mathematical framework supporting the method. This strong theoretical basis enhances the credibility and rigor of the proposed approach. WHile i have not gone through all the proofs, i think  the results seem sound.\n\nSufficient empirical evidence: Further, the paper demonstrates the practical relevance of the scale-invariant optimal subsampling function through numerical experiments using both simulated and real-world datasets.The results show non-trivial improvements in prediction error minimization and estimation efficiency.\n\nClarity: The paper is largely well-written,, provides clear and detailed explanations of the methodology, assumptions, and implications of the scale-invariant optimal subsampling function, making it accessible even though there is a lot of notation and theory.", "weaknesses": "I think this is a good paper that should be accepted. Minor gripe:\n\nIt is not clear to me how exactly the methods in section 1 are affected due to x (and thereby making them NOT scale invariant). I understand this would require going through the earlier papers e..g [21], which i did not do. It would make the paper complete if the authors could adda discussion to this effect in the appendix in the final version.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful and positive review of our paper. We appreciate your constructive suggestion on providing a detailed discussion on why the A-OS\nand L-OS in Section 1 are affected by the scaling of the covariate\n$\\mathbf{x}$. Yes, we also believe the point would be clarified by the formulas\nfor the A-OS and L-OS probabilities derived from [21] for producing Figure 1. In\nthe final version, we will add these details along with a comprehensive\ndiscussion of the example in Section 1 to the appendix.\n\nPlease see our responses to your specific questions below.\n\n**Q1. Sufficiency of minimizing $tr(\\mathbf{V}_{w(\\mathcal{A})})$.**\n\n**R1.** The subsample estimator is asymptotically unbiased (Theorem 1), implying that lower variance indicates a better subsample estimator. Optimal subsampling aims to improve estimation efficiency by finding subsampling probabilities that minimize variance. However, due to the lack of complete ordering in variance matrices, the trace is often used as a criterion for minimization. For an asymptotically unbiased estimator, the trace ${tr}(\\mathbf{V}_{w(\\mathcal{A})})$ is\ninterpreted as the asymptotic mean squared error (MSE), making its usage appreciated\nin practice.\n\n**Q2. Discussions on non-asymptotic behaviors.**\n\n**R2.** Non-asymptotic behaviors, such as prediction error, are indeed of\ninterest in variable selection. Our numerical experiments actually assess these non-asymptotic performance. However, non-asymptotic theoretical results are typically expressed as error bounds, which may be less applicable for defining optimal subsampling probabilities. The asymptotic distribution, on the other hand, provides a direct\napproximation to the distributional behavior of the subsample estimator, which makes it more suitable for defining optimal subsampling probabilities. We do not intend to downplay the importance of non-asymptotic behaviors. In fact, we may investigate the theoretical non-asymptotic behaviors of our estimator in future research.\n\n**Q3. Applications in overparameterized/high-dimensional problems.**\n\n**R3.** For overparameterized/high-dimensional problems, if the model is sparse\nenough that the true model is low-dimensional, we conjecture that results such\nas the scale-invariant probabilities and the theoretical properties of the MSCL\nestimator still hold under some additional assumptions. However, if the model is\nnot sparse, our results do not carry over. One reason is that asymptotic\nnormality may no longer hold. While it may be possible to derive some\nnon-asymptotic error bounds, it is unclear how to use them to define better\nsubsampling probabilities for the problem in consideration.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "IOSzOQF4Rh", "reviewer": "Reviewer_8W9n", "summary": "This paper studies the problem of learning a sparse model under rare events and scale-invariance. They provide optimal subsampling function to handle scale-invariance for a Lasso-regularized model. They finally provide experimental results for their subsampling scheme.", "questions": "* In Proposition 1, why is minimizing $tr(M_{w(A)})$ of interest? It seems like all one would care about is minimizing the asymptotic variance (which as far as I understand is $V_{w(A)}$).\n* In the preamble of section 3.1 (lines 156-158), the MSPE is motivated by the desire for a scale-invariant objective. Is the problem with the previous approach(es) that scaling $x_j$ would inversely scale the corresponding entries of $V_{w(A)}$?\n* In the preamble of section 4 (lines 183 and 184) that discusses the inefficiency of the IPW estimator, the paper says IPW gives less weight to more informative points. Is this a formal (or at least, formalizable) claim, or is it just a heuristic? Either way, can this be elaborated upon?\n* The optimality of this work seems to be the optimal subsampling function towards the objective of MSPE. Is it possible to argue that this estimator is information-theoretically the best that one can do?", "limitations": "The authors addressed their limitations.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* Well-written\n* Looking at MSPE rather than asymptotic variance seems like a good way to get around the scaling issue.", "weaknesses": "* This problem setting feels very specific: rare events data, sparse models, and scale-invariance. It would be good to motivate why we should care about all of these simultaneously.\n* None of the theoretical results given in the main body seem to have any discussion about proof strategy or novel techniques. It would be interesting to know what technical challenges were posed in developing these results, as well as how they were resolved.\n* In section 4, the paper compares the proposed penalized MSCL estimator to the MSCL estimator. As far as I can tell, this MSCL estimator has not been defined in the paper and so these comparisons and discussions are hard to make sense of.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful review of our paper.\n\nWe would like to point out that rare events data and sparse models are common in\npractice, and scale-dependence is a frequent and crucial issue in subsampling\nthat has not been addressed in existing literature. Following your great\nsuggestion, we will provide more details on the motivation behind our paper's\nsettings in the revised version.\n\nWe did not compare the penalized MSCL estimator to the MSCL estimator. The\npenalized MSCL estimator is defined in equation (9), while the MSCL estimator is\nsimply this without the last penalty term. We did not separately define the MSCL\nas it is not implemented in our paper. We will clarify this in the final\nversion. Furthermore, your comment has made us realize that using the subscript\n\"mscl\" instead of \"lik\" would more clearly indicate the penalized MSCL\nestimator. We will make this change in the revision.\n\nPlease see the following for our response to your specific comments/questions.\n\n**Q1. The reason why minimizing $tr(M_{w(A)})$ is of interest in Proposition 1.**\n\n**R1.** Minimizing $tr(M_{w(A)})$ is of interest because it is the asymptotic\nvariance of estimating $M_{(A)}\\theta_{t(A)}$, a linearly transformed\nparameter. One advantage of this criterion is it often leads to optimal\nsubsampling probabilities that are easier to calculate. We will provide further\nclarification on this point in the final revision.\n\n**Q2. The problem of previous approaches.**\n\n**R2.** Thank you for your insightful question. Yes, the problem with previous\napproaches is essentially caused by the fact that scaling $x_j$ would scale some\nentries of $V_{w(A)}$. However, the resulting effects on the optimal\nprobabilities are more complex than might be expected. For example, Figure 1\ndemonstrates that scaling a single covariate has markedly different effects on\nthe A-OS and L-OS probabilities.\n\n**Q3. Discussion of the inefficiency of the IPW estimator.**\n\n**R3.** The statement that the IPW estimator is not the most efficient due to\nits assignment of lower weights to more informative data points is a heuristic\nexplanation for its inefficiency. More specifically, optimal sampling\nprobabilities are constructed so that more informative data points have higher\nprobabilities $\\pi(\\mathbf{x}_i^{\\mathrm{sub}}, y_i^{\\mathrm{sub}})$'s of being\nincluded in the subsample. However, the IPW estimator assigns smaller weights of\n$1/\\pi(\\mathbf{x}_i^{\\mathrm{sub}}, y_i^{\\mathrm{sub}})$'s. This inverse\nrelationship suggests that the IPW estimator does not fully utilize the\ninformation contained in the optimal subsample, which intuitively explains its\nreduced efficiency.\n\nFurthermore, we can also formally prove that the IPW estimator is less efficient\nthan the penalized MSCL estimator. Here is a sketch of the proof. Let\n\n$$ h = 1 + c \\{\\varphi(\\mathbf{x})\\}^{-1} e^{f(\\mathbf{x}; \\beta_t)}, $$\n\n$$ \\mathbf{v} = \\sqrt{e^{f(\\mathbf{x}; \\beta_t)}} \\dot{g}_{(A)}(\\mathbf{x}; \\theta_t), $$\n\n$$ \\mathbf{f} = h^{1/2} \\mathbf{v}, $$ \n\n$$ \\mathbf{g} = h^{-1/2} \\mathbf{v}. $$\n\nWe have that\n\n$$ \\mathbb{E}(\\mathbf{g} \\mathbf{f}^{T}) = \\mathbb{E}(\\mathbf{f} \\mathbf{g}^{T}) = \\mathbb{E}(\\mathbf{v} \\mathbf{v}^{T}) = M_{(A)}, $$\n\n$$ \\mathbb{E}(\\mathbf{f} \\mathbf{f}^{T}) = M_{w(A)}, $$\n\n$$ \\mathbb{E}(\\mathbf{g} \\mathbf{g}^{T}) = \\Lambda_{mscl(A)}. $$\n\nThen, the matrix form of Cauchy-Schwarz's inequality [ref1],\n\n$$ \\mathbb{E}(\\mathbf{g} \\mathbf{g}^{T}) \\geq \\mathbb{E}(\\mathbf{g} \\mathbf{f}^{T}) \\{\\mathbb{E}(\\mathbf{f} \\mathbf{f}^{T})\\}^{-1} \\mathbb{E}(\\mathbf{f} \\mathbf{g}^{T}), $$\n\nleads to\n\n$$ \\Lambda_{mscl(A)} \\ge V_{w(A)}^{-1}, $$\n\nwhich implies that\n\n$$ V_{mscl(A)} \\le V_{w(A)}. $$\n\nIn the revision, we will elaborate on this point, provide a detailed proof, and\nfurther clarify why the IPW estimator is not the most efficient.\n\n[ref1] Tripathi, Gautam. \"A matrix extension of the Cauchy-Schwarz inequality.\"\nEconomics Letters 63.1 (1999): 1-3.\n\n**Q4. Theoretical optimality of the estimator**\n\n**R4.** The proposed penalized MSCL estimator is indeed the best in the sense\nthat it is the most efficient estimator among a large class of asymptotically\nunbiased estimators. Specifically, the asymptotic variance $V_{mscl(A)}$ attains\nthe Cramer-Rao lower bound for this class of estimators. This result formally\nestablishes the estimator's optimality within the class. We will present this\nfinding in detail, along with a detailed proof, in the future version of this\npaper.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We sincerely thank all the reviewers for their time and effort in assessing our
work. Their insightful comments and questions are valuable in helping us improve
the paper. Please find our individual responses to the reviewers' comments and
questions below.

Meta-review: This paper provides scale-invariant optimal subsampling probabilities that minimize the predictor error of the adaptive lasso for rare binary responses. The proposed algorithm involves two steps: (i) to calculate approximate optimal sampling probabilities using a pilot sample and (ii) to compute an optimally sub-sampled adaptive lasso estimator. For the latter, the paper proposes two options: the inverse probability weighted (IPW) estimator and the more accurate maximum sampled conditional likelihood (MSCL) estimator. It would make the paper more influential if the author(s) can improve the following aspects among other things: (i) the motivation for the estimation problem; (ii) discussion on what sense scale invariance matters for optimal subsampling probabilities; (iii) discussion of the inefficiency of the IPW estimator as well as the optimality of the MSCL estimator; (iv) limitations in terms of finite-sample bounds and the effect of different pilot estimation procedures.
Decision: Accept (poster)