Optimal Sampling for Generalized Linear Models under Measurement Constraints

TL;DR

This paper introduces OSUMC, a novel response-free sampling method for generalized linear models under measurement constraints, optimizing sampling probabilities to improve statistical efficiency without requiring responses for most data points.

Contribution

The paper develops a response-free sampling procedure for GLMs, establishes an unconditional asymptotic distribution, and proposes an algorithm to approximate the optimal sampling distribution under measurement constraints.

Findings

OSUMC outperforms existing methods in empirical studies.

Theoretical guarantees are established for the asymptotic distribution of estimators.

The method effectively balances statistical efficiency and measurement cost.

Abstract

Under "measurement constraints," responses are expensive to measure and initially unavailable on most of records in the dataset, but the covariates are available for the entire dataset. Our goal is to sample a relatively small portion of the dataset where the expensive responses will be measured and the resultant sampling estimator is statistically efficient. Measurement constraints require the sampling probabilities can only depend on a very small set of the responses. A sampling procedure that uses responses at most only on a small pilot sample will be called "response-free." We propose a response-free sampling procedure \mbox{(OSUMC)} for generalized linear models (GLMs). Using the A-optimality criterion, i.e., the trace of the asymptotic variance, the resultant estimator is statistically efficient within a class of sampling estimators. We establish the unconditional asymptotic distribution of a general class of response-free sampling estimators. This result is novel compared with the existing conditional results obtained by conditioning on both covariates and responses. Under our unconditional framework, the subsamples are no longer independent and new martingale techniques are developed for our asymptotic theory. We further derive the A-optimal response-free sampling distribution. Since this distribution depends on population level quantities, we propose the Optimal Sampling Under Measurement Constraints (OSUMC) algorithm to approximate the theoretical optimal sampling. Finally, we conduct an intensive empirical study to demonstrate the advantages of OSUMC algorithm over existing methods in both statistical and computational perspectives.

Equations

$$P(Y|X,{\beta }_0,\sigma )\propto \exp \left\{\frac{Y\cdot X^T{\beta }_0-b(X^T{\beta }_0)}{c(\sigma )}\right\},$$

$${\hat{\beta }}_{\mathrm{mle}}\in \underset{\beta \in {\mathbb{R}}^p}{\mathrm{argmin}}\left[-\frac{1}{n}\sum _{i=1}^n\left\{Y_i\cdot X_i^T\beta -b(X_i^T\beta )\right\}\right].$$

$${\Psi }_n(\beta ):=\frac{1}{n}\sum _{i=1}^n{\psi }_{\beta }(X_i,Y_i)=\frac{1}{n}\sum _{i=1}^n\left\{b^{\prime }(X_i^T\beta )-Y_i\right\}\cdot X_i=0.$$

$${\Psi }_n^{\ast }(\beta ):=\frac{1}{r}\sum _{i=1}^r\frac{b^{\prime }(X_i^{\ast T}\beta )-Y_i^{\ast }}{n{\pi }_i^{\ast }}\cdot X_i^{\ast }$$

$${\Psi }_n^{\ast }({\beta }_0)=-\Phi ({\hat{\beta }}_n-{\beta }_0)+o_p(||{\hat{\beta }}_n-{\beta }_0||).$$

$$M$$

$$V(T{)}^{-\frac{1}{2}}\Phi ({\hat{\beta }}_n-{\beta }_0)\stackrel{d}{⟶}N(0,I).$$

$${\pi }_j=\frac{\sqrt{b^{{}^{\prime \prime }}(X_j^T{\beta }_0)}||{\Phi }^{-1}{X}_j||}{{\sum }_{i=1}^n\sqrt{b^{{}^{\prime \prime }}(X_j^T{\beta }_0)}||{\Phi }^{-1}{X}_j||},$$

$${\pi }_j\propto \sqrt{b^{{}^{\prime \prime }}(X_j^T{\tilde{\beta }}_n)}||{\tilde{\Phi }}_n^{-1}{X}_j||,$$

$$P(Y=1|X,{\beta }_0)=\frac{\exp (X^T{\beta }_0)}{1+\exp (X^T{\beta }_0)},$$

$${\pi }_j\propto ||{\tilde{\Phi }}_n^{-1}{X}_j||,$$

$$\sum _{i=1}^{n}E\left[\frac{\{b^{\prime }(X_i^T\beta )-Y_i{\}}^2}{{\pi }_i}{x}_{ij}^2\right]=o(n^{2}r),$$

$$E{\Psi }_n^{\ast }(\beta )$$

$$=E\left[E\left\{\frac{1}{r}\sum _{i=1}^r\frac{b^{{}^{\prime }}(X_i^{\ast T}\beta )-Y_i^{\ast }}{n{\pi }_i^{\ast }}\cdot X_i^{\ast }|(X_i,Y_i{)}_{i=1}^n\right\}\right]$$

$$=\zeta .$$

$$({\Psi }_n^{\ast }{)}_j$$

$$E[({\Psi }_n^{\ast }{)}_j-{\zeta }_j{]}^2$$

$$E\left[E\left\{{\left[({\Psi }_n^{\ast }{)}_j-{\zeta }_j\right]}^2|(X_i,Y_i{)}_{i=1}^n\right\}\right]$$

$$=E\left[\frac{1}{r}\sum _{i=1}^n{\pi }_i{\left\{\frac{b^{{}^{\prime }}(X_i^T\beta )-Y_i}{n{\pi }_i}\cdot x_{ij}-{\zeta }_j\right\}}^2\right]$$

$$=\frac{1}{r{n}^2}\sum _{i=1}^{n}E\left[\frac{{\left\{b^{{}^{\prime }}(X_i^T\beta )-Y_i\right\}}^2}{{\pi }_i}\cdot x_{ij}^2\right]-\frac{1}{r}{\zeta }_j^2,$$

$$E[({\Psi }_n^{\ast }{)}_j-{\zeta }_j{]}^2=\frac{1}{r{n}^2}\sum _{i=1}^{n}E\left[\frac{(b^{{}^{\prime }}\left(X_i^T\beta \right)-Y_i{)}^2}{{\pi }_i}\cdot x_{ij}^2\right]-\frac{1}{r}{\zeta }_j^2=o(1).$$

$$\left|[{\Psi }_n^{\ast }({\beta }_1)-\zeta ({\beta }_1)]-[{\Psi }_n^{\ast }({\beta }_2)-\zeta ({\beta }_2)]\right|$$

$$=$$

$$\le$$

$$:=$$

$$E\left[\frac{1}{r}\sum _{i=1}^r\frac{X_i^{\ast }{X}_i^{\ast T}}{n{\pi }_i^{\ast }}\right]$$

$$=$$

$$\underset{\beta \in \mathcal{B}}{\sup }\Vert {\Psi }_n^{\ast }(\beta )-\zeta (\beta )\Vert \stackrel{p}{⟶}0.$$

$$\sum _{i=1}^{n}E\left\{\frac{b^{{}^{\prime \prime }}(X_i^T{\beta }_0{)}^2}{{\pi }_i}(x_{ik}{x}_{ij}{)}^2\right\}=o(n^{2}r).$$

$$E{\dot{\Psi }}_n^{\ast }({\beta }_0)$$

$$\displaystyle=E\left[E\left\{\frac{1}{r}\sum_{i=1}^{r}\frac{b^{{}^{\prime\prime}}(X_{i}^{*T}\beta_{0})}{n\pi_{i}^{*}}\cdot X_{i}^{*}X_{i}^{*T}\bigg{|}(X_{i},Y_{i})_{i=1}^{n}\right\}\right]=\Phi.$$

$$({\dot{\Psi }}_n^{\ast }{)}_{kj}$$

$${\Phi }_{kj}$$

$$E{\left[({\dot{\Psi }}_n^{\ast }{)}_{kj}-{\Phi }_{kj}\right]}^2$$

$$\displaystyle E\left[E\left\{\left[(\dot{\Psi}_{n}^{*})_{kj}-\Phi_{kj}\right]^{2}\bigg{|}(X_{i},Y_{i})_{i=1}^{n}\right\}\right]$$

Full text

Keywords: Generalized linear models, Measurement constraints, Unconditional asymptotic distribution, Martingale central limit theorem, A-optimality

1 Introduction

1.1 Motivation and Contribution

Measurement constrained datasets (WangY2017) where only a small portion of the data points have known a response, $Y$, but the covariates, $X$, are available for all data points, are common in practice. Datasets of this type can happen when the response is more expensive or time-consuming to collect than covariates. We will present two motivating examples. For more real-world examples, we refer readers to semi-supervised learning literature (e.g., Zhu2005; Chapelle2006; Chakrabortty2018).

Example 1. (Critical Temperature of Superconductors)  Critical temperature, which is sensitive to chemical composition, is an important property of superconducting materials. Since no scientific model for critical temperature prediction is available (Hamidieh2018), a data-driven prediction model is desirable to guide researchers synthesizing superconducting materials with higher critical temperature. Due to the cost in both money and time for material synthesis, only a small portion of the thousands of potential chemical compounds can be manually tested. So selecting representative compositions to build a statistical model with maximum efficiency is important.

Example 2. (Galaxy Classification)  Galaxy classification is an important task in astronomy (Banerji2010). Visual classification is time-consuming and expensive and is becoming infeasible because the size of astronomical datasets is growing rapidly as more advanced telescopes enter operation. The size of modern galaxy datasets is often of millions or even billions (Reiman2019). It is important to select a representative subsample of galaxies which can be classified accurately by humans, so that an effective classification model can be built.

In addition to many responses being missing, another characteristic of the datasets in these examples is the extreme size which brings huge challenges to statistical computing, data storage and communication. Sampling is a popular approach to super-large datasets where a small portion of the dataset is sampled and used as a surrogate of the entire dataset. There is a large literature on this problem, e.g., Wang2017, but many of the proposed sampling methods assume that the response is known on the entire dataset so are not applicable under measurement constraints. Moreover, they might not increase statistical efficiency even if they were applicable.

The problem addressed in the paper is fitting a generalized linear model (GLM) efficiently to massive measurement-constrained datasets using only a relatively-small subsample obtained by sampling. Measurement constrained datasets require sampling methods to be response-free meaning that the sampling probabilities can only depend on responses of a small pilot sample. Though a huge literature has been devoted to sampling methods for datasets where all records have both responses and covariates, efficient sampling algorithms under measurement constraints, which selects subsamples in a (nearly) response-independent way, are less well-studied. We bridge this gap by proposing a statistically efficient sampling method under measurement constraint (OSUMC) for GLMs which does not require knowledge of responses (except in the pilot sample). Specifically, OSUMC selects subsamples with replacement according to a response-free A-optimal sampling distribution and computes the estimator by solving a weighted score equation based on the selected subsamples. We defer the details to Section 3 and Section 4.

This paper contributes to both theoretical and computational perspectives of the growing literature on statistical sampling. On the theoretical side, we prove asymptotic normality of a general class of response-free sampling estimators without conditioning on the data. Our asymptotic results are significantly different from conditional asymptotics in the traditional sampling literature which condition on both covariates and responses. In particular, the conditional independence of subsamples in the conditional asymptotics is no longer true under our framework. To deal with the correlation in our sampling estimator, we develop novel martingale techniques for our asymptotic theory. Another significant difference is that our asymptotic results compare our sampling estimator to the true parameter rather than to the maximum likelihood estimator as done by Wang2017. Since the true parameter is the object of interest, our results are more informative. Based on the asymptotic theory, we derive the optimal sampling probabilities by minimizing the A-criterion (Khuri2006) i.e., the trace of asymptotic variance matrix, which is equivalent to minimizing the sum of the asymptotic mean squared errors. On the computational side, we propose the OSUMC algorithm to approximate the theoretical A-optimal sampling method. Our sampling algorithm achieves an optimal design by assigning higher probabilities to data points that achieve estimators with higher efficiency. We show in our numerical study that the performance of OSUMC is comparable to that of sampling algorithms which use complete responses to calculate sampling probabilities. This shows that, provided that an efficient algorithm such as OSUMC is used, there is little or no loss in accuracy due to not having the responses available for sampling because of measurement constraints.

1.2 Related Work

A large literature provide numerous variants of subsampling algorithms for linear regression (Drineas2006; Drineas2011; Ma2015; wang2019information). One traditional approach is the leverage sampling which defines the sampling probabilities based on the empirical leverage scores of the design matrix (Huber2011). In more recent literature, leverage sampling serves as an important basis for more refined sampling methods of linear regression (Drineas2006; Drineas2011; Ma2015). These papers fall in a paradigm termed algorithmic leveraging, which is fundamental to randomized numerical linear algebra (RandNLA) which aims at developing fast randomized algorithms for large-scale matrix-based problems. We refer readers to Ma2015 and drineas2016randnla for more references. However, most of the algorithmic leveraging literature is concerned with the numerical performance of algorithms and only a few papers, for example, Ma2015 and raskutti2016statistical, provide statistical guarantees. Ma2015 derives bias and variance for their sampling methods while raskutti2016statistical provide statistical error bounds for their sketching method. In contrast, we will focus on the asymptotic efficiency which has not been considered in the algorithmic leveraging or RandNLA literature. Because leverage sampling for the least-square estimator is response-free it can adapt to measurement constraints, but this is not true for other GLMs.

For other GLMs, Wang2017 developed an A-optimal sampling procedure for logistic regression. Though a similar optimality criterion is adopted, our paper is very different from Wang2017. Most importantly, the sampling procedure in Wang2017 requires the complete set of responses and therefore cannot be used under measurement constraints, while OSUMC is tailored for that setting. From a theoretical perspective, the unconditional asymptotics derived in this paper are much more challenging than the conditional asymptotic setting considered in Wang2017, which conditions on both covariates and responses. The essential conditional independence assumption in the proof of Wang2017 is no longer valid and the correlated structure of the sampling estimator must be treated in our asymptotic theory. New martingale techniques are developed which are significantly different from the techniques used in Wang2017. Besides, unconditional framework allows us to prove the asymptotic results between our sampling estimator and the true parameter, instead of the conditional MLE as in Wang2017. A recent paper by ai2018optimal generalized the results in Wang2017 to other GLMs under a similar framework, so again not applicable to the problem of measurement constraints. In another paper by Ting2018, the authors studied optimality of sampling for asymptotic linear estimators. Their conclusions will reduce to exactly the same results in Wang2017 for logistic regression, and hence cannot deal with measurement constraints neither.

One new research area also dealing with measurement constraints is semi-supervised learning (SSL). SSL attempts to use the unlabeled $X$ ($X$ without $Y$) to improve statistical performance. Though a huge SSL literature is devoted to algorithmic aspects (Zhu2005; Chapelle2006), only a few recent papers studied statistical estimation under the semi-supervised setting (zhang2016semi; cai2018semi; Chakrabortty2018). For example, Chakrabortty2018 estimate linear regression coefficients by regressing imputations of unobserved responses on the corresponding covariates. However, their method is computationally prohibitive for large datasets due to the nonparametric imputation approach while our method is computationally affordable even for massive datasets as demonstrated in out empirical study.

Another closely related area is optimal experiment design which determines the settings of covariates that yield estimators with optimal properties (Khuri2006). Whereas the design is determined freely in classical experiment design (Pukelsheim2006), the design points of this paper must be selected from the original dataset. Also, the solutions of the traditional optimal design problem are often combinatorial, which are computationally infeasible for even moderate size datasets. WangY2017 propose sampling algorithms based on the convex relaxation of the traditional combinatorial problem. However, WangY2017 mainly focus on algorithms and only a mean square error bound is provided for the statistical guarantee. Their results are also proved under traditional conditional framework which is different from the unconditional analysis of this paper. Additionally, independence assumptions for with-replacement sampling in WangY2017 cannot be justified in many real situations, for instance, the two motivating examples. Such assumptions are avoided in our theory. On the practical side, our method offers simple closed forms of the optimal sampling weights while WangY2017’s procedure involves solving a semi-definite programming problem which can be computationally intensive for large datasets.

1.3 Organization of the Paper

The remainder of the paper is organized as follows. In Section 2, we give a brief overview of GLMs. In Section 3, we formulate a general response-free sampling scheme and define the class of sampling estimators of interest. In Section 4, we derive the asymptotic normality of our sampling estimator and define the A-optimal sampling distribution. OSUMC algorithm is proposed to approximate the theoretical optimal sampling procedure. Section 5 compares OSUMC algorithm with several existing sampling methods on both synthetic and real data, where linear and logistic regression models are mainly discussed. We summarize the paper in Section 6. Proofs are in the supplementary materials.

2 Background and Setup

We begin with the background on GLMs. Assume $n$ independent and identically distributed data couples, $(X_{1},Y_{1}),\dots,(X_{n},Y_{n})\sim(X,Y)$, where $X\in\mathbb{R}^{p}$ is a covariate vector, $Y\in\mathbb{R}$ is the response, and $Y$ given $X$ satisfies a GLM with the canonical link:

[TABLE]

Here $b(\cdot)$ is a known function, $\sigma$ is a known dispersion parameter, and $\beta_{0}\in\mathbb{R}^{p}$ is the unknown parameter of interest and assumed to be in a compact set $\mathcal{B}\subseteq\mathbb{R}^{p}$. Without loss of generality, we take $c(\sigma)=1$. The standard estimator of $\beta_{0}$ is the MLE

[TABLE]

Equivalently, we could solve the score equation to obtain the MLE

[TABLE]

Iterative methods such as Newton’s method and its variants are usually adopted to solve such problems numerically (mccullagh1989generalized; aitkin2005statistical). If the sample size $n$ is very large, the computational cost for just one iteration will be huge. Therefore, a sampling approach can be used to reduce computational cost.

In addition, we assume a measurement constraint setting where only a small portion of responses are available initially. As mentioned in the introduction, this setting is common in practice. The main purpose of this paper is to develop a unified and statistically efficient sampling procedure for GLMs under measurement constraints.

3 General Sampling Scheme

We first present in Algorithm 1 a general response-free sampling scheme for GLMs. The class of response-free sampling estimators is defined accordingly.

We emphasize that Algorithm 1 is tailored for the measurement constraints setting, as the sampling weight $\pi_{i}$ only depends on $(X_{1},\dots,X_{n})$ not $(Y_{1},\dots,Y_{n})$ which may not be completely observable. In practice, once the subsample of size $r$ are selected, one need measure only those $r$ responses. This may bring huge economic savings as response measurements are usually expensive in measurement constrained situations.

An additional benefit of sampling is cost savings. In the last step of Algorithm 1, Newton’s method and its equivalent variants are usually adopted (mccullagh1989generalized; aitkin2005statistical). Given the sampling probabilities $\pi$, Algorithm 1 dramatically reduces the computational and storage costs by making them scale in $r$ instead of $n$ which could be much larger than $r$. More concretely, if $n=10^{6}$ and $p=20$, the computational time of one iteration of Newton’s method for the full sample MLE is $O(np^{2})=O(4\times 10^{8})$. In addition, if each data point occupies 1MB of storage space, then the original dataset would occupy around 1TB space. In contrast, for Algorithm 1 with $r=1,000$, the computational time for each iteration is $O(4\times 10^{5})$ and the storage space is less than 1 GB, which substantially lower the computational and storage cost.

The performance of Algorithm 1 depends crucially on the choice of the sampling weights $\pi_{i}$ and the subsample size $r$. Under measurement constraints, subsample size $r$ is determined by the cost of measuring the response and, perhaps, by the availability of the computational or storage resources. The more important question is how to determine the sampling distribution in a data-driven approach such that the resultant sampling estimator achieves optimal efficiency. We will answer this question in the next section by defining the A-optimal sampling distribution based on the asymptotic results therein.

4 Optimal Sampling Procedure and Asymptotic Theory

In this section, we assume the classical asymptotic setting in which $n\to\infty$ and $p$ is fixed. We first show the asymptotic normality of the sampling estimator defined in Algorithm 1, and then find the A-optimal sampling distribution by minimizing the asymptotic mean squared error.

4.1 Notation

The $j$th entry of the covariate vector $X_{i}$ is denoted by $x_{ij}$. For $X\in\mathbb{R}^{p}$, $||X||$ is the Euclidean norm of $X$. We also define tuple notations: $X_{1}^{n}:=(X_{1},X_{2},\dots,X_{n})$ and $Y_{1}^{n}:=(Y_{1},Y_{2},\dots,Y_{n})$. $V(X)$ and $E(X)$ denote the variance and expectation of $X$, respectively.

4.2 Consistency of $\hat{\beta}_{n}$

We first show the statistical consistency of $\hat{\beta}_{n}$.

Theorem 1** (Consistency of $\hat{\beta}_{n}$).**

Assume the following conditions

1. (i)

*Either $(ia)$ $b^{{}^{\prime\prime}}(\cdot)$ is bounded or $(ib)$ X is bounded. * 2. (ii)

$EXX^{T}$* is finite and $\zeta(\beta):=E\left[\{b^{\prime}(X^{T}\beta)-Y\}X\right]$ is finite for any $\beta\in\mathcal{B}$.* 3. (iii)

$\sum_{i=1}^{n}E\left[\frac{\{b^{\prime}(X_{i}^{T}\beta)-Y_{i}\}^{2}}{\pi_{i}}x_{ij}^{2}\right]=o(n^{2}r)$* for $1\leq j\leq p$ and $\beta\in\mathcal{B}$.* 4. (iv)

$\inf_{\beta:||\beta-\beta_{0}||\geq\epsilon}||\zeta(\beta)||>0$* for any $\epsilon>0$.*

Then $\hat{\beta}_{n}\mathop{\to}\limits^{p}\beta_{0}$.

Condition (i$a$) is satisfied for most of GLMs except for Poisson regression. For Poisson regression, our theory is still applicable if Condition (ib) is satisfied. Condition (iii) - (iv) are needed to apply a consistency theorem for M-estimators (Van2000). In particular, condition (iii) ensures the uniform convergence of $\Psi_{n}^{*}$ while condition (iv) is a common well-separated condition for consistency proofs which is satisfied if $\zeta(\cdot)$ has unique minimizer and if, as we have assumed, the parameter space $\mathcal{B}$ is compact.

4.3 Asymptotic Normality of $\hat{\beta}_{n}$

To establish asymptotic normality of $\hat{\beta}_{n}$, we start with an important asymptotic representation.

Lemma 1** (Asymptotic linearity).**

Assume the following conditions,

1. (i)

$\Phi=E\left\{b^{{}^{\prime\prime}}(X^{T}\beta_{0})XX^{T}\right\}$* is finite and non-singular.* 2. (ii)

$\sum_{i=1}^{n}E\left\{\frac{b^{{}^{\prime\prime}}(X_{i}^{T}\beta_{0})^{2}}{\pi_{i}}(x_{ik}x_{ij})^{2}\right\}=o(n^{2}r)$, for $1\leq k,j\leq p$. 3. (iii)

$b(x)$* is three-times continuously differentiable for every $x$ within its domain.* 4. (iv)

Every second-order partial derivative of $\psi_{\beta}(x)$ w.r.t $\beta$ is dominated by an integrable function $\ddot{\psi}(x)$ independent of $\beta$ in a neighborhood of $\beta_{0}$.

If $\Psi_{n}^{*}(\hat{\beta}_{n})=0$ for all large $n$ and if $\hat{\beta}_{n}$ is consistent for $\beta_{0}$, then

[TABLE]

The proof of Lemma 1 is based on the asymptotic linearity of M-estimators, e.g., Theorem 5.41 in Van2000. It turns out to be important for the establishment of asymptotic normality.

We will apply a multivariate martingale central limit theorem (Lemma 4 in the supplementary materials) to the above asymptotic linear representation and show the asymptotic normality of $\hat{\beta}_{n}$.

We first define a filtration $\{\mathcal{F}_{n,i}\}_{i=1}^{r(n)}$ adaptive to our sampling procedure: $\mathcal{F}_{n,0}=\sigma(X_{1}^{n},Y_{1}^{n})$; $\mathcal{F}_{n,1}=\sigma(X_{1}^{n},Y_{1}^{n})\lor\sigma(*_{1})$; $\cdots$; $\mathcal{F}_{n,i}=\sigma(X_{1}^{n},Y_{1}^{n})\lor\sigma(*_{1})\lor\cdots\lor\sigma(*_{i})$; $\cdots$, where $\sigma(*_{i})$ is the $\sigma$-algebra generated by $i$th sampling step, which can be interpreted as the smallest $\sigma$-algebra that contains all the information in $i$th step. In the following, we always assume subsample size $r$ is increasing with $n$. Based on the filtration, we define the martingale

[TABLE]

where $\{M_{i}\}_{i=1}^{r}$ is a martingale difference sequence adapt to $\{\mathcal{F}_{n,i}\}_{i=1}^{r}$. In addition, we define: $Q:=\frac{1}{n}\sum_{j=1}^{n}(b^{\prime}(X_{j}^{T}\beta_{0})-Y_{j})\cdot X_{j}$; $T:=\Psi_{n}^{*}(\beta_{0})=M+Q$; $\xi_{ni}:=V(T)^{-\frac{1}{2}}M_{i}$; $B_{n}:=V(T)^{-\frac{1}{2}}V(M)V(T)^{-\frac{1}{2}}$, which is the variance of the normalized martingale $V(T)^{-\frac{1}{2}}M$.

Theorem 2** (Asymptotic normality of $\hat{\beta}_{n}$).**

Under the conditions in Lemma 1 and we further assume

1. (i)

$\lim\limits_{n\to\infty}\sum\limits_{i=1}^{r}E[||\xi_{ni}||^{4}]=0$, 2. (ii)

$\lim\limits_{n\to\infty}E[||\sum\limits_{i=1}^{r}E[\xi_{ni}\xi_{ni}^{T}|\mathcal{F}_{n,i-1}]-B_{n}||^{2}]=0$,

we have

[TABLE]

Condition (i) and (ii) are martingale version integrability conditions similar to Lindeberg-Feller conditions. In fact, such conditions are common in martingale central limit theorems (hall2014martingale).

4.4 Optimal Sampling Weights under Measurement Constraints

In this section, we will derive the A-optimal sampling distribution for our general response-free sampling procedure.

Theorem 2 shows that for sufficiently large subsample size, the distribution of $\hat{\beta}_{n}-\beta_{0}$ can be well approximated by $N\left(0,\mathbf{V}\right)$ with $\mathbf{V}:=\Phi^{-1}V(T)\Phi^{-1}$. If $\beta_{0}$ is univariate, we can optimize the sampling probability by minimizing the asymptotic variance $\Phi^{-1}V(T)\Phi^{-1}$ which results in highest statistical efficiency. For multi-dimensional $\beta_{0}$, we adopt the A-optimality criterion of experiment design (Kiefer1959) and minimize the trace of the covariance matrix. Minimization of the trace of $\mathbf{V}$, i.e. ${\rm tr}(\Phi^{-1}V(T)\Phi^{-1})$, is equivalent to minimization of the asymptotic mean squared error. The following theorem specifies A-optimal sampling distribution.

Theorem 3**.**

If for $1\leq j\leq n$, the sampling probability is

[TABLE]

then ${\rm tr}(\mathbf{V})$ will attain its minimum, i.e., $\{\pi_{j}\}_{j=1}^{n}$ is the A-optimal sampling distribution.

The optimal weights cannot be calculated directly in practice, since they depend on population level quantities $\Phi^{-1}$ and $\beta_{0}$. Therefore, to implement response-free sampling, we need pilot estimates of $\Phi$ and $\beta_{0}$. The details are shown in Algorithm 2.

Remarks

• •

Step 1 in Algorithm 2 is designed for pilot estimation when very few or even none of the responses are available initially, but expensive responses collection is possible. If a moderate number of responses are accessible in the initial data pool, we could use them to calculate pilot estimators. Otherwise, a small random sample can be taken with uniform sampling probabilities. The size of the pilot sample, $r_{0}$ could be fairly small compared with total sample size $n$ and even of the same order of magnitude as the dimension $p$ of the problem. In our empirical study, we set $r_{0}=500$ with the complete dataset of size $n=10^{5}$ and $p$ ranging from $20$ to $100$, and the performance of Algorithm 2 is satisfactory.

• •

(Computational complexity) When we use Newton’s method, or one of its variants such as Fisher scoring, to compute the root of the score equation, the computation requires $O(\zeta np^{2})$ computational time, where $\zeta$ is the number of iterations needed for the algorithm to converge. In our empirical study, $\zeta$ varies from $10$ to $30$ under different models. For OSUMC algorithm the first step requires $O(\zeta_{1}r_{0}p^{2})$ computation time where $\zeta_{1}$ is the number of iterations. In the second step, $O(np^{2}+\zeta_{2}rp^{2})$ computation time is required where $\zeta_{2}$ is the number of the iterations in this step. Hence, OSUMC algorithm has complexity of order $O(np^{2}+\zeta_{1}r_{0}p^{2}+\zeta_{2}rp^{2})$. If $n$ is extremely large such that $p$, $r_{0}$, $r$, $\zeta_{1}$ and $\zeta_{2}$ are all much smaller than $n$, the computation complexity of the algorithm is $O(np^{2})$. Therefore, the computational advantage of using OSUMC algorithm compared with full-sample MLE is still huge if the scale of the problem is large, i.e., $np^{2}$ is large and $\zeta>10$. The intensive numerical study in the following provides strong evidence for such advantage. To further reduce the computational complexity, one may use modified Newton’s method for large-scale computation; see, for example, xu2016sub.

5 Numerical Examples

5.1 Simulation Results

In this section, we evaluate the performance of the OSUMC algorithm on synthetic datasets. Due to page limitation, we will show the numerical results for logistic and linear regression and defer Poisson regression to Section A6 in the supplementary material. All the results are obtained in the R environment with one Intel Xeon processor and 8 Gbytes RAM over Red Hat OpenStack Platform.

5.1.1 Logistic Regression

We generate datasets of size $n=100,000$ from the logistic regression model,

[TABLE]

where $\beta_{0}$ is a $20$ dimensional vector with all entries $1$. We consider four different scenarios to generate $X$ as in Wang2017.

• •

mzNormal. $X$ follows the multivariate normal distribution $N(0,\Sigma)$ with $\Sigma_{ij}=0.5^{I(i\neq j)}$. In this case, we have a balanced dataset, i.e., the number of 1’s and the number of 0’s in the responses are almost equal.

• •

nzNormal. $X$ follows the multivariate normal distribution $N(0.5,\Sigma)$. In this case, we have an imbalanced dataset where about $75\%$ of the responses are 1’s.

• •

unNormal. $X$ follows the multivariate normal distribution with mean zero but different variances. To be more specific, $X$ follows the multivariate normal distribution $N(0,\Sigma_{1})$ with $\Sigma_{1}=U_{1}\Sigma U_{1}$, where $U_{1}=diag(1,1/2,\dots,1/20).$

• •

mixNormal $X\sim 0.5N(0.5,\Sigma)\ +\ 0.5N(-0.5,\Sigma)$.

In each case, we compare our optimal sampling procedure (OSUMC) with the mMSE method in Wang2017 (OSMAC), uniform sampling (Unif), and the benchmark full data MLE under different subsample sizes ($r$). In our procedure, we set the subsample size of uniform sampling in the first step equal to $r_{0}=500$. For uniform sampling, we directly subsample $r$ points and calculate the subsample MLE.

We repeat simulations $S=500$ times, and calculate the empirical MSE as $S^{-1}\sum_{s=1}^{S}||\hat{\beta}_{n}^{(s)}-\beta_{0}||^{2}$ where $\hat{\beta}_{n}^{(s)}$ is the estimate from the $s$th repetition. The comparison of the empirical MSE is presented in Figure 1.

From Figure 1, both OSUMC method and the OSMAC method in Wang2017 uniformly dominate the uniform sampling method in all four scenarios, which agrees with Theorem 3. In most of the simulation settings (except for unNormal), our sampling procedure performs similarly to the OSMAC in Wang2017. This is because both methods adopt the A-optimality criterion in respective framework to derive the sampling weights. However, we note that OSMAC requires the response of each data point, which is infeasible under measurement constraints, while our method can be implemented as long as a moderate number of responses are available for the pilot estimators.

We also compare the average computational time for each method under all scenarios and the computational time plot can be found in Section A5.1 in the supplementary materials. Our simulation reveals that the computation time is not very sensitive to the subsample size for all the four methods. In most cases, OSUMC and OSMAC perform similarly and require significantly less computational time compared with the full-sample MLE.

To see whether the asymptotic normality in our theory holds under the previous four different design generation settings, we plot the chi-square Q-Q plot of the resultant estimator $\hat{\beta}_{n}$ based on $1000$ simulations with fixed subsample size $r=5000$ for each considered setting. The plots are presented in the supplementary materials, Section A5.2. Q-Q plots reveal that the resultant sampling estimator $\hat{\beta}_{n}$ is approximately normal with sufficiently large sample size $n$ and subsample size $r$ in the four considered design generation settings.

5.1.2 Linear Regression

We generate datasets of size $n=100,000$ and dimension $p=100$ from the following linear regression model: $Y=X\beta_{0}+\epsilon$, where $\beta_{0}=(\underbrace{0.1,\dots,0.1}_{5},\underbrace{10,\dots,10}_{90},\underbrace{0.1,\dots,0.1}_{5})^{T}$ and $\epsilon\sim N(0,9I_{n})$.

We note that in linear regression model, OSUMC algorithm is equivalent to the following Algorithm 3. Algorithm 3 is similar to the general least-squares sampling meta-algorithm in Ma2015, which is adopted in Drineas2006; Drineas2011; Drineas2012.

We consider the following design generation settings from Ma2015. Similar settings are also investigated in WangY2017.

GA. The $n\times p$ design matrix $\mathbf{X}$ is generated from multivariate normal $N(1_{p},\Sigma_{2})$ with $\Sigma_{2}=U_{2}\Sigma U_{2}$, where $U_{2}={\rm diag}(5,5/2,\dots,5/30)$. 2. 2.

$\mathbf{T_{3}.}$ Design matrix $\mathbf{X}$ is generated from multivariate t-distribution with 3 degrees of freedom and covariance $\Sigma_{2}$ as GA. 3. 3.

$\mathbf{T_{1}.}$ Design matrix $\mathbf{X}$ is generated from multivariate t-distribution with 1 degrees of freedom and covariance $\Sigma_{2}$ as GA.

We compare our method (OSUMC) with uniform sampling (Unif), leverage sampling (Leverage) and shrinkage leveraging (SLEV) in Ma2015 over different subsample size in the three design generation settings above. For the SLEV method, the shrinkage parameter $\alpha$ is set to be $0.9$ as in Ma2015. Again, we repeat the simulation $500$ times and report the empirical MSE and computational time in Figures 2 and A3, respectively.

For all three design generation settings, our method always results in smaller MSE than the other three methods, which again is consistent with our theoretical results. The advantage of our method becomes more obvious when the design generation distribution is more heavy-tailed by noting that $GA$ has moments of arbitrary orders while $T_{k}$ only has moments up to order $k-1$. It is interesting to see that our method outperforms other methods even in the $T_{1}$ design setting where the moment assumptions imposed in our theory are violated. The performance of both leverage sampling and shrinkage leverage sampling improves with heavier-tailed design generation distributions. This has been well understood in the literature on leverage sampling and outlier diagnosis (Rousseeuw2011; Ma2015). As expected, all the three sampling methods above yield smaller MSEs than uniform sampling.

The average computational times of the four methods are reported in Figure A3 in the supplementary materials. Again, the results show the insensitivity of the computational time to increasing subsample sizes. Our method requires the second smallest computing time, being inferior only to the uniform sampling. Both leverage related methods take more than double the computational time of our method due to the intensive computation of leveraging score of each data point.

Again, we provide Q-Q plots of the resultant estimator $\hat{\beta}_{n}$ for each considered design setting and the results can be found in Section A5.4 in the supplementary materials.

5.2 Superconductivty Dataset

In this section, we analyze the superconductivty dataset (Hamidieh2018), which is available from the Machine Learning Repository at: https://archive.ics.uci.edu/ml/datasets/Superconductivty+Data#. The purpose of Hamidieh2018 is to build a statistical model to predict the superconducting critical temperature of superconducting materials based on their chemical formulas. In the dataset, 21,263 different superconductors’ critical temperatures are collected along with $81$ features extracted from the chemical formulas the superconductors. Multiple linear model is considered in Hamidieh2018 and regression coefficients are calculated based on the full sample, which is treated as “true” parameters ($\beta_{0}$) in the following analysis.

We compare our sampling method (OSUMC) with the three other sampling methods in the simulation studies for linear models. Besides the estimation accuracy which is the main focus before, prediction performance of the sampling algorithm will also be evaluated. We randomly select 19,000 data points as the training set and use the rest as the test set for prediction purpose. Then we implement the sampling method on the training dataset and obtain the coefficient estimator $\hat{\beta}$. We now measure the estimation and prediction performance by relative mean squared error: $||\hat{\beta}-\beta_{0}||^{2}/||\beta_{0}||^{2}$ and prediction relative squared error: $||X\hat{\beta}-Y||^{2}/||X\beta_{0}-Y||^{2}$ which is calculated over test dataset, respectively. We repeat the process $500$ times for different subsample sizes and the median of the two criteria are recorded for each subsample size. The results are presented in Figure 3. We also report the median running times of the four sampling methods over different subsample sizes in Figure A5.

Figure 3 shows that our method consistently achieves the best performance on both estimation and prediction. In addition, our sampling procedure outperforms both leverage-based methods in computational time. Though SLEV method achieves similar prediction accuracy as our method, it takes more than twice the computational time. Therefore, the proposed sampling procedure maintains a good balance between statistical efficiency and computational cost.

6 Conclusion

In this paper, we propose a novel sampling procedure, OSUMC, to address measurement constraints when estimating GLMs. We show unconditional asymptotic normality of a general class of response-free sampling estimators by using newly developed martingale techniques. Our unconditional asymptotic results obtained without conditioning on the data are different from existing conditional results which condition on both covariates and responses. Building on the asymptotic theory, we derive the A-optimal sampling distribution, which depends on population level quantities. For practical applications, we propose OSUMC algorithm to approximate the theoretical optimal sampling scheme. Additionally, we conduct extensive numerical studies which show that the performance of OSUMC is comparable to that of sampling algorithms which use complete responses to calculate sampling probabilities. This indicates that OSUMC successfully prevents the loss of statistical efficiency due to measurement constraints.

A number of interesting extensions and open areas remain. For example, it would be interesting to consider the high-dimensional scenario where the dimension $p$ could be much larger than the subsample size $r$. Under such setting, techniques like regularization and debiasing in the high-dimensional statistics would be need which is out of the scope of this paper. We will leave it for future investigations.

Supplementary Materials

Appendices A supplementary PDF file with proofs of all the theoretical results in Section 4 and additional plots and simulation results for Section 5.

R code Code for the simulations and real data analysis in Section 5.

Dataset The Superconductivity dataset used for real data analysis.

Acknowledgements

Ning was supported in part by NSF grant DMS-1854637.

A1 Extra Notation

For a matrix $A\in\mathbb{R}^{p\times p}$, $\lambda_{\rm max}(A)$ denotes the largest eigenvalue of matrix $A$ and $||A||$ is the Frobenius norm of $A$. For two positive definite matrices, $B$ and $C$, we write $B>C$ if and only if $B-C$ is positive definite (p.d.). In particular, $B$ is p.d. if and only if $B>0$.

A2 Proof of Theorem 1

To show the consistency of $\hat{\beta}_{n}$, we will need the following lemma.

Lemma 2**.**

For any $\beta\in\mathcal{B}$, if $\zeta(\beta):=E\left[\{b^{\prime}(X^{T}\beta)-Y\}X\right]$ is finite and the following condition holds

[TABLE]

then we have $\Psi_{n}^{*}(\beta)-\zeta(\beta)\mathop{\longrightarrow}\limits^{p}0$ for any $\beta\in\mathcal{B}$.

Proof.

Observe that

[TABLE]

We will use Chebyshev’s inequality to show convergence in probability. We denote the $j$-th element in the vector $(\Psi_{n}^{*})$ by

[TABLE]

and the $j$-th coordinate of $\zeta$ as $\zeta_{j}$.

By using Chebyshev’s inequality, it suffices to show $E[(\Psi_{n}^{*})_{j}-\zeta_{j}]^{2}=o(1)$, where

[TABLE]

For this expectation

[TABLE]

where the first equality is based on the fact that after conditioning on the $n$ data points, the $r$ repeating sampling steps should be independent and distributionally identical in each step. And we use $E_{*}$ to denote expectation with respect to sampling randomness. Hence, we have

[TABLE]

The second equality is due to the assumption (A2.1). ∎

We now show the consistency of $\hat{\beta}_{n}$.

We will verify the conditions in Theorem 5.9 in Van2000 and apply the theorem to prove the consistency of $\hat{\beta}_{n}$.

First of all, for any $\beta_{1}$ and $\beta_{2}$ in $\mathcal{B}$,

[TABLE]

The first step is due to mean value theorem with $\tilde{\beta}_{1}$ and $\tilde{\beta}_{2}$ lying on the segment between $\beta_{1}$ and $\beta_{2}$.

We now show $L_{n}=O_{p}(1)$. By assumption (i), it suffices to show $\frac{1}{r}\sum\limits_{i=1}^{r}\frac{X_{i}^{*}X_{i}^{*T}}{n\pi_{i}^{*}}=O_{p}(1)$. This is true because we have

[TABLE]

and it follows from Markov inequality. Now we apply Lemma 2.9 in NEWEY19942111 to conclude $\Psi^{*}_{n}(\beta)-\zeta(\beta)$ is stochastic equicontinuous. Again, by Theorem 21.9 in Davidson1994, Lemma 2 and stochastic equicontinuity imply

[TABLE]

This uniform convergence condition together with condition (iv) in the theorem yield the desired conclusion by applying Theorem 5.9 in Van2000. ∎

A3 Proof of Theorem 2

In this section, we will establish the asymptotic normality of $\hat{\beta}_{n}$. Let us start will some auxiliary lemmas.

A3.1 Proof of Lemma 1

We first prove following lemma.

Lemma 3**.**

Assume that $\Phi=E\left\{b^{{}^{\prime\prime}}(X^{T}\beta_{0})XX^{T}\right\}$ is finite and non-singular. Further assume for $1\leq k,j\leq p$,

[TABLE]

Then, we have $\dot{\Psi}_{n}^{*}(\beta_{0})=\frac{1}{r}\sum_{i=1}^{r}\frac{b^{{}^{\prime\prime}}(X_{i}^{*T}\beta_{0})}{n\pi_{i}^{*}}\cdot X_{i}^{*}X_{i}^{*T}\mathop{\longrightarrow}\limits^{p}\Phi$.

Proof.

First we note

[TABLE]

To show the convergence in probability, we use Chebyshev’s inequality. Consider each element in the matrix

[TABLE]

By using Chebyshev’s inequality, it suffices to show $E\left[(\dot{\Psi}_{n}^{*})_{kj}-\Phi_{kj}\right]^{2}=o(1)$.

[TABLE]

For this expectation

[TABLE]

where the first equality is based on the fact that after conditioning on the $n$ data points, the $r$ repeating sampling steps should be independent and distributionally identical in each step, and we use $E_{*}$ to denote expectation with respect to sampling randomness. Hence, we have

[TABLE]

The second equality is due to the assumption (A3.2). ∎

Now we prove Lemma 1.

By Taylor’s Theorem:

[TABLE]

where $\tilde{\beta}_{n}$ is on the line segment between $\beta_{0}$ and $\hat{\beta}_{n}$. $\ddot{\Psi}_{n}^{*}$ is a k-vector of $(k\times k)$ matrices.

We now show that $\left\|\ddot{\Psi}^{*}_{n}(\tilde{\beta}_{n})\right\|=O_{p}(1)$. By assumption (iv)

[TABLE]

The last equality is because of the fact

[TABLE]

and application of Markov inequality.

Therefore,

[TABLE]

This implies the conclusion

[TABLE]

∎

A3.2 Multivariate martingale CLT

Now, we prove a multivariate extension of the martingale central limit theorem stated in Ohlsson1989 Theorem A.1, which will be appropriate for our with replacement sampling setting.

Lemma 4** (Multivariate version of martingale CLT).**

For $k=1,2,3,\dots$, let $\{\xi_{ki};i=1,2,\dots,N_{k}\}$ be a martingale difference sequence in $\mathbb{R}^{p}$ relative to the filtration $\{\mathcal{F}_{ki};i=0,1,\dots,N_{k}\}$ and let $Y_{k}\in\mathbb{R}^{p}$ be an $\mathcal{F}_{k0}$-measurable random vector. Set $S_{k}=\sum\limits_{i=1}^{N_{k}}\xi_{ki}$. Assume the following conditions.

1. (i)

$\lim\limits_{k\to\infty}\sum\limits_{i=1}^{N_{K}}E\left[\left\|\xi_{ki}\right\|^{4}\right]=0$** 2. (ii)

$\lim\limits_{k\to\infty}E\left[\left\|\sum\limits_{i=1}^{N_{k}}E\left[\xi_{ki}\xi_{ki}^{T}|\mathcal{F}_{k,i-1}\right]-B_{k}\right\|^{2}\right]=0$* for some sequence of positive definite matrices $\{B_{k}\}_{k=1}^{\infty}$ with $\sup\limits_{k}\lambda_{max}(B_{k})<\infty$ i.e. the largest eigenvalue is uniformly bounded.* 3. (iii)

For some probability distribution $L_{0}$, $*$ denotes convolution and $L(\cdot)$ denotes the law of random variables:

[TABLE]

Then we have

[TABLE]

Proof.

We use Cramer-Wold device to deduce it from the univariate case. For any $a\in\mathbb{R}^{p}$, by Cramer-Wold device, it suffices to show

[TABLE]

We check the conditions of Theorem A.1 in Ohlsson1989.

$\sum\limits_{i=1}^{N_{k}}E\left[\left(a^{T}\xi_{ki}\right)^{4}\right]\leq\sum\limits_{i=1}^{N_{k}}||a||^{4}\cdot E\left[||\xi_{ki}||^{4}\right]=||a||^{4}\sum\limits_{i=1}^{N_{k}}E[||\xi_{ki}||^{4}]\longrightarrow 0$ The inequality is due to Cauchy-Schwarz inequality. 2. 2.

[TABLE] 3. 3.

[TABLE]

Here we use $\phi_{*}(t)$to denote the characteristic function. Hence

[TABLE]

From above verification, we use Theorem A.1 in Ohlsson1989 to obtain

[TABLE]

And by Cramer-Wold device, this finishes the proof. ∎

A3.3 More Auxiliary Results

Lemma 5**.**

$\{M_{i}\}_{i=1}^{r}$* is a martingale difference sequence relative to the filtration $\{\mathcal{F}_{n,i}\}_{i=1}^{r}$.*

Proof.

The $\mathcal{F}_{n,i}$-measurability follows from the definition of $M_{i}$ and the definition of the filtration $\{\mathcal{F}_{n,i}\}_{i=1}^{r}$. And we also have

[TABLE]

Combine these two results, we finish the proof. ∎

With Lemma 4, we could easily get the following result.

Corollary 1**.**

$V(T)=V(M)+V(Q).$**

Lemma 6**.**

$\sup\limits_{n}\lambda_{max}(B_{n})\leq 1.$**

Proof.

Since $B_{n}$ is symmetric, it suffices to show for any $n$, $I-B_{n}$ is positive definite.

[TABLE]

Therefore, $I-B_{n}$ is congruent to matrix $V(Q)$ which is positive definite. Hence $I-B_{n}$ is also positive definite and this finishes the proof. ∎

Lemma 7** (Asymptotic normality of $\Psi_{n}^{*}(\beta_{0})$).**

Assume the following conditions

1. (i)

$\lim\limits_{n\to\infty}\sum\limits_{i=1}^{r}E\left[||\xi_{ni}||^{4}\right]=0.$** 2. (ii)

$\lim\limits_{n\to\infty}E\left[\left\|\sum\limits_{i=1}^{r}E[\xi_{ni}\xi_{ni}^{T}|\mathcal{F}_{n,i-1}]-B_{n}\right\|^{2}\right]=0.$**

Then we will have

[TABLE]

Proof.

We verify the conditions in Lemma 4 with

[TABLE]

By Lemma 5, conditions (i) and (ii), we can easily see the first two conditions of Lemma 4 are satisfied. It suffices to show the third condition in Lemma 4 holds. We first note the following conclusion

[TABLE]

This because $Q$ is a sum of i.i.d mean zero random variables, $(b^{\prime}(X_{j}^{T}\beta_{0})-Y_{j})\cdot X_{j}$, which have finite variance and a simple application of central limit theorem will give the above conclusion.

Now, we verify the third condition. For any $t\in\mathbb{R}^{p}$

[TABLE]

The first equality is due to Lemma 8 in the following. Therefore, we have verified the third condition in Lemma 4. And by that lemma we have

[TABLE]

∎

Now we state the following lemma that has been used in the proof of previous lemma.

Lemma 8**.**

Under conditions in Lemma 7 For any $t\in\mathbb{R}^{p}$,

[TABLE]

as $n\to\infty$ and $A_{0}\sim N(0,I)$.

Proof.

Since $V(Q)^{-\frac{1}{2}}Q\mathop{\longrightarrow}\limits^{d}N(0,I)$, for any $\xi\in\mathbb{R}^{p}$,

[TABLE]

as $n\to\infty$. And the convergence is uniform in any finite set of $\xi$. (see Chapter 6 of Chung2001). By setting $\xi=V(Q)^{\frac{1}{2}}V(T)^{-\frac{1}{2}}t^{T}$, to prove the lemma, it suffices to show

[TABLE]

for any fixed $t$. Also we note that

[TABLE]

Hence, it is enough to show $\lambda_{max}\left(V(Q)^{\frac{1}{2}}V(T)^{-\frac{1}{2}}\right)\leq 1$. For notation simplicity, we denote $A=V(Q)$ and $B=V(T)$ in the following proof of the lemma.Note the following equation holds

[TABLE]

That is $A^{\frac{1}{2}}B^{-\frac{1}{2}}$ is similar to $B^{-\frac{1}{4}}A^{\frac{1}{2}}B^{-\frac{1}{4}}$. Therefore, we only need to show $\lambda_{max}\left(B^{-\frac{1}{4}}A^{\frac{1}{2}}B^{-\frac{1}{4}}\right)\leq 1$. This is implied by the fact

[TABLE]

In fact

[TABLE]

That is, $I-B^{-\frac{1}{4}}A^{\frac{1}{2}}B^{-\frac{1}{4}}$ is congruent to $B^{\frac{1}{2}}-A^{\frac{1}{2}}$. Therefore, it suffices to show $B^{\frac{1}{2}}-A^{\frac{1}{2}}$ is positive definite.

Note $B>A$. This is because $B-A=V(M)>0$. Now we use Theorem 1.1 in Zhan2004, then we will have $B^{\frac{1}{2}}-A^{\frac{1}{2}}>0$ which finishes proof. ∎

Now we are able to prove our theorem 2 which shows the asymptotic normality of the sampling estimator $\hat{\beta}_{n}$.

A3.4 Proof of Theorem 2

Proof.

By Lemma 1

[TABLE]

Now we normalize both sides with $V(T)^{-\frac{1}{2}}$

[TABLE]

By Lemma 7

[TABLE]

∎

A4 Proof of Theorem 3

First of all, we condition on $X_{1}^{n}$(or consider $X_{1}^{n}$ is fixed). We now find out $V(T|X_{1}^{n})=V(\Psi_{n}^{*}(\beta_{0})|X_{1}^{n})$. We have

[TABLE]

Here $E_{Y}$ means we take expectation w.r.t randomness of $Y$ after we conditioning on $X$. After some simple calculation, we could get

[TABLE]

Hence, we have

[TABLE]

We now minimize $tr(\Phi^{-1}V(T|X_{1}^{n})\Phi^{-1})$

[TABLE]

where in the last step we use Cauchy-Schwarz inequality and the equality holds iff $\pi_{j}\propto\sqrt{b^{{}^{\prime\prime}}(X_{j}^{T}\beta_{0})}\left\|\Phi^{-1}X_{j}\right\|$.

Now we consider $V(T)$ under random design.

[TABLE]

However, we could verify that in GLM

[TABLE]

Therefore, we have $V(T)=E[V(T|X_{1}^{n})]$. From this we have

[TABLE]

A5 Additional Plots of Section 5

A5.1 Computational Time Plots for Logistic Regression

Figure A1 reveals that the computation time is not very sensitive to the subsample size for all the four methods. All of the three sampling methods outperform the full sample MLE. It is not surprising to see the uniform sampling always takes the least computation time, since it does not involve the computation of the sampling probability to compensate for the loss of efficiency. In most cases, OSUMC and OSMAC require significantly less computational time compared with full-sample MLE.

A5.2 Q-Q Plots for Logistic Regression

To see whether the asymptotic normality in our theory implies approximate finite-sample normality under the previous four different design generations in logistic regression model, we plot the chi-square Q-Q plot of the resultant estimator $\hat{\beta}_{n}$ for each considered setting. Here, we replace the approximated optimal sampling weight in Algorithm 2 with the oracle optimal weights to calculate the estimator $\hat{\beta}_{n}$, i.e., the true $\beta_{0}$ is used in the calculation of optimal sampling weights. Experiments are repeated $1000$ times under each setting and corresponding Q-Q plots are presented in the Figure A2. Nearly all the points lie on the straight line in each plot, which is consistent with $\hat{\beta}_{n}$ being approximately normally distributed in the four considered design generation settings.

A5.3 Computational Time Plots for Linear Regression

From Figure A3, Again, the results show the insensitivity of the computational time to increasing subsample sizes. Our method requires the second smallest computing time, being inferior only to the uniform sampling. Both leverage related methods take more than double the computational time of our method due to the intensive computation of leveraging score of each data point.

A5.4 Q-Q Plots for Linear Regression

To explore further how sensitive the approximate finite-sample normality is to the moment condition of the design distribution in linear regression setting, we show chi-square Q-Q plots for several design generation distributions with different orders of moment. To be more specific, GA, $T_{9}$, $T_{3}$ and $T_{1}$ distributions are considered. Experiments are repeated $1000$ times under each setting and results are presented in the following.

As shown in Figure A4, the resultant sampling estimator $\hat{\beta}_{n}$ is approximately normal in GA, $T_{9}$ and $T_{3}$ settings where we should note that the multivariate t-distribution with 3 degrees of freedom doesn’t even have a third moment. This indicates that the normality of $\hat{\beta}_{n}$ in linear models holds under very weak moment conditions for the design-generation distribution. One surprising fact is that OSUMC outperforms the other sampling methods in the $T_{1}$ setting, even though normality fails to hold.

A5.5 Computational Time Plot for Superconductivity Data Set

A6 Simulation for Poisson Regression

We generate datasets of size $n=100,000$ from the following Poisson regression model,

[TABLE]

where $\beta_{0}$ is a $100$ dimensional vector with all entries $0.5$. We consider two different scenarios to generate $X$.

• •

Case 1. each covariate of $X$ follows independent uniform distribution over $[-0.5,0.5]$.

• •

Case 2. First half of covariates of $X$ follow independent uniform distribution over $[-0.5,0.5]$ while the other half follow independent uniform distribution over $[-1,1]$.

In each case, we compare our optimal sampling procedure (OSUMC) with uniform sampling (Unif), and the benchmark full data MLE under different subsample sizes. In our procedure, $r_{0}$, the subsample size in the first step uniform sampling, equals $500$. For uniform sampling, we directly subsample $r$ points and calculate the subsample MLE. Again, we repeat the simulation $500$ times and report the empirical MSE and computational time in Figures A6 and A7, respectively.

From Figure A6, OSUMC method uniformly dominates the uniform sampling method in both scenarios in terms of mean squared errors, which supports the A-optimality of OSUMC. For average computational time, our simulation reveals that the computation time is not very sensitive to the subsample size. In both cases, OSUMC requires significantly less computational time compared with full-sample MLE.

To see whether the asymptotic normality in our theory holds under the both design generation settings, we plot the chi-square Q-Q plot of the resultant estimator $\hat{\beta}_{n}$ for each considered setting and the results are presented in Figure A8. Q-Q plots reveal that $\hat{\beta}_{n}$ is approximately normal with sufficiently large sample size $n$ and subsample size $r$ in the both considered design generation settings, which again provides empirical support for our theoretical results, especially the discussion of Poisson regression after Theorem 1.