\begin{document}


\def\spacingset#1{\renewcommand{\baselinestretch}
{#1}\small\normalsize} \spacingset{1}


\if1\blind
{
  \title{\bf Optimal Subsampling Algorithms for Big Data Generalized Linear Models}
  \author{Mingyao Ai, Jun Yu, Huiming Zhang and HaiYing Wang\footnote{
      Mingyao Ai is Professor,  LMAM, School of Mathematical Sciences and Center for Statistical Science, Peking University, Beijing 100871, China (myai@math.pku.edu.cn);
      Jun Yu is PhD student, School of Mathematical Sciences and Center for Statistical Science, Peking University, Beijing 100871, China (yujunstd@pku.edu.cn);
      Huiming Zhang is PhD student, School of Mathematical Sciences and Center for Statistical Science, Peking University, Beijing 100871, China (zhanghuiming@pku.edu.cn);
      and HaiYing Wang is Assistant Professor, Department of Statistics, University of Connecticut, Storrs, Mansfield, CT 06269, USA (haiying.wang@uconn.edu).
      Correspondence should be addressed to HaiYing Wang.
    The authors would like to thank Dr. Jinzhu Jia for his helpful suggestions and discussions. The work is partially supported by NSFC grants 11271032 and 11331011, and LMEQF.}
}

  \maketitle
} \fi

\if0\blind
{
  \bigskip
  \bigskip
  \bigskip
  \begin{center}
    {\LARGE\bf Optimal Subsampling Algorithm for Big Data Generalized Linear Models}
\end{center}
  \medskip
} \fi

\bigskip
\begin{abstract}
  To fast approximate the maximum likelihood estimator with massive data, Wang \emph{et al.} (\emph{JASA}, 2017) proposed 
an {\bf o}ptimal {\bf s}ubsampling {\bf m}ethod under the {\bf A}-optimality {\bf c}riterion (OSMAC) for in logistic regression. 
This paper extends the scope of the OSMAC framework to include generalized linear models with canonical link functions.
The consistency and asymptotic normality of the estimator from a general subsampling algorithm are established, and optimal subsampling probabilities under the A- and L-optimality criteria are derived. Furthermore, using Frobenius norm matrix concentration inequality, finite sample properties of the subsample estimator based on optimal subsampling probabilities are derived. Since the optimal subsampling probabilities depend on the full data estimate, an adaptive two-step algorithm is developed. Asymptotic normality and optimality of the estimator from this adaptive algorithm are established.
The proposed methods are illustrated and evaluated through numerical experiments on simulated and real datasets.
\end{abstract}

\noindent
{\it Keywords:}  generalized linear models; massive data; matrix concentration inequality; optimal subsampling; OSMAC.
\vfill

\newpage
\spacingset{1.45} 
\section{Introduction} \label{sec:intro}

Nowadays, massive data sets are ubiquitous in many scientific fields and practices such as in astronomy, economics, and industrial problems. Extracting useful information from these large data sets is a core challenge for different communities including computer science, machine learning, and statistics. Over last decades, progresses have been make through various investigations to meet this challenge. However, computational limitations still exist due to the faster growing pace of data volumes. Faced with this, subsampling is a popular technique to extract useful information from massive data. This paper focus on this technique and will develop optimal subsampling strategies for generalized linear models. 


Some commonly used generalized linear models include Poisson regression models for count data, logistic regression models for binary data, and binomial regression models for categorical data. To be precise, consider the  generalized linear model introduced in  \cite{Mccullagh1989Generalized}.  Let $(\bm x, y)$ be a pair of random variables with $\bm x\in \mathbb{R}^{p}$ being the covariate variable and $y \in \mathbb{R}$ being the response variable. 
 The conditional distribution of the response given $\bm x$ is $f(y|\bm\beta, \bm x)=h(y)\exp(y{\bm\beta}^{T}\bm x-\psi({\bm \beta}^{T}\bm x))$, where ${\bm \beta}^{T}\bm x$ satisfies that $\int h(y) \exp(y{\bm \beta}^{T}\bm  x)\mu(dy)<\infty $  with $\psi$ being the normalized function and $\mu(\cdot )$ being the dominating measure (Lebesgue measure for continuous distribution and counting measure discrete distribution).
 The unknown parameter vector, $\bm\beta$, is typically estimated with maximum likelihood techniques, namely, the maximum likelihood estimator (MLE) of $\bm\beta$ is the maximizer of the following log-likelihood function,  
 \begin{equation}\label{eq:likelihood}
   L(\beta)=\sum_{i=1}^n [y_i{\bm\beta}^{T}\bm x_i-\psi({\bm \beta}^{T}\bm x_i)].
 \end{equation}

 The MLE are often found numerically by using a Newton-Raphson method with updates of the form:
 \begin{equation}\label{eq:newton}
 {\bm{\beta}}^{(t+1)}={\bm {\beta }}^{(t)}+{\mathcal {J}}^{-1}({\bm {\beta }}^{(t)})\frac{\partial L({\bm {\beta }}^{(t)})}{\partial\bm \beta},
 \end{equation}
where ${{\mathcal {J}}({\bm {\beta }}^{(t)})}$  is the Hessian matrix (observed information matrix). However, fitting a generalized linear model on massive data is not an easy task through the iterative Newton-Raphson method, since it requires $O(np^2)$ time in each iteration of the optimization procedure.

When statistical methods are no longer applicable due to limited computing resources, a common practice to extract useful information from data is the subsampling method \citep[see][as an example]{drineas2006sampling}.
 This approach uses the estimator based on a subsample set that is taken randomly from the full data to approximate the full data estimator. 
 
 \cite{Drineas2011Faster} proposed to make a randomized Hadamard transform on data and then use uniform subsampling to take random subsamples to approximate ordinary least square estimators in linear regression models. \cite{Ma2015A} developed effective subsampling methods for linear regression models, which use normalized statistical leverage scores of the covariate matrix as non-uniform subsampling probabilities. \cite{Jia2014Influence} studied leverage sampling for generalized linear models based on generalized statistical leverage scores. 
 
 
 
 
 
 
 \cite{Wang2017Optimal} developed an optimal subsampling procedure to minimize the asymptotic  mean squared error (MSE) of the resultant subsample-estimator given the full data which is based on $A$- or $L$-optimality criteria in the language of optimal design. 


\cite{Wang2017Information} proposed a new algorithm called information-based optimal subdata selection method for linear regressions on big data.
The basic idea is to select the most informative data points deterministically based on $D$-optimality without relaying on random subsampling. 


To the best of our knowledge, theoretical and methodological investigations with statistical guarantees on subsampling methods for massive data are still limited. 


Most of the existing results concern linear regression models such as in \cite{Ma2015A} and  \cite{Wang2017Information}. The optimal subsampling method in \cite{Wang2017Optimal} is designed specifically for logistic regression models. In addition, they did not consider finite sample properties of the subsample-estimator. 

  
  
  In this paper, we fill these gaps by deriving optimal subsampling probabilities for generalized linear regression models. Furthermore, we will derive 
 finite-sample upper bounds for approximation errors that can be practically
used to make the balance between the subsample size and prediction accuracy. 
 This is an important aspect of subsampling method that \cite{Wang2017Optimal} didn't consider.


The rest of this paper is organized as follows. Section \ref{sec:model} introduces the model setup and derives asymptotic properties for the general subsampling estimator. 
 Section \ref{sec:appr-optim-subs} derives  optimal subsampling strategies based on $A$- and $L$-optimality criteria  for generalized linear models. Finite-sample error bounds are also derived in this section. 
Section \ref{sec:two-step} designs a two-step algorithm to approximate the optimal subsampling procedure and  
obtains asymptotic properties of the resultant estimator. 
Section \ref{sec:experiments} illustrates our methodology through numerical simulations. 
Real data applications are presented in Section \ref{sec:realdata}. Section \ref{sec:summary} concludes this paper with some discussions.  
All proofs are given in the Appendix.

\section{Preliminaries}\label{sec:model}

In this section, we present the generalized linear regression models and extend the general subsampling algorithm mentioned in \cite{Wang2017Optimal} to this framework.


\subsection{Models and Assumptions}\label{sec:glm}
Recall the definition of one parameter exponential family of distributions
\begin{equation}\label{eq:expfamily}
\left\lbrace f(y|\theta)=h(y)\exp(\theta y -\psi(\theta)), \theta \in \Theta \right\rbrace,
\end{equation}
as in (5.50) of \cite{Efron2016}. Here $f(\cdot|\theta)$ is a probability density function for the continuous case or a probability mass function for the discrete case; $h(\cdot)$ is a specific function that does not depend on $\theta$; and $\Theta$ is the parameter space defined as $\Theta:= \{ \theta \in \mathbb{R}:\: \int h(x) \exp(\theta x)\mu(dx)<\infty \}$ with $\mu$ being the dominating measure. In the formulation in (\ref{eq:expfamily}), $\theta$ is called the canonical parameter, and $\Theta$ is called the natural parameter space. 
 The exponential family includes most of the commonly used distributions such as normal, gamma, Poisson, and binomial distributions \citep[see][Chapter~8]{Efron2016}.

A key tactic for a generalized linear regression model is to express $\theta$ in form of a linear function of regression coefficients. 
Let $(\bm x,y)$ be a pair of random variables where $y \in\mathbb{R}$ and $\bm x\in \mathbb{R}^{p}$. A generalized linear regression model assumes that the conditional distribution of $y$ given $\bm x$ is determined by $\theta={\bm\beta}^{T}\bm x$. Specifically for exponential family, it assumes that the distribution of $y\vert \bm x$ is
\begin{equation}\label{eq:1}
  f(y\vert\bm\beta,\bm x)=h(y)\exp(y{\bm\beta} ^{T}\bm x-\psi({\bm\beta} ^{T}\bm x)),
  \quad\text{ with }\quad {\bm\beta}^{T}\bm x \in \Theta. 
\end{equation}
The problem of interest is to estimate the unknown $\bm\beta$ from observed data. 
 To facilitate the presentation, denote the full data matrix by $\mathcal{F}_n=(\bm X_d, \bm y_d)$, where
$\bm X_d=(\bm x_1, \dots,\bm x_n)^T$ 
is the covariate matrix and $\bm y_d=(y_1,y_2,\dots,y_n)^T$ consist of the corresponding responses. In this paper, we consider the case that $(\bm x_i, y_i)$'s are independent observations. 

\subsection{General Subsampling Algorithm and its Asymptotic Properties}\label{sec:subsamplingalg}
In this subsection, the general subsampling algorithm developed in \cite{Wang2017Optimal} is extended  to generalized linear models, and then  some asymptotic results for the resultant estimator are obtained.

Let $S$ be a set of subsample with $r$ data points, and define the sampling distribution $\pi_i$, for all data points $i=1,2,...n$ as $\bm \pi$. A general subsampling algorithm is presented in Algorithm \ref{alg:1}.
\begin{algorithm}[h]
  \begin{itemize}
  \item[ ] \textbf{Input:} data matrix $\mathcal{F}_n=(\bm X_d, \bm y_d)$, sampling distribution $\bm \pi$, and subsample size $r$.
  \item[ ]  \textbf{Output:} approximated estimator $\tilde{\bm\beta}$.
    
  \item[ ]  {Initialization:} ${i^*}=0,S=\varnothing$.
  \item[ ]  {For $t=1,\ldots,r$ do:}
  \begin{itemize}
    \item sample ${i^*}\sim\bm\pi$;
    \item update $S_{i^*}=S_{{i^*}-1}\cup\{(y_{i}^*,\bm {x}_{i}^*,\pi_{i}^*)\}$, where $\bm x^{*}_i$, $y^{*}_i$, and $\pi_i^{*}$ stand for covariates, responses, and subsampling probabilities in the subsample, respectively.
  \end{itemize}
  \item[ ] {Estimation:}
  maximize the following weighted log-likelihood function to obtain the estimate $\tilde{\bm\beta}$ based on the subsample set $S$,
    \begin{equation}\label{eq:reweigt2}
      L^*(\bm\beta)=\frac{1}{r}\sum^r_{t=1}\frac{1}{{\pi}_{i}^*}[y_{i}^*{\bm\beta}^{T}\bm x_{i}^*-\psi({\bm \beta}^{T}\bm x_{i}^*)].
    \end{equation}
    

  Finally, output the weighted MLE $\tilde{\bm\beta}$ derived from maximizing (\ref{eq:reweigt2}).
  \end{itemize}
  \caption{General subsampling algorithm}
  \label{alg:1}
\end{algorithm}
    The weighted MLE in Algorithm \ref{alg:1} can be implemented by Newton-Raphson method, i.e., iteratively applying the
    following formula until $\tilde{\bm\beta}^{(t+1)}$ and $\tilde{\bm\beta}^{(t)}$ are close enough,
    \begin{equation}\label{eq:21}
      \tilde{\bm\beta}^{(t+1)}=\tilde{\bm\beta}^{(t)}
      +\left\{\sum_{i=1}^r\frac{\ddot{\psi}(\bm x_i^{*T}\tilde{\bm\beta}^{(t)})
          \bm x_i^*(\bm x_i^*)^T}{{\pi}_i^*}\right\}^{-1}
      \sum^r_{i=1}\frac{\big\{y_i^*-\dot{\psi}(\bm x_i^{*T}\tilde{\bm\beta}^{(t)})\big\}
        \bm x_i^*}{{\pi}_i^*},
    \end{equation}
    where $\dot\psi(t)$ and $\ddot\psi(t)$ denote the first and the second derivatives of $\psi(t)$, respectively. 

An important feature of Algorithm \ref{alg:1} is that subsample estimator is essentially a weighted MLE and the corresponding weights are inverses of subsampling probabilities. 
 This is analogous to the Hansen-Hurwitz estimator \citep{Hansen1943On} in classic sampling techniques.
 For an overview see \cite{Sarndal1992Model}. 
 Although \cite{Ma2015A} showed that the unweighted subsample estimator is asymptotically unbiased for $\bm\beta$ in leveraging sampling, an unweighted subsample estimator is in general biased if the sampling distribution $\bm \pi$ depends on the responses. The inverse-probability weighting scheme is to remove bias, and we restrict our analysis on the weighted estimator here.

 In order to characterize asymptotic properties of subsampling estimators, we need some regularity assumptions listed below. 
 Note that we have no intent to make our assumptions the weakest possible. 

 \begin{itemize}
\item [\textbullet] (H.1): The variable $\bm x$ is almost surely bounded by a constant $L$, i.e., there exists a constant $L>0$ such that $\parallel\bm x\parallel \le L$ a.s., where $\|\cdot\|$ denote the euclidean norm.
\item [\textbullet] (H.2): For all $\bm x \in \left[-L,L \right]^{p} $, ${\bm\beta}^{T}\bm x \in \textrm{Int}(\Theta)$, where $\textrm{Int}(\Theta)$ denotes the interior of the parameter space $\Theta$.
\item [\textbullet] (H.3): The regression coefficient $\bm\beta$ lies in a compact domain:
  $$\Lambda _B=\left\lbrace  \bm\beta \in \mathbb{R}^{p}:\: \forall \bm x \in \left[-L,L \right]^{p},\bm\beta^{T}\bm x \in \Theta, \|\bm\beta\| \le B \right\rbrace$$
  for some large constant $B$. 


\item [\textbullet] (H.4): As $n\rightarrow\infty$, the observed information matrix
  $\mathcal{J}_X=n^{-1}\sum_{i=1}^n \ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\bm x_i\bm x_i^T$
  goes to a positive-definite matrix in probability. 
  
  
  
\item [\textbullet] (H.5): $ \max_{i=1,\ldots,n}(n\pi_i)^{-1}=O_P(1)$.
  
\end{itemize}

Here, assumptions (H.1)-(H.3) are the set of assumptions used in \cite{Blazere2014Oracle} and \cite{zhang2017elastic}. We listed them here for completeness. 
The bounded covariate assumption (H.1) is commonly used in literature for generalized linear models, see page 46 of \cite{fahrmeir2001multivariate}, and Section 2.2.7 of \cite{Chen2011Quasi} as examples. This assumption may be strong, and we relax it so that $\bm x$ can be unbounded in Section~\ref{sec:discuss-case-unbound} in the appendix. 


Assumption (H.3) is required for the proofs of consistency and asymptotic normality of subsampling estimators. 


Condition (H.4) is a typical assumption in the literature of likelihood theory and it is also assumed in \cite{Wang2017Optimal}. 
Assumption (H.5) imposes a restriction on the subsampling distribution. 
It is to protect the likelihood function from being inflated by data points with extremely small subsampling probabilities.

 The theorem below presents the consistency of the estimator from the subsampling algorithm to the full data MLE.
 \begin{theorem}\label{thm:as-general-alg}
  If Assumptions (H.1)--(H.5) hold, then as $n\rightarrow\infty$ and $r\rightarrow\infty$, $\tilde{\bm\beta}$ is consistent to $\hat{\bm\beta}_{\rm MLE}$ in conditional probability
  given $\mathcal{F}_n$. Moreover, the rate of convergence is $r^{-1/2}$.
  That is, with probability approaching one, for any $\epsilon>0$, there exist  finite $\Delta_\epsilon$ and $r_\epsilon$ such that
  \begin{equation}\label{eq:18}
    P(\|\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}\|\ge
    r^{-1/2}\Delta_\epsilon|\mathcal{F}_n)<\epsilon
  \end{equation}
  for all $r>r_\epsilon$.
\end{theorem}

\begin{remark}
 Note that if a sequence of random variables is bounded in conditional probability, then it is also bounded in unconditional probability \citep{xiong2008some}. Thus, \eqref{eq:18} implies that $\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}=O_P(r^{-1/2})$, meaning that subsampling estimator is $\sqrt{r}$-consistent to the full data estimator under the unconditional distribution as well.
\end{remark}

Besides consistency, we derive the asymptotic distribution of the approximation error, and prove that the approximation error,
$\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}$, is asymptotically normal. 


 
\begin{theorem}\label{thm:CLT}
  If Assumptions (H.1)--(H.5) hold, then as $n\rightarrow\infty$ and $r\rightarrow\infty$, conditional on $\mathcal{F}_n$ in probability,
  \begin{equation}\label{normal}
    V^{-1/2}(\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE})
    \longrightarrow N(0,I)
  \end{equation}
  in distribution, where
  \begin{equation}\label{varorder}
    V=\mathcal{J}_X^{-1}V_c\mathcal{J}_X^{-1}=O_{p}(r^{-1}),
  \end{equation}
  and
  \begin{equation}
    V_c=\frac{1}{rn^2}
    \sum_{i=1}^n\frac{\{y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^2\bm x_i\bm x_i^T}{\pi_i}.
    \label{eq:vc}
  \end{equation}
\end{theorem}


\begin{remark}
Result~\eqref{normal} shows that the distribution of
$\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}$ given $\mathcal{F}_n$ can be approximated by a normal random variable with distribution $N({\bf 0},V)$. For convenience, we denote the corresponding random variable as $\bm u$. It worth to mention that this does not necessarily imply that
$E(\|\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}\|^2|\mathcal{F}_n)$ is close to $E(\|\bm u\|^2|\mathcal{F}_n)$ unless $\|\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}\|^2$ is uniformly integrable.
\end{remark}

\section{Optimal Subsampling Strategies}\label{sec:appr-optim-subs}

In this section, we will consider how to specify subsampling distribution $\bm\pi=\{\pi_i\}_{i=1}^{n}$. Obviously, the uniform subsampling strategy, $\bm \pi^{\mathrm{UNI}}=\{\pi_i=n^{-1}\}_{i=1}^{n}$, satisfies Assumption (H.5) in pervious section and the resultant estimator will possess the asymptotic properties described in Theorems \ref{thm:as-general-alg} and \ref{thm:CLT}.
  However, Algorithm \ref{alg:1} with the uniform subsampling probabilities may not be optimal in terms of parameter approximation and a subsampling with unequal probabilities may have better performance. 
  To better approximate $\hat{\bm\beta}_{\rm MLE}$, our goal is to minimize the asymptotic MSE or its variant of $\tilde{\bm\beta}$ in approximating $\hat{\bm\beta}_{\rm MLE}$ 

with the same size of the subsamples. 
This is to adopt the same idea as the $A$-optimality criterion in the theory of design of experiments 

\citep[see][]{pukelsheim2006optimal}. 

\subsection{Minimum Asymptotic MSE of $\tilde{\bm\beta}$}\label{sec:mV}

From Theorem~\ref{thm:CLT}, the asymptotic MSE of $\tilde{\bm\beta}$ is equal to the trace of $V$, namely,
\begin{equation} \label{eq:A-opt}
  \mathrm{AMSE}(\tilde{\bm\beta})=E(\|\bm u\|^2|\mathcal{F}_n)=\tr(V),
\end{equation}
where $\rm tr(\cdot)$ denotes the trace function.

From \eqref{varorder}, $V$ depends on $\{\pi_i\}_{i=1}^n$, and clearly, $\bm \pi^{\mathrm{UNI}}$ may not produce the smallest value of $\rm tr(V)$.
 The key idea of optimal subsampling in \cite{Wang2017Optimal} is to choose subsampling probabilities such that the $\mathrm{AMSE}(\tilde{\bm\beta})$ in \eqref{eq:A-opt} is minimized and we can the resulting subsampling strategy mV-optimal.


\begin{theorem}\label{thm:3}
  In Algorithm 1, if the subsampling probability is chosen such that
  \begin{equation}\label{eq:pi-amse-w}
    \pi_i^{\mathrm{mV}}=
    \frac{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm
      x_i)|\|\mathcal{J}_X^{-1}\bm x_i\|}
    {\sum_{j=1}^n|y_j-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm
      x_j)|\|\mathcal{J}_X^{-1}\bm x_j\|},\;
    i=1,2,...,n,
  \end{equation}
  then  $\tr(V)$ attains its minimum.
\end{theorem}


As observed in (\ref{eq:pi-amse-w}), the optimal subsampling probability $\bm{\pi}^{\mathrm{mV}}=\{\pi_{i}^{\mathrm{mV}}\}_{i=1}^{n}$ depends on data through both the covariates and the responses directly.
 For the covariates, subsampling probabilities are higher for larger $\|\mathcal{J}_X^{-1}\bm x_i\|$, which are square roots of diagonal elements of the matrix  ${\bm X_d}\mathcal{J}_X^{-2}{\bm X_d}^T$.
 For linear regression, diagonal elements of ${\bm X_d}\mathcal{J}_X^{-1}{\bm X_d}^T$ are statistical leverage scores, usually denoted by $h_{ii}$ for $i=1, \ldots, n$.

Let $\sigma_{\max}(A)$ and $\sigma_{\min}(A)$ denote the maximum and minimum non-zero singular values of matrix $A$, respectively. The condition number of $A$ is $\kappa(A):=\sigma_{\max}(A)/\sigma_{\min}(A)$. 
 Note that $\sigma_{\min}(\mathcal{J}_{X}^{-1/2})h_{ii}\le \|\mathcal{J}_X^{-1}\bm x_i\| \le \sigma_{\max}(\mathcal{J}_{X}^{-1/2})h_{ii}$. 
 
 It can be seen that $\|\mathcal{J}_X^{-1}\bm x_i\|$'s measure effects of individual observations in terms of covariates in a way that is similar to statistical leverage scores. In other words, the optimal subsampling probability $\bm{\pi}^{\mathrm{mV}}$ give higher preference to high-leverage observations, which share the same property as leveraging subsampling. 
 
  The direct effect of the response on the optimal subsampling probability is through the term $|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|$. This term measures how surprise each data points are observed. Note that $\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)$ is the estimated mean of $y_i$ given $\bm x_i$ using the full data MLE. Thus, data points with large values of $|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|$ are harder to observe as they are in the tail region of the conditional distribution of the response. The optimal subsampling probability prefers to include more surprising data points in the subsample. This also improves the robustness of the subsample estimator. For example, if $y_i$ follows the Bernoulli distribution, i.e., the support of $y_i$ is \{0, 1\}, then the MLE uniquely exists if and only if the convex cone generated by $\bm x_i$ in the sample with $y_i=0$ overlaps with the convex cone generated by $\bm x_i$ in the sample with $y_i=1$ \citep{silvapulle1981existence}. 

 
 
 
 

  The optimal subsampling probability $\bm{\pi}^{\mathrm{mV}}$ also has a meaningful interpretation from the view-point of optimal design of experiments \citep{pukelsheim2006optimal}.
  
  
  
 
 


    Note that under mild condition the  ``empirical information matrix'' ${\mathcal{J}^e_X}=\frac{1}{n}\sum_{i=1}^n(y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i))^2\bm x_i\bm x_i^T$ and $\mathcal{J}_X$ converge to the same limit, the Fisher information matrix of model~\eqref{eq:1}. This means that ${\mathcal{J}^e_X}-\mathcal{J}_X=o_P(1)$. Thus, $\mathcal{J}_X$ can be replaced by ${\mathcal{J}^e_X}$ in $\bm{\pi}^{\mathrm{mV}}$, because Theorem~\ref{thm:CLT} still holds if $\mathcal{J}_X$ is replaced by ${\mathcal{J}^e_X}$ in \eqref{varorder}. 
    Let $\eta_{\bm x_i}=(y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i))^2\bm x_i\bm x_i^T$, the contribution of the $i$-th observation to the empirical information matrix, and ${\mathcal{J}^e_{X\bm x_i\alpha}}=(1-\alpha){\mathcal{J}^e_X}+\alpha\eta_{\bm x_i}$, which can be interpreted as a movement of the information matrix in a direction determined by the $i$-th observation. The directional derivative of $\text{tr}({\mathcal{J}^e_X}^{-1})$ through the direction determined by the $i$th observation is $F_i=\lim_{\alpha\rightarrow0+} \alpha^{-1}\{\text{tr}({\mathcal{J}^e_X}^{-1})-\text{tr}({\mathcal{J}^e_{X\bm x_i\alpha}}^{-1})\}$. This directional derivative is used to measure the relative gain in estimation efficiency under the $A$-optimality by adding the $i$-th observations into the sample. 


Thus, the optimal subsampling strategy prefers to select the data points with large values of directional derivatives, i.e., data points that will result in a larger gain under the $A$-optimality. 

\subsection{Minimum Asymptotic MSE of $\mathcal{J}_X\tilde{\bm\beta}$}\label{sec:mVc}

The optimal subsampling strategy  derived in the previous section requires the calculation of $\|\mathcal{J}_X^{-1}\bm x_i\|$ for $i=1,2,...,n$, which takes $O(np^2)$ time.
 To reduce the calculating time, \cite{Wang2017Optimal}  proposed to minimize a modified optimality criterion, i.e., minimize $\tr(V_c)$, instead.
 This criterion essentially is the linear optimality (L-optimality) criterion in optimal experimental design \citep[see][]{Atkinson2007Optimum}, which is to improve the quality of some linear combinations for unknown parameters.

 
Here are additional rationale for this optimality. Note that $\mathcal{J}_X$ and $V_c$ are nonnegative definite, and $V=\mathcal{J}_X^{-1}V_c\mathcal{J}_X^{-1}$. Simple matrix algebra yields that ${\rm{tr}}(V) = {\rm{tr}}({V_c}{\cal J}_X^{ - 2}) \le {\sigma _{\max }}({\cal J}_X^{ - 2}){\rm{tr}}({V_c})$. 
Since $\sigma _{\max }({\cal J}_X^{ - 2})$ does not depend on $\bm{\pi}$, minimizing $\tr(V_c)$ minimizes an upper bound of $\tr(V)$. 
 In fact $V$ depends on $\bm{\pi}$ through $V_c$, and $\mathcal{J}_X$ does not depend on $\bm{\pi}$.
  For two given subsampling strategies  $\bm{\pi}^{(1)}$ and $\bm{\pi}^{(2)}$, $V(\bm{\pi}^{(1)})\le V(\bm{\pi}^{(2)})$ in the sense of Loewner-ordering if and only if $V_c(\bm{\pi}^{(1)})\le V_c(\bm{\pi}^{(2)})$.
 This alternative  optimality criterion is to further reduce the computing time without scarifying the estimation efficiency too much.


The following theorem gives the optimal subsampling probability that minimizes the trace of $V_c$.
\begin{theorem}\label{thm:5}
  In Algorithm 1, if the subsampling probability is chosen such that
  \begin{equation}\label{eq:optimalBpi-w}
    \pi_i^{\mathrm{mVc}}=\frac{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\|\bm x_i\|}
    {\sum_{j=1}^n|y_j-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_j)|\|\bm x_j\|},\ i=1,2, ..., n,
  \end{equation}
  then $\tr(V_c)$ attains its minimum.
\end{theorem}


It turns out that the alternative optimality criterion greatly reduces the computing time without scarifiing too much accuracy. From Theorem~\ref{thm:5}, the effect of the covariates on
$\bm{\pi}^{\mathrm{mVc}}=\{\pi_i^{\mathrm{mVc}}\}_{i=1}^{n}$ is presented by $\|\bm x_i\|$, instead of $\|\mathcal{J}_X^{-1}\bm x_i\|$ as in $\bm{\pi}^{\mathrm{mV}}$. The computational benefit is obvious: it needs $O(nd)$ time to evaluate $\|\bm x_i\|$ for $i=1,2,...,n$, which is significantly less than the required $O(nd^2)$ time to calculate
$\|\mathcal{J}_X^{-1}\bm x_i\|$ for $i=1,2,...,n$. 

Note that the score function for the log-likelihood function is
\[U(\bm\beta)=\sum_{i=1}^n \{ y_i - \dot\psi(\bm\beta^T\bm x_i)\}\bm x_i.\]
It is interesting that the $\pi_i^{\mathrm{mVc}}$'s in Theorem \ref{thm:5} are proportional to $\|\{y_i - \dot\psi(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}\bm x_i\|$, norms of gradients of the log-likelihood at individual data points evaluated at the full data MLE. This is trying to find the subsample that best approximate the full data score function at the full data MLE. 


\subsection{Non-asymptotic Properties 
}\label{sec:r-size}


 
 


 


 
 
 
 
 


 
  


 
 
 


  


 
 
 
 

 
 
 
 
 

 


In this section, we derive some finite sample properties of the subsample estimator based on optimal subsampling probabilities $\bm{\pi}^{\mathrm{mV}}$ and $\bm{\pi}^{\mathrm{mVc}}$. The results are presented in the form of excess risk for approximating the mean responses and they hold for fixed $r$ and $n$ without requiring any quantity to go to infinity. These results shows the factors that affect the approximation accuracy.  


Since $\dot{\psi}(\bm x_i^T\bm\beta)$ is the conditional expectation of the response $y_i$ given $\bm x_i$, we aim to characterize the quantity of $\tilde{\bm\beta}$ in prediction by examining $\|\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})-\dot{\psi}(\bm X_d^T\tilde{\bm\beta})\|$. This quantity is the distance between the estimated conditional mean responses based on the full data and that based on the subsample. Intuitively, it measures the goodness of fit in using subsample estimator to predict the mean responses.
Note that we can always improve the accuracy of the estimator by increasing the subsample size $r$. Here we want to have a closer look at the effects of different quantities such as covariate matrix and data dimension and the effect of subsample size $r$ on approximation accuracy.

 
 
 
 
 


 Recall that $\sigma_{\max}(A)$ and $\sigma_{\min}(A)$ are the maximum and minimum non-zero singular values of matrix $A$, respectively, and  $\kappa(A):=\sigma_{\max}(A)/\sigma_{\min}(A)$. 
  For the estimator $\tilde{\bm \beta}$ obtained from Algorithm \ref{alg:1} based on the subsampling probabilities, $\bm{\pi}^{\mathrm{mV}}$ and $\bm{\pi}^{\mathrm{mVc}}$,
  the following theorem holds.
 \begin{theorem}\label{thm:r-size}
 Let $\bm{\tilde X}_d$ denotes the design matrix consisted by the subsamples and each sampled element has been rescaled by $1/\sqrt{r\pi_i^*}$.  Assume that both $\sigma_{\max}(\bm X_d)/\sqrt{n}$ and $\sigma_{\min}(\bm X_d)/\sqrt{n}$ are bounded 
 
  and $\sigma_{\min}^2(\bm{\tilde{ X}_d})\ge 0.5\sigma^2_{\min}(\bm X_d)$. For any given $\epsilon\in(0,1/3)$,   with probability at least $1-\epsilon$:
 \begin{equation}\label{eq:r-size}
   \|\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})-\dot{\psi}(\bm X_d^T\tilde{\bm\beta})\|\le \frac{8\alpha\kappa^2(\bm X_d)\sqrt{p\log(1/\epsilon)}}{\sqrt{r}}\|\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})\|,
 \end{equation}
 where 
 $\alpha=\kappa(\mathcal{J}_X^{-1})$ for $\bm{\pi}^{\mathrm{mV}}$ and $\alpha=1$ for $\bm{\pi}^{\mathrm{mVc}}$.
 \end{theorem}

 Theorem \ref{thm:r-size} not only indicates that the accuracy increase with subsample size $r$, which agrees with the results in Theorem \ref{thm:as-general-alg}, but also enables us to have a closer look at the effects of different quantities such as covariate matrix and data dimension and the effect of subsample size $r$ on approximation accuracy.
 
 Heuristically, the condition number of $\bm X_d$ measures the collinearity of covariates in the full data covariate matrix; $p$ shows the curse of dimensionality; and $\|\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})\|$ measures the goodness of fit of the underlying model on the full data.


The result in (\ref{eq:r-size}) also indicates that we should choose $r\propto p$ to control the error bound, hence it seems reasonable to choose the subsample size as $r=cp$. This agrees with the 

recommendation of choosing a sample size as large as 10 times of the number of covariates in \cite{Chapman1994Arctic} 
and \cite{Loeppky2009Choosing} for designed experiments.  
However, for designed experiments, covariate matrices are often orthogonal or close to be orthogonal, so $\kappa(\bm X_d)$ is equal or close to 1 in these cases. For the current paper, full data may not obtained from well designed experiments so $\bm X_d$ may vary a lot. Thus, it should also be considered in determining required subsample size for a given level of prediction accuracy.  
  


 

 The particular constant $0.5$ in condition $\sigma_{\min}^2(\bm{\tilde{ X}_d})\ge 0.5\sigma^2_{\min}(\bm X_d)$ of Theorem \ref{thm:r-size} { can be replaced by any constant between 0 and 1. Here we follow the setting of \cite{Drineas2011Faster} and choose 0.5 for convenience. 
 }
 This condition indicates that the rank of $\bm{\tilde{X}_d}$ is the same as that of $\bm X_d$. More details and interpretations about this condition can be found in \cite{Mahoney2012Randomized}.
 Using similar argument as in the proof of Theorem \ref{thm:r-size}, it is proved that this condition holds with high probability.

 \begin{theorem}\label{thm:r-size1}
 Let $\bm{\tilde X}_d$ denote the design matrix consisted by the subsamples and each sampled element has been rescaled by $1/\sqrt{r\pi_i^*}$.
 For any given $\epsilon\in(0,1/3)$, let $r>64{c_d}\alpha^2\log(1/\epsilon)p^2/\kappa^4(\bm X_d)$ where ${c_d} \le 1$ is a constant depending on full data matrix $\bm{X}_d$, and assume that $\sigma_{\max}(\bm X_d)$, $\sigma_{\min}(\bm X_d)$ are bounded, then with probability at least $1-\epsilon$:
 \[\sigma_{\min}^2(\bm{\tilde{ X}_d})\ge 0.5\sigma^2_{\min}(\bm X_d),\]
 where $\alpha=\kappa(\mathcal{J}_X^{-1})$ for $\bm{\pi}^{\mathrm{mV}}$ and $\alpha=1$ for $\bm{\pi}^{\mathrm{mVc}}$.
 \end{theorem}


 \section{Practical Consideration and Implementation}\label{sec:two-step}

 
 
  
 
  
  
 


 For practical implementation, the optimal subsampling probabilities $\pi_i^{\mathrm{mV}}$'s and $\pi_i^{\mathrm{mVc}}$'s cannot be used directly because they depend on the unknown full data MLE $\hat{\bm\beta}_{\rm MLE}$. 
  
 
 As suggested in \cite{Wang2017Optimal}, in order to calculate $\bm \pi^{\mathrm{mV}}$ or $\bm\pi^{\mathrm{mVc}}$, a pilot estimator of $\hat{\bm\beta}_{\rm MLE}$ has to be used. Let $\tilde{\bm\beta}_0$ be a pilot estimator based on a subsample of size $r_0$. It can be used to replace $\hat{\bm\beta}_{\rm MLE}$ in $\bm\pi^{\mathrm{mV}}$ or $\bm\pi^{\mathrm{mVc}}$, which then can be used to sample more informative subsamples.
 
 From the expression of $\bm\pi^{\mathrm{mV}}$ or $\bm\pi^{\mathrm{mVc}}$, the approximated optimal subsampling probabilities are both proportional to $|y_i-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)|$, so data points with $y_i\approx\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)$ have very small probabilities to be selected and data points with  $y_i=\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)$ will never be included in a subsample. 
 On the other hand, if these data points are included in the subsample, the weighted log-likelihood function in (\ref{eq:reweigt2}) may be dominated by them. As a result, the subsample estimator may be sensitive to these data points. \cite{Ma2015A} also noticed the problem that some extremely small subsampling probabilities may inflate the variance of the subsampling estimator in the context of leveraging sampling.

 
 To protect the weighted log-likelihood function  from being inflated by these data points in practical implementation, we propose to set a threshold, say $\delta$, for $|y_i-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)|$, i.e., use $\max\{|y_i-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)|,\delta\}$ to replace $|y_i-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)|$. Here, $\delta$ is a small positive number, say $10^{-6}$ as an example. 
 Setting a threshold $\delta$ in subsampling probabilities results in a truncation in the weights for the subsample weighted log-likelihood. 
 Truncating the weight function is a commonly used technique in practice for 
 
 robust estimation. 
 Note that in practical application, an intercept should always be included in a model, so it is typical that $\|\bm x_i\|$ and $\|\mathcal{J}_X^{-1}\bm x_i\|$ are bounded away from 0 and we do not need to set a threshold for them. 


 
 With threshold $\delta$ and pilot estimate $\tilde{\bm\beta}_0$, the practical subsampling probabilities corresponding to $\bm\pi^{\mathrm{mV}}$ and $\bm\pi^{\mathrm{mVc}}$ are 
\begin{align}
    \tilde\pi_i^{\mathrm{mV}}&=
    \frac{\max(|y_i-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)|,\delta)\|\tilde{\mathcal{J}}_X^{-1}\bm x_i\|}
    {\sum_{j=1}^n\max(|y_j-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_j)|,\delta)\|\tilde{\mathcal{J}}_X^{-1}\bm x_j\|},\label{eq:2}\\
    \tilde\pi_i^{\mathrm{mVc}}&=
    \frac{\max(|y_i-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)|,\delta)\|\bm x_i\|}
    {\sum_{j=1}^n\max(|y_j-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_j)|,\delta)\|\bm x_j\|},\label{eq:3}
\end{align}
respectively. They are minimizers of $\tr(\tilde{V}^\delta)$ and $\tr(\tilde{V}_c^\delta)$, respectively, where  $\tilde{V}^\delta=\tilde{\mathcal{J}}_X^{-1}\tilde{V}_c^\delta\tilde{\mathcal{J}}_X^{-1}$,


\[ \tilde{V}_c^\delta=\frac{1}{rn^2}
  \sum_{i=1}^n\frac{\max(\{y_i-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)\}^2,\delta^2)\bm x_i\bm x_i^T}{\pi_i},\]
and $\tilde{\mathcal{J}}_X=(nr_0)^{-1}\sum_{i\in \tilde{S}_{r_0}}(\pi_i^*)^{-1}\ddot{\psi}(\tilde{\bm\beta}_0^T\bm x_i^*)\bm x_i^*\bm x_i^{*T}$. Here, we can also approximate $\mathcal{J}_X$ by using $\tilde{\bm\beta}_0$ to replace $\hat{\bm\beta}_{\rm MLE}$ in its expression. However, this approach requires an additional round of calculation on the full data. 


Let $\tilde{V}$ be the version of $V$ with $\hat{\bm\beta}_{\rm MLE}$ substituted by $\tilde{\bm\beta}_0$. It can be shown that
  
\[\tr(\tilde{V})\le \tr(\tilde{V}^\delta)\le \tr(\tilde{V})+\frac{\delta^2}{n^2r}\sum_{i=1}^n\frac{1}{\pi_i}\|\tilde{\mathcal{J}}_X^{-1}\bm x_i\|^2. \]
Thus, minimizing $\tr(\tilde{V}^\delta)$ is close to minimizing $\tr(\tilde{V})$ if $\delta$ is sufficiently small. The threshold $\delta$ is 

to make our subsampling estimator more robust without scarifying the estimation efficiency too much. 


  


 
 For transparent presentation, we combine all the aforementioned practical considerations in this section and present a two-step algorithm in Algorithm~\ref{alg:3}.

\begin{algorithm}[htp]
  \caption{Two-step Algorithm}\label{alg:3}
  \begin{itemize}
  \item[ ] \textbf{Input:} covariate matrix $\bm X_d\in \mathbb{R}^{n\times p}$, response vector $\bm y_d\in \mathbb{R}^p$, pilot subsampling distribution $\bm\pi_0$, pilot subsample size $r_0$, subsample size $r$, and threshold $\delta$.
  \item[ ] \textbf{Output:} approximated estimator ${\breve{\bm\beta}}$ based on the total subsample of size $r_0+r$.
  \item[ ]  {Initialization:} $i^*=0,S=\varnothing$.
  \item[ ]  Run Algorithm~\ref{alg:1} with subsample size $r_0$ and pilot sampling distribution $\bm\pi_0$ to obtain a subsample set $\tilde{S}_{r_0}$ and a pilot estimate $\tilde{\bm\beta}_0$. 
  \item[ ] Using $\tilde{\bm\beta}_0$ to calculate approximated subsampling probabilities $\tilde{\bm\pi}^{\rm opt}=\{\tilde\pi_i^{\mathrm{mV}}\}_{i=1}^n$ or $\tilde{\bm\pi}^{\rm opt}=\{\tilde\pi_i^{\mathrm{mVc}}\}_{i=1}^n$, where $\tilde\pi_i^{\mathrm{mV}}$ and $\tilde\pi_i^{\mathrm{mVc}}$ are defined in \eqref{eq:2} and \eqref{eq:3}, respectively.
    
    
  \item[ ]  {For $i^*=1,\ldots,r$ do:}
  \begin{itemize}
    \item sample $i^* \sim\tilde{\bm\pi}^{\rm opt}$;
    \item update $S_{i^*}=S_{{i^*}-1}\cup\{(y_{i}^*,\bm {x}_{i}^*,\pi_{i}^*)\}$ with $S_0=\tilde{S}_{r_0}$.
  \end{itemize}
  \item[ ] {Estimation:}
  maximize the following weighted log-likelihood function 
  (\ref{eq:reweigt1}) to obtain the 
  
  estimator ${\breve{\bm\beta}}$
  \begin{equation}\label{eq:reweigt1}
      L^*(\bm\beta)=\frac{1}{r+r_0}\sum_{i\in S_{r^*}}\frac{1}{{\pi}_{i}^*}[y_{i}^*{\bm\beta}^{T}\bm x_{i}^*-\psi({\bm \beta}^{T}\bm x_{i}^*)].
  \end{equation}
  \end{itemize}
\end{algorithm}

The pilot sampling distribution $\bm\pi_0$ in Algorithm~\ref{alg:3} can be simply selected as the discrete uniform distribution, i.e., $\pi_{i0}=n^{-1}$ for $i=1,\ldots,n$ or some other informative sampling distribution satisfying condition (H.5).

As shown in Theorem \ref{thm:as-general-alg}, the pilot $\tilde{\bm\beta}_0$ approximates $\hat{\bm\beta}_{\rm MLE}$ accurately as long as $r_0$ is not too small.
 On the other hand, the efficiency of the two-step algorithm would decrease, if $r_0$ gets close to the total subsample size $r_0+r$ and $r$ is relatively small. 
   Empirically from simulations in various settings, we find that the
   algorithm works well when the proportion $r_0/r$ is around 0.2.


  Comparing with the full data approach, Algorithm~\ref{alg:3} greatly reduces the computational cost.
  The major computing time is to approximate the optimal subsampling probabilities which does not require iterative calculations on the full data. Once the approximately optimal subsampling probabilities are available, the time to obtain $\breve{\bm\beta}$ in the second step is $O\{\zeta(r_0+r)d^2\}$ where ${\zeta}$ is the number of iterations of the iterative procedure in optimizing the log-likelihood in {\red (\ref{eq:reweigt1})}. For the estimator $\breve{\bm\beta}$ obtained from Algorithm~\ref{alg:3} based on the subsampling probabilities $\tilde{\pi}^{\mathrm{mVc}}_i$'s or $\tilde{\pi}^{\mathrm{mV}}_i$'s, we establish its asymptotic properties as follows. 
  
\begin{theorem}\label{thm:asy-2step-alg}
Under Assumptions~(H.1)--(H.4), if 
$r_0r^{-1}\rightarrow0$, then for the estimator $\breve{\bm\beta}$ obtained from Algorithm~\ref{alg:3},  
  as $r_0\rightarrow\infty, {r}\rightarrow\infty$ and $n\rightarrow\infty$, with probability approaching one, for any $\epsilon>0$, there exist   finite $\Delta_\epsilon$ and $r_\epsilon$ such that
  \begin{equation*}
  P(\|{\breve{\bm\beta}}-\hat{\bm\beta}_{\rm MLE}\|\ge {r}^{-1/2}\Delta_\epsilon|\mathcal{F}_n)<\epsilon
\end{equation*}
for all ${r}>r_\epsilon$.


\end{theorem}

The result of asymptotic normality is presented in the following theorem.
\begin{theorem}\label{thm:CLT-2step}
Under assumptions (H.1)--(H.4), if $r_0r^{-1}\rightarrow0$, then as $r_0\rightarrow\infty$, $r\rightarrow\infty$ and $n\rightarrow\infty$, conditional on $\mathcal{F}_n$,
  \begin{equation*}
    V_{opt}^{-1/2}(\breve{\bm\beta}-\hat{\bm\beta}_{\rm MLE}) \rightarrow N(0,I),
  \end{equation*}
 where $V_{opt}=\mathcal{J}_X^{-1}V_{c,opt}\mathcal{J}_X^{-1}$,
  \begin{equation*}
    V_{c,opt}=\frac{1}{r}\onen
    \sum_{i=1}^n \frac{\{y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^2\bm x_i\bm x_i^T}{ {\max(|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|,{\delta})\|\bm x_i\|}}
    \times\onen\sumn\max(|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|,{\delta})\|\bm x_i\|,
  \end{equation*}
  for the estimator obtained from Algorithm~\ref{alg:3} based on the subsampling probabilities $\tilde{\pi}^{\mathrm{mVc}}_i$'s, and
  \begin{equation*}
    V_{c,opt}=\frac{1}{r}\onen
    \sumn\frac{\{y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^2\bm x_i\bm x_i^T}{{\max(|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|,{\delta})\|\mathcal{J}_X^{-1}\bm x_i\|}}
    \times\onen\sumn\max(|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|,{\delta})\|\mathcal{J}_X^{-1}\bm x_i\|,
  \end{equation*}
  
  
  
  
  for the estimator obtained from Algorithm~\ref{alg:3} based on the subsampling probabilities $\tilde{\pi}^{\mathrm{mV}}_i$'s.
  
 

\end{theorem}


In order to get standard error of the corresponding estimator, we adopt the method of moments to estimate the variance-covariance matrix of $\breve{\bm\beta}$. To be precise, the variance-covariance matrix of $\breve{\bm\beta}$ is estimated by 
 \begin{equation}\label{eq:24}
   \breve{V}=\breve{\mathcal{J}}_X^{-1}\breve{V}_c\breve{\mathcal{J}}_X^{-1},
 \end{equation}
where
\begin{equation*}
\breve{\mathcal{J}}_X=\frac{1}{n(r_0+r)}\left\{\sum_{i=1}^{r_0}\frac{\ddot{\psi}(\breve{\bm\beta}^T\bm x_i^*)\bm x_i^*(\bm x_i)^{*T}}{\pi_{i0}^*}+\sum_{i=1}^{r} \frac{\ddot{\psi}(\breve{\bm\beta}^T\bm x_i^*)\bm x_i^*(\bm x_i)^{*T}}{\tilde\pi_i^*}\right\},
\end{equation*}
  \begin{equation*}
    \breve{V}_c=\frac{1}{n^2(r_0+r)^2}\left\{\sum_{i=1}^{r_0}\frac{\{y_i-\dot{\psi}(\breve{\bm\beta}^T\bm x_i^*)\}^2\bm x_i^*(\bm x_i^*)^T}{(\tilde\pi_{i0}^*)^2}
    +\sum_{i=1}^{r}\frac{\{y_i-\dot{\psi}(\breve{\bm\beta}^T\bm x_i^*)\}^2\bm x_i^*(\bm x_i^*)^T}{(\tilde\pi_i^*)^2}\right\},
  \end{equation*}
 and $\tilde\pi_i^*=\tilde\pi_i^{\mathrm{mV}*}$ or $\tilde\pi_i^{\mathrm{mVc}*}$ for $i=1,\ldots,r$.
 


\section{Numerical Studies}\label{sec:experiments}

\subsection{Simulation Studies}


In this section, we use simulation to evaluate the finite sample performance of 
the proposed method in Poisson regression. 
Full data of size $n=100,000$ are generated from the model $y \sim \mathcal{P}(\exp(\bm\beta^{T}\bm x))$
with the true value of $\bm\beta$, $\bm\beta_0$, being a $7\times1$ vector
of 0.5. We consider the following four cases to generate the covariates $\bm x_i=(x_{i1}, ..., x_{i7})^{T}$.
\begin{enumerate}[{Case} 1:]
\item The seven covariates are independent and identically distributed (i.i.d) from the standard uniform distribution, namely, $x_{ij}\overset{\text{i.i.d}}{\sim} U(\left[0,1 \right])$ for $j=1, ..., 7$.
\item The first two covariates are highly correlated. Specifically, $x_{ij}\overset{\text{i.i.d}}{\sim} U(\left[0,1 \right])$ for all $j$ except for $x_{i2}=x_{i1}+\varepsilon_{i}$ with $\varepsilon_{i}\overset{\text{i.i.d}}{\sim} U(\left[0,0.1 \right] )$. For this setup, the correlation coefficient between the first two covariates are about 0.8.
 
 
 
 
\item
This case is the same as the second one except that $\varepsilon_{i}$ $\overset{\text{i.i.d}}{\sim} U(\left[0,1 \right] )$. For this case, the correlation between the first two covariates is close to $0.5$.

\item This case is the same as the third one except that  $x_{ij}\overset{\text{i.i.d}}{\sim} U(\left[-1, 1 \right] )$ for $j=6, 7$. For this case, the bounds for each covariates are not all the same.
\end{enumerate}

\iffalse
\begin{example}\label{ex1}
  Consider the ideal case that $\hat{\bm\beta}_{\rm MLE}$ is used to calculate the subsampling probabilities. Note that this case is not practical because if $\hat{\bm\beta}_{\rm MLE}$ is available, there is no need to perform subsampling anymore. \red Let's talk about this.
  We mainly consider five scenarios here.
 The first scenario is that subsampling method is taken as ``uniform" sampling, i.e., $\pi_i=1/n$ for all $i$.
 The second is that $\pi_i=\pi_i^{\mathrm{mV}}$ with $\delta=10^{-6}$.
 Third is that $\pi_i=\pi_i^{\mathrm{mVc}}$ with $\delta=10^{-6}$.
 The last two scenarios are the same as the third and fourth except $\delta=1$.
 $\pi_i^{\mathrm{mV}}$ and $\pi_i^{\mathrm{mVc}}$ is calculated as discussed in Section \ref{sec:appr-optim-subs}.

 For comparison, we also provide the results for the leverage score subsampling strategy considered in \cite{Ma2015A}, see also \cite{Jia2014Influence} for the generalized linear models case. Each row of the design matrix $\bm X_d$, $\bm x_i$, is sampled  according to an importance sampling distribution that is proportional to the statistical leverage scores $\bm x_i(\bm X_d^T\bm X_d)^{-1}\bm x_i$.
 It worth to mention that this strategy is designed for subsampling in linear regression models and for generalized linear models, the leverage score is adjusted by replacing $\bm x_i,\bm X_d$ with $\tilde{\bm x_i}$ and $\tilde{\bm X_d}$, where $\tilde{\bm x_i}$ is defined as
 \[\tilde{\bm x_i}=\sqrt{-E\frac{\partial^2\log f(y_i|\check\theta_i)}{\partial\theta^2}}\bm x_i,\]
 $\check\theta_i=\bm x_i^T\check{\bm\beta}$ with an initial estimate $\check{\bm\beta}$ available,
 and $\tilde{\bm X_d}$ is the design matrix for the transformed data \citep[see][]{Lee1987Diagnostic}.
 In this example, simple algebra yields $\tilde{\bm x_i}=\sqrt{e^{\bm x_i^T\hat{\bm\beta}_{\rm MLE}}}\bm x_i$.
 For the leverage score subsampling, we considered two scenario here. The first is that $\pi_i\propto \bm x_i(\bm X_d^T\bm X_d)^{-1}\bm x_i$ and the second is adjusted version with $\bm x_i,\bm X_d$ replaced by $\tilde{\bm x_i}$ and $\tilde{\bm X_d}$, respectively. We denote the seven cases as ``uniform", ``mV", ``mVc", ``mV-LD", ``mVc-LD", ``Lev-O" and ``Lev-A" respectively and present the results in Figure \ref{fig:1}.

\begin{figure}[!htp]
  \centering
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplot1-case1-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 1 (independent covariates)}
  \end{subfigure}
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplot1-case2-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 2 (highly correlated covariates)}
  \end{subfigure}\\[5mm]
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplot1-case3-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 3 (weakly correlated covariates)}
  \end{subfigure}
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplot1-case4-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 4 (unequal bounds of covariates)}
  \end{subfigure}\\[5mm]
  \caption{MSEs for different size ${r}$ with different distributions of covariates which are listed in the front of Section 5.}
  \label{fig:1}
\end{figure}
\end{example}

From Figure \ref{fig:1}, it is clear to see that for all the four data sets, $\bm\pi^{\mathrm{mV}}$ and  $\bm\pi^{\mathrm{mVc}}$ always result in smaller MSE than the uniform subsampling, which agrees with the theoretical result that they aim to minimize the asymptotic MSEs of the resultant estimator. If the components of $\bm x$ are independent, $\bm\pi^{\mathrm{mV}}$ and  $\bm\pi^{\mathrm{mVc}}$ have similar performances while for some highly correlated covariates this is not true. The reason is obvious, $MVc$-optimality reduces the impact of the data correlation structure, since the $\|\mathcal{J}_X^{-1}\bm x_i\|^2 $ in $MV$-optimality has been replaced by $\|\bm x_i\|^2$.
 It is worth noting that when $r$ is relative small, the subsampling strategy with slightly larger disturbance has a better performance than the smaller one.
 We guess the  disturbance in the criterion protects the estimator against the leverage points which lead the estimator much more stable. This can be regard as a windfall for adding $\delta$ in the criterion.

 Also note that Cases 1,3,4 perform quite well since the condition number of $\bm X_d$ is quite small and around at 5, hence when we choose $r=100$ can achieve satisfying results.
 However, for Case 2, the condition number is around at 40, which indicates 1000 subsamples are needed to well approximate $\hat{\bm\beta}_{\rm MLE}$. And we can easily see from Figure \ref{fig:1} that $\bm\pi^{\mathrm{mV}}$ performs better than  $\bm\pi^{\mathrm{mVc}}$ since $\kappa(\mathcal{J}_X^{-1})$ is usually great than 1 which echoes the results in subsection \ref{sec:r-size}. To see the influence of condition number clearly,  we first remove the effect of the pilot estimator and  did additional experiments by observing  the size of subsamples from 20 to 200. From the rule of thumb proposed in section \ref{sec:r-size}, the recommend sample size is $800 \sim 1000$. Hence we also report the MSE under the second case  when $r$ is selected from $200\sim 1000$.    These results are presented in Figure \ref{fig:5}.
 \begin{figure}[!htp]
  \centering
   \begin{subfigure}{0.32\textwidth}
    \includegraphics[width=\textwidth]{Rplot5-case2-1-8-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 2 ($\kappa\approx39$)}
  \end{subfigure}
  \begin{subfigure}{0.32\textwidth}
    \includegraphics[width=\textwidth]{Rplot5-case2-8-16-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 2 ($\kappa\approx39$)}
  \end{subfigure}
  \begin{subfigure}{0.32\textwidth}
    \includegraphics[width=\textwidth]{Rplot5-dim7-4-915-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 4 ($\kappa\approx5$)}
  \end{subfigure}
  \caption{MSEs for different size ${r}$ with different distributions of covariates.  }
  \label{fig:5}
\end{figure}

To see the effects of dimensionality, we also conduct the experiments for the case $p=2$ and $p=7$ when the covariates are generated from case 1. It worth mentioning that the condition number will decrease a little as the dimension decreases. Regardless of this, from Figure \ref{fig:6} it is clear to see that the case with $p=2$ performs better than $p=7$ under the same subsample size.
\begin{figure}[!htp]
  \centering
   \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\textwidth]{Rplot5-dim7-iid-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 1 ($\kappa\approx5, p=7$)}
  \end{subfigure}
  \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\textwidth]{Rplot5-dim2-iid-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 1 ($\kappa\approx3, p=2$)}
  \end{subfigure}
  \caption{MSEs for different size ${r}$ with different dimension of covariates.  }
  \label{fig:6}
\end{figure}
\fi


We consider both $\tilde{\pi}_i^{\mathrm{mV}}$ and $\tilde{\pi}_i^{\mathrm{mVc}}$ 
 
 
 and choose the value of $\delta$ to be $\delta=10^{-6}$.  
 
 
 For comparison, we also consider uniform subsampling, i.e., $\pi_i=1/n$ for all $i$ and
 
 the leverage subsampling strategy in \cite{Ma2015A} and \cite{Jia2014Influence}.
 For linear regression with covariate matrix $\bm X_d=(\bm x_1, ..., \bm x_n)^T$, leverage scores are defined as $h_i=\bm x_i(\bm X_d^T\bm X_d)^{-1}\bm x_i$. Leverage sampling defines the subsampling probabilities to be $\pi_i=h_i/\sum_{j=1}^nh_i=h_i/p$ \citep{Ma2015A}.
 
 
 For generalized linear models, leverage scores are defined by using the adjusted covariate matrix, namely, $\tilde{h}_i=\tilde{\bm x}_i(\tilde{\bm X}_d^T\tilde{\bm X}_d)^{-1}\tilde{\bm x}_i$, where $\tilde{\bm X}_d=(\tilde{\bm x}_1, ..., \tilde{\bm x}_n)^T$,
 
 \[\tilde{\bm x_i}=\sqrt{-E\Big[\frac{\partial^2\log f(y_i|\tilde\theta_i)}{\partial\theta^2}}\Big]\bm x_i,\]
 and  $\tilde\theta_i=\bm x_i^T\tilde{\bm\beta}$ with an initial estimate $\tilde{\bm\beta}$ 
 
 \citep[see][]{Lee1987Diagnostic}.
 In this example, simple algebra yields $\tilde{\bm x_i}=\sqrt{e^{\bm x_i^T\tilde{\bm\beta}_0}}\bm x_i$.
 For the leverage score subsampling, we considered both $h_i$ and $\tilde{h}_i$. 
 Here is a summary for the methods to be compared:
 UNIF, uniform subsample;
mV, $\pi_i=\tilde{\pi}_i^{\mathrm{mV}}$;
mVc, $\pi_i=\tilde{\pi}_i^{\mathrm{mVc}}$;
Lev, leverage sampling;
Lev-A, adjust leverage sampling.

 


We first consider the case with the first step sample size fixed. We let $r_0=200$, and second step sample size $r$ be 300, 500, 700, 1000, 1200 and 1400, respectively. Since the subsample probabilities do not depend on unknown parameters and no pilot subsamples are required for uniform and Lev methods, they are  implemented with subsample size $r+r_0$ for fair comparisons.

Here, we evaluate the performance of a sampling strategy $\bm\pi$ by the empirical mean squared error (MSE) of the resultant estimator:
$\text{MSE}=K^{-1}\sum_{k=1}^K\|{\bm\beta}_{\bm\pi}^{(k)}-\hat{\bm\beta}_{\rm MLE}\|$, 
where ${\bm\beta}_{\bm\pi}^{(k)}$ is the estimator from the $k$-th subsample with subsampling probability $\bm\pi$ and $\hat{\bm\beta}_{\rm MLE}$ is the MLE calculated from the whole dataset. We set $K = 1000$ throughout this section.
Computations are performed by {\verb"R"} \citep{Rpackage2018}. 

\begin{figure}[!htp]
  \centering
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplot1-case1-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 1 (independent covariates)}
  \end{subfigure}
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplot1-case2-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 2 (highly correlated covariates)}
  \end{subfigure}\\[5mm]
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplot1-case3-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 3 (weakly correlated covariates)}
  \end{subfigure}
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplot1-case4-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 4 (unequal bounds of covariates)}
  \end{subfigure}\\[5mm]
  \caption{MSEs for different second step subsample size ${r}$ with
    the first step subsample size being fixed at $r_0=200$.
     The different distributions of covariates which are listed in the front of Section 5.}
  \label{fig:3}
\end{figure}


 


Figure \ref{fig:3} gives the empirical MSEs. It is seen that for all the four data sets, subsampling methods based on $\tilde{\bm\pi}^{\mathrm{mV}}$ and  $\tilde{\bm\pi}^{\mathrm{mVc}}$ always result in smaller MSE than the uniform subsampling, which agrees with the theoretical result that they aim to minimize the asymptotic MSEs of the resultant estimator. 
If the components of $\bm x$ are independent, $\tilde{\bm\pi}^{\mathrm{mV}}$ and  $\tilde{\bm\pi}^{\mathrm{mVc}}$ have similar performances, while they may perform differently if some covariates are highly correlated. The reason is that $\tilde{\bm\pi}^{\mathrm{mVc}}$ reduces the impact of the data correlation structure since $\|\tilde{\mathcal{J}}_X^{-1}\bm x_i\|$ in $\tilde{\bm\pi}^{\mathrm{mV}}$ are replaced by $\|\bm x_i\|$ in $\tilde{\bm\pi}^{\mathrm{mVc}}$.
 
 

 For Cases 1, 3 and 4, empirical MSEs are small. This is because the condition number of $\bm X_d$ is quite small ($\approx$ 5) and a small subsample size $r=100$ produces satisfactory results.
 However, for Case 2, the condition number is large ($\approx$ 40), so a larger subsample size is needed to approximate $\hat{\bm\beta}_{\rm MLE}$ accurately. This agrees with the conclusion in Theorem~\ref{thm:r-size}.


 Another contribution of Theorem \ref{thm:CLT-2step} is to enable us to do inference on $\bm\beta$. Note that in subsampling setting, $r$ is much smaller than the full data size $n$. If $r=o(n)$, then $\hat{\bm\beta}_{\rm MLE}$ in Theorem \ref{thm:CLT-2step} can be replaced by the true parameter. 
 As an example, we take $\beta_2$ as a parameter of interest and construct 95\% confidence intervals for it. For this, the estimator given in \eqref{eq:24} is used to estimate variance-covariance matrices based on selected subsamples. For comparison, uniform subsampling method is also implemented.

 Table \ref{tab:tabl} reports the empirical coverage probabilities and average lengths over the four synthetic data sets with 
 the first step subsample size being fixed at $r_0=200$.
It is clear that $\tilde{\bm\pi}^{\mathrm{mV}}$ and  $\tilde{\bm\pi}^{\mathrm{mVc}}$ have similar performances and are uniformly better than the uniform subsampling method. As $r$ increases,  the lengths of confidence intervals decrease uniformly which echos the results of Theorem \ref{thm:CLT-2step}. The length of confidence interval in Case 2 is longer than those in other cases with the same subsample sizes. This is due to the fact that the condition number of $\bm X_d$ in Case 2 is bigger than that of $\bm X_d$ in other cases. This indicates that we should select a larger subsample when the condition number of the full dataset is bigger. This echos the results discussed in Section \ref{sec:r-size}.

\begin{table}[htbp]

  \centering
  \caption{Empirical coverage probability and average lengths of confidence intervals for 
    $\beta_2$. 
    The first step subsample size is fixed at $r_0=200$. }
    \setlength{\tabcolsep}{3.5mm}{
    \begin{tabular}{cccccccc}
    \hline
         & \multicolumn{1}{c}{method} & \multicolumn{2}{c}{mV} &  \multicolumn{2}{c}{mVc} &  \multicolumn{2}{c}{UNIF} \\


      & \multicolumn{1}{c}{r} & \multicolumn{1}{l}{cover rate} & \multicolumn{1}{l}{Length} & \multicolumn{1}{l}{Coverage} & \multicolumn{1}{l}{Length} &  \multicolumn{1}{l}{Coverage} & \multicolumn{1}{l}{Length} \\
    
    \hline
    \multirow{3}[2]{*}{case 1} & 300   & 0.954 & 0.2037    & 0.955 & 0.2066   & 0.952 & 0.2275  \\
          & 500   & 0.954 & 0.1684    & 0.945 & 0.1713   & 0.942 & 0.1924  \\
          & 1000  & 0.946 & 0.1254    & 0.938 & 0.1281    & 0.953 & 0.1471  \\
    
    \multirow{3}[2]{*}{case 2} & 300   & 0.961 & 1.9067    & 0.946 & 2.0776    & 0.950  & 2.2549  \\
          & 500   & 0.958 & 1.5470    & 0.948 & 1.7263   & 0.947 & 1.9082  \\
          & 1000  & 0.954 & 1.1379    & 0.948 & 1.2919    & 0.945 & 1.4559  \\
    
    \multirow{3}[2]{*}{case 3} & 300   & 0.959 & 0.1770    & 0.953  & 0.1816    & 0.939  & 0.2000  \\
          & 500   & 0.942 & 0.1451   & 0.949  & 0.1507    & 0.942  & 0.1693  \\
          & 1000  & 0.954 & 0.1082    & 0.954  & 0.1132    & 0.939  & 0.1291  \\
    
    \multirow{3}[2]{*}{case 4} & 300   & 0.955 & 0.2097   & 0.951 & 0.2179    & 0.953 & 0.2402  \\
          & 500   & 0.951 & 0.1721    & 0.956 & 0.1803    & 0.942 & 0.2033  \\
          & 1000  & 0.957 & 0.1276    & 0.960  & 0.1347    & 0.943 & 0.1552  \\
    
    \hline
    \end{tabular}}
  \label{tab:tabl}
\end{table}


 \iffalse
 To see the influence of condition number clearly, we consider the ideal case that $\hat{\bm\beta}_{\rm MLE}$ is used to calculate the subsampling probabilities and conduct Algorithm \ref{alg:1}  by using the optimal subsampling probabilities with $r$ ranges from 20 to 200. From the rule of thumb proposed in section \ref{sec:r-size}, the recommend sample size is $800 \sim 1000$. Hence we report the MSE under the second case  when $r$ is selected from $200\sim 1000$.    These results are presented in Figure \ref{fig:5}.
 \begin{figure}[!htp]
  \centering
   \begin{subfigure}{0.32\textwidth}
    \includegraphics[width=\textwidth]{Rplot5-case2-1-8-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 2 ($\kappa\approx39$)}
  \end{subfigure}
  \begin{subfigure}{0.32\textwidth}
    \includegraphics[width=\textwidth]{Rplot5-case2-8-16-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 2 ($\kappa\approx39$)}
  \end{subfigure}
  \begin{subfigure}{0.32\textwidth}
    \includegraphics[width=\textwidth]{Rplot5-dim7-4-915-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 4 ($\kappa\approx5$)}
  \end{subfigure}
  \caption{MSEs for different size ${r}$ with different distributions of covariates.  }
  \label{fig:5}
\end{figure}

To see the effects of dimensionality, we also conduct the experiments for the case $p=2$ and $p=7$ when the covariates are generated from case 1. To remove the impact of the pilot estimator, we also consider the case that $\hat{\bm\beta}_{\rm MLE}$ is known at beginning. It worth mentioning that the condition number will decrease a little as the dimension decreases. Regardless of this, from Figure \ref{fig:6} it is clear to see that the case with $p=2$ performs better than $p=7$ under the same subsample size.
\begin{figure}[!htp]
  \centering
   \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\textwidth]{Rplot5-dim7-iid-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 1 ($\kappa\approx5, p=7$)}
  \end{subfigure}
  \begin{subfigure}{0.45\textwidth}
    \includegraphics[width=\textwidth]{Rplot5-dim2-iid-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 1 ($\kappa\approx3, p=2$)}
  \end{subfigure}
  \caption{MSEs for different size ${r}$ with different dimension of covariates.  }
  \label{fig:6}
\end{figure}
\fi


Now we investigate the effect of different sample size allocations between the two steps. We calculate MSEs for various proportions of first step subsamples with fixed total subsample sizes. The results are given in Figure \ref{fig:4} with total subsample size $r_0+r$ = 800 and 1200, respectively. { Since all the cases have similar performance, we only demonstrate the performance of Case 4 here.}

It is worth noting that the two-step method outperforms the uniform subsampling method for all the four data sets, when $r_0/r\in[0.1,0.9]$. 
This indicates that the two-step approach is more efficient than the uniform subsampling.  
The two-step approach works the best when $r_0/r$ is around 0.2.

\begin{figure}[!htp]
  \centering
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplot2-case4-800-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 4 ($r_0+r=800$)}
  \end{subfigure}
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplot2-case4-1200-eps-converted-to.pdf}\\[-1cm]
    \caption{Case 4 ($r_0+r=1200$)}
  \end{subfigure}\\[5mm]
  \caption{MSEs vs proportions of the first step subsample with fixed
    total subsample sizes $r+r_0$.}
  \label{fig:4}
\end{figure}

To further evaluate the computational efficiency of the subsampling
strategies, we record the computing time  of five subsampling strategies ( uniform, $\pi^{\rm mV}$, $\pi^{\rm mVc}$, leverage score and adjust leverage score) by using \verb"microbenchmark()" function in R package \verb"microbenchmark" \citep{Olaf2018microbenchmark} and evaluating each subsampling strategy 50 times.
For fair comparison, R programming language is used to implement each strategy.  Computations were carried out on a
desktop running Window 10 with an Intel I7 processor and 32GB
memory.   Table~\ref{tab:tab2} shows the results
for Case 4 synthetic dataset with different  r and a fixed $r_0=200$. The computing time for using the full data
is also given  for comparisons.
It is not surprising to observe  that the uniform subsampling algorithm requires the
least computing time because it does not require an additional step to
calculate the subsampling probability. The algorithm based on $\bm{\pi}^{\mathrm{mV}}$
requires longer computing time than the algorithm based on
$\bm{\pi}^{\mathrm{mVc}}$, which agrees with the theoretical analysis
in Section~\ref{sec:two-step}. 


Table \ref{tab:tab2} reports the computing time with $p=7$ and the Leverage score is computed directly by the definition since $p=7$ is not that big to use the fast computing method mentioned in \cite{Drineas2011Fast}.
For fairness, we also consider the case with $p=40,n=100,000$, which is suitable to use the fast computing method for the Lev and Lev-A methods.
The first seven variables are generated as the case 4 and the rest are simply generated independently from $U(\left[0,1 \right])$. To ensure that the pilot estimator exists, $r_0$ is selected as 800.
The results are reported in Table \ref{tab:tab3}.

From Table \ref{tab:tab3} it is clear that all the subsampling algorithms take significantly less computing time compared to using the full data approach.
As the dimension increases, the computational advantage of $\bm\pi^{\rm mVc}$ is even more prominent.

\begin{table}[htbp]
  \centering
  \caption{Average Computing  time (in millisecond) for  synthetic dataset Case 4 with  different $r$ and a fixed $r_0=200$. }
  \setlength{\tabcolsep}{5.5mm}{
    \begin{tabular}{ccccccc}
    \hline
     r     & FULL  & UNIF  & mV    & mVc   & Lev   & Lev-A \\
    \hline
    \iffalse
    \multicolumn{1}{c}{\multirow{4}[0]{*}{10,000}} & 500   & 19.54022 & 1.936729 & 4.451563 & 4.202017 & 15.89786 & 18.52855 \\
          & 1000  & 19.86290 & 2.899916 & 5.591547 & 5.185855 & 16.94539 & 19.47171 \\
          & 1500  & 19.68008 & 3.895450 & 6.456192 & 6.160542 & 18.05423 & 20.75646 \\
          & 2000  & 20.33609 & 4.841714 & 7.529642 & 7.227572 & 19.01052 & 21.90701 \\
    \hline
    \multicolumn{1}{c}{\multirow{4}[0]{*}{50,000}} & 500   & 128.87450 & 2.642706 & 15.21282 & 11.93374 & 92.86035 & 105.2449 \\
          & 1000  & 129.41000 & 3.868140 & 15.70686 & 13.19126 & 91.98165 & 104.8732 \\
          & 1500  & 120.45070 & 5.006979 & 15.99610 & 13.6724 & 83.62497 & 95.64575 \\
          & 2000  & 120.03480 & 5.925759 & 16.69876 & 15.20175 & 85.32086 & 96.50285 \\
    \hline
    \fi
     500   & 279.358 & 2.690 & 28.263 & 21.621 & 190.222 & 215.302 \\
     1000  & 274.692 & 4.158 & 29.759 & 23.516 & 192.733 & 219.766 \\
     1500  & 277.668 & 5.705 & 30.902 & 24.287 & 193.545 & 218.556 \\
     2000  & 283.681 & 7.352 & 33.288 & 26.158 & 205.934 & 231.874 \\
    \hline
    \end{tabular}}
  \label{tab:tab2}
\end{table}


\begin{table}[htbp]
  \centering
  \caption{Average Computing  time (in millisecond) for synthetic dataset with full sample size $n=100,000$, dimension $p=40$, different $r$, and a fixed $r_0=800$.}
  \setlength{\tabcolsep}{5.5mm}{
    \begin{tabular}{ccccccc}
    \hline
     r     & FULL  & UNIF  & mV   & mVc  & Lev   & Lev-A \\
    \hline
    
    1000  & 5582.989  & 74.399  & 256.677  & 144.516  & 524.909  & 619.627  \\
           1500  & 5564.029  & 90.893  & 275.119  & 163.675  & 538.146  & 638.587  \\
           2000  & 5586.085  & 111.808  & 294.143  & 180.168  & 561.614  & 658.294  \\
           2500  & 5611.224  & 135.920  & 318.774  & 204.375  & 586.907  & 682.043  \\
    \hline
    \iffalse
    \multirow{4}[0]{*}{80} & 1000  & 4358.178  & 64.209  & 609.892  & 152.668  & 694.588  & 820.866  \\
          & 1500  & 4422.596  & 81.968  & 627.799  & 168.246  & 709.031  & 836.640  \\
          & 2000  & 4400.645  & 98.044  & 656.457  & 185.230  & 729.920  & 868.055  \\
          & 2500  & 4358.141  & 116.116  & 660.989  & 201.622  & 736.895  & 870.618  \\
    \hline
    \fi
    \end{tabular}}
  \label{tab:tab3}
\end{table}


\subsection{Real Data Studies} \label{sec:realdata}

\iffalse

\subsubsection{The Physicochemical Properties of Protein Tertiary Structure Data Set}
In the following, we illustrate our methods described in Section \ref{sec:two-step} by applying them to a data set from bioinformatics.
This data set is about physicochemical properties of protein tertiary structure taken from Critical Assessment of protein Structure Prediction (CASP),  which is a community-wide, worldwide experiment for protein structure prediction taking place every two years. The data set is from the fifth to ninth experiment and it is available from the UCI database \citep{Dua:2017} at \url{http://archive.ics.uci.edu/ml/datasets/Physicochemical+Properties+of+Protein+Tertiary+Structure}. 
The objective is to model the RMSD, a measure commonly used to evaluate the structure similarity from CASP.

 In order to predict the RMSD, we fit a linear regression model by using the following nine factors provided in the dataset as covariates: total surface area ($x_1$), non polar exposed area ($x_2$), fractional area of exposed non polar residue ($x_3$), fractional area of exposed non polar part of residue ($x_4$), molecular mass weighted exposed area ($x_5$),
average deviation from standard exposed area of residue ($x_6$),  Euclidean distance ($x_7$), secondary structure penalty ($x_8$), spacial distribution constraints ($x_9$).
Since the variables in the data set are on different scales, we normalize the dataset first and draw the inference by  simple linear model.
 
 A total of $n = 45,730$ decoys are recorded, with no missing values. 
Full data least-squares produces  
the following model:
\[Y=1.04X_1+0.34X_2+0.18X_3-0.98X_4-0.37X_5-0.27X_6-0.04X_7+0.13X_8-0.11X_9+\varepsilon.\]

Similarly to synthetic datasets, we also compare our method with uniform subsampling method and the leverage score subsampling method and report the results actual subset under various measurement budget ($r$) varying from 500 to 2000.
In our experiments we compare four subsampling strategies ( uniform, $\pi^{mV}$, $\pi^{mVc}$ and leverage score) when the actual subset under various measurement budget ($r$) varying from 700 to 2000.

The  mean-square prediction error $\frac{1}{n}\|\bm y-\bm X_d\tilde{\bm\beta}\|^2$ is reported in Figure \ref{fig:7}.


\begin{figure}[!htp]
  \centering
    \includegraphics[scale=0.6]{Rplotrd-case1new-eps-converted-to.pdf}
  \caption{Prediction MSEs for different second step subsample size ${r}$ with
    the first step subsample size being fixed at $r_0=300$.}
  \label{fig:7}
\end{figure}

{\blue From Figure \ref{fig:7}, it is clear to see that  the $\bm{\pi}^{\mathrm{mV}}$, $\bm{\pi}^{\mathrm{mVc}}$ and leverage score  subsampling method are uniformlly better than the uniform subsampling method. As $r$ increases,  our methods  have a competitive performance as the leverage score  subsampling method which is the best among the five subsampling strategies in this example. It worth to mention that $\pi^{mVc}$ has quite lower computational complexity compare with the leverage score method. Thus $\pi^{mVc}$ is also a feasible solution especially for super large dataset. }
\fi


\subsubsection{The Demand for Medical Care Data Set}
In the following, we illustrate our methods by applying them to a cross-sectional data set from health economics.
\cite{Deb1997DEMAND} analyzed this dataset on 4406 individuals, aged 66 and over, who are covered
by Medicare, a public insurance program.
The data are available from the data archive  of the Journal of Applied Econometrics at \url{http://qed.econ.queensu.ca/jae/1997-v12.3/deb-trivedi/}.
The objective is to model the demand for medical care by the covariates
available for the patients, where the demand is captured by the
number of physician/non-physician office and hospital outpatient visits. Here, we use the number of physician office visits \verb"ofp" as the
dependent variable and  take the following six factors as covariates: number of hospital stays ($x_1$), self-perceived health status
($x_2$), number of chronic conditions ($x_3$),  gender ($x_4$), number of years of education ($x_5$), and private
insurance indicator ($x_6$).
An intercept term is also included in the regression model. 
In order to capture the relationship between the number of physician office visits
and all covariates, we use Poisson regression to model the data.
We implement our method with $r_0=200$ in first step and let $r$ vary 
in the second stage. 
For comparison, the uniform subsampling method and the leverage score subsampling methods are also implemented. The  mean-square prediction errors (MSPEs) $n^{-1}\|\bm y-\exp(\bm X_d\tilde{\bm\beta})\|^2$ and MSEs are reported in Figure \ref{fig:8}. It is seen that $\bm{\pi}^{\mathrm{mV}}$ dominates other methods in terms of both MSE and MSPE. 


\begin{figure}[!htp]
  \centering
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplotrd1-mse-eps-converted-to.pdf}\\[-1cm]
    \caption{MSE}
  \end{subfigure}
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplotrd1-mspe-eps-converted-to.pdf}\\[-1cm]
    \caption{$\log_{10}({\rm MSPE})$}
  \end{subfigure}\\[5mm]
  \caption{MSEs and MSPEs for different second step subsample sizerwith the first step subsample sizebeing fixed at $r_0= 200$.}
  \label{fig:8}
\end{figure}

\iffalse
\begin{table}[htbp]
  \centering
  \caption{MSEs for different second step subsample size ${r}$ with
    the first step subsample size being fixed at $r_0=200$}
    \begin{tabular}{cccccc}
    \hline
    r     & UNIF  & mV  & mVc  & Lev   & Lev-A \\
    \hline
    300   & 0.146  & 0.089  & 0.120  & 0.114  & 0.137  \\
    500   & 0.105  & 0.053  & 0.085  & 0.081  & 0.097  \\
    700   & 0.081  & 0.038  & 0.066  & 0.059  & 0.073  \\
    1000  & 0.062  & 0.029  & 0.050  & 0.045  & 0.053  \\
    1200  & 0.052  & 0.024  & 0.041  & 0.041  & 0.046  \\
    1400  & 0.046  & 0.020  & 0.036  & 0.033  & 0.040  \\
    \hline
    \end{tabular}
  \label{tab:tab4}
\end{table}
\fi


\subsubsection{The Gas Sensor Array Drift Data Set }
In the following, we illustrate our methods described in Section \ref{sec:two-step} by applying them to a data set from chemistry.
This data set is about 16 chemical sensors exposed to the mixture of Ethylene and CO at varying concentrations in air.
A total of $n = 4,208,261$ readings are recorded for each sensor, with no missing values.
Readers may refer to \cite{Fonollosa2015Reservoir} for more information about this data set. 
The objective is to predict the observations from last sensor by using the other sensors as covariates.
Using domain knowledge \citep{Goodson2011Mathematical}, a log transformation is taken to the sensors readings.
In addition, we drop out the first 20,000 data points which correspond to the system run-in time. 
In accordance with the irrational state of data acquisition, the second sensor's readings are excluded due to the reason that $20\%$ of the original reading for gas concentrations are negative.

A linear regression model is used to model the last sensor's readings with $n=4,188,261$ data points and $p=14$ covariates. 

  


Similar to synthetic datasets, we also compare our method with uniform subsampling and the leverage score subsampling methods, and report the results  for $r$ varying from 500 to 2000.

The  MSPEs $n^{-1}\|\bm y-\bm X_d\tilde{\bm\beta}\|^2$ and MSEs are presented in Figure \ref{fig:9}. It is seen that in terms of MSE $\bm{\pi}^{\mathrm{mV}}$ is the best; in terms of MSPE, $\bm{\pi}^{\mathrm{mV}}$ and $\bm{\pi}^{\mathrm{mVc}}$ perform similarly, and they are better than the leverage score subsampling for larger values of $r$.

\begin{figure}[!htp]
  \centering
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplotrd2-mse-eps-converted-to.pdf}\\[-1cm]
    \caption{MSE}
  \end{subfigure}
  \begin{subfigure}{0.49\textwidth}
    \includegraphics[width=\textwidth]{Rplotrd2-mspe-eps-converted-to.pdf}\\[-1cm]
    \caption{$\log_{10}({\rm MSPE})$}
  \end{subfigure}\\[5mm]
  \caption{MSEs and MSPEs for different second step subsample size $r$ with the first step subsample size being fixed at $r_0= 200$.}
  \label{fig:9}
\end{figure}


\iffalse
\begin{figure}[!htp]
  \centering
    \includegraphics[scale=0.6]{Rplotrd_chem_mspe-eps-converted-to.pdf}
  \caption{Log prediction MSEs for different second step subsample size ${r}$ with
    the first step subsample size being fixed at $r_0=200$.}
  \label{fig:8}
\end{figure}


\begin{table}[htbp]
  \centering
  \caption{MSEs for different second step subsample size ${r}$ with
    the first step subsample size being fixed at $r_0=200$}
     \setlength{\tabcolsep}{10mm}{
    \begin{tabular}{ccccc}
    \hline
    r     & UNIF  & mV  & mVc  & Lev    \\
    \hline
    500   & 0.122  & 0.079  & 0.096  & 0.086  \\
    700   & 0.097  & 0.057  & 0.066  & 0.064  \\
    1000  & 0.072  & 0.038  & 0.046  & 0.050  \\
    1200  & 0.064  & 0.032  & 0.039  & 0.042  \\
    1500  & 0.052  & 0.026  & 0.029  & 0.034  \\
    1700  & 0.046  & 0.023  & 0.027  & 0.031  \\
    2000  & 0.041  & 0.019  & 0.022  & 0.026  \\
    \hline
    \end{tabular}}
  \label{tab:tab5}
\end{table}
\fi


\section{Concluding remarks} \label{sec:summary}


  In this paper, we study statistical properties of a general subsampling algorithm for generalized linear models.


Based on the asymptotic distribution, we derive optimal subsampling probabilities under the A-optimality criterion. 

{
For practical implementation, a two-step algorithm 
is developed, for which we also derive the theoretical properties. As shown in numerical experiments, our algorithm is computationally feasible for super large datasets and it produces more efficient inference on parameter of interest compared with the uniform subsampling method. }


Subsampling of big data is a cutting-edge problem and more research is needed for more complicated model setups and data sets, such as Cox model, panel count data, and time series data. For these examples, the score function is a sum of dependent random variables, so the optimal subsample probabilities have more complicated structures.


\section{Appendix}

\subsection{Some Lemmas}

To prove Theorem \ref{thm:as-general-alg}, we begin with the following lemma which gives upper bound of the absolute moments of exponential family random variables.


\begin{lemma}\label{lem-moment} (Lemma III.2 in \cite{Blazere2014Oracle}) Let $Y$ be a random variable with density or mass function $f(y|\eta ) = c(y)\exp (y\eta  - A(\eta ))$ and assume $\left| \eta  \right| \le T$ are fulfilled. For all $t \in \mathbb{N}^{*}$ there exists a constant $C_{T}$ (which depends only on $T$) such that ${\rm{E}}|Y{|^t} \le t!C_T^t$.
\end{lemma}

With Assumption (H.3), $\left| \eta  \right| \le T$ holds automatically by observing that
\[\left| \eta  \right| = |{\bm \beta ^T}x| \le BL: = T.\]
Thus we could apply the above Lemma to the proof of next Lemma.

\begin{lemma}\label{lem:lem2}
If Assumptions (H.1)--(H.5) hold, then as $n\rightarrow\infty$,
\[\frac{1}{n^2}\sum^n_{i=1}\frac{\{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\}^k\bm x_i\bm x_i^T}{\pi_i}=O_P(1),\]
for $k=1,2,3$.
\end{lemma}

\begin{proof}
Direct calculation shows that for any component of $\frac{1}{n^2}\sum^n_{i=1}\frac{\{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\}^k\bm x_i\bm x_i^T}{\pi_i}$,
\[\frac{1}{n^2}\sum^n_{i=1}\frac{\{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\}^k x_{ij_1} x_{ij_2}}{\pi_i}\le \left(\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^2}{n\pi_i}\right)\left(\sum^n_{i=1}\frac{\{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\}^k}{n}\right).\]

From Lemma~\ref{lem-moment}, let $Z_n := \sum^n_{i=1}\frac{\{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\}^k}{n}$, we have
\[\begin{split}
E{Z_n} &= \frac{1}{n}\sum\limits_{i = 1}^n {E|} {y_i} - \dot \psi (\widehat{ \bm \beta} _{{\rm{MLE}}}^T{\bm x_i}){|^k} \le \frac{1}{n}\sum\limits_{i = 1}^n {E(|} {y_i}| + |\dot \psi (\widehat{ \bm \beta} _{{\rm{MLE}}}^T{\bm x_i})|{)^k} \le \frac{1}{n}\sum\limits_{i = 1}^n {E(|} {y_i}| + M{)^k}\\
& \le \frac{1}{n}\sum\limits_{i = 1}^n {\sum\limits_{t = 0}^k {\left( {\begin{array}{*{20}{c}}
k\\
t
\end{array}} \right)} E} |{y_i}|{M^{k - t}} \le \frac{1}{n}\sum\limits_{i = 1}^n {\sum\limits_{t = 0}^k {\left( {\begin{array}{*{20}{c}}
k\\
t
\end{array}} \right)} t!C_T^t} {M^{k - t}} = \sum\limits_{t = 0}^k {\left( {\begin{array}{*{20}{c}}
k\\
t
\end{array}} \right)} t!C_T^t{M^{k - t}},
\end{split} \]
where $M := \mathop {\sup }\limits_{\bm x} |\dot \psi (\widehat{ \bm \beta}_{{\rm{MLE}}}^T\bm x)|$. It is clear that $M < \infty $ by applying (H.1)-(H.3) and Remark~\ref{rm:analytic}.

By Chebyshev's inequality, we know that $Z_n =O_p(1)$ as $n\rightarrow\infty$. See also the similar proof of Theorem 3.1 in \cite{zhang2017elastic}. And from (H.1) and (H.5), $M_n:=\max_{i=1,\ldots,n}\|\bm x_i\|^2/(n\pi_i)=O_P(1)$. Then we are ready to show that $Z_n M_n$ is $O_P(1)$.

Note that both $Z_n$ and $M_n$ are positive and $O_P(1)$. For any given ${\varepsilon _1},{\varepsilon _2}>0$, there exist $\tau_1$ and $\tau_2$ such that $P\left(Z_n<\tau_1\right)\ge 1- \varepsilon _1$, $P\left(M_n<\tau_2\right)\ge 1- \varepsilon _2$. Then,
\[\begin{split}
P\left( {{Z_n}{M_n} \ge {\tau _1}{\tau _2}} \right) &= P\left( {{Z_n}{M_n} \ge {\tau _1}{\tau _2},{Z_n} < {\tau _1}} \right) + P\left( {{Z_n}{M_n} \ge {\tau _1}{\tau _2},{Z_n} \ge {\tau _1},{M_n} \ge {\tau _2}} \right) + \\
& + P\left( {{Z_n}{M_n} \ge {\tau _1}{\tau _2},{Z_n} \ge {\tau _1},{M_n} < {\tau _2}} \right)\\
& \le P\left( {{M_n} \ge {\tau _2}} \right) + \min \{ P\left( {{M_n} \ge {\tau _2}} \right),P\left( {{Z_n} \ge {\tau _1}} \right)\}  + P\left( {{Z_n} \ge {\tau _1}} \right)\\
& \le {\varepsilon _2} + \min \{ {\varepsilon _2},{\varepsilon _1}\}  + {\varepsilon _1}.
\end{split} \]
Let $\tau  = {\tau _1}{\tau _2},\varepsilon  = {\varepsilon _2} + \min \{ {\varepsilon _2},{\varepsilon _1}\}  + {\varepsilon _1}$. Thus $P\left( {Z_n}{M_n} < \tau \right) \ge 1 - \varepsilon $, namely, $Z_n M_n=O_P(1)$.
\end{proof}

\begin{remark}\label{rm:analytic}
By the fact that $\lambda (\theta ): = \exp (\psi (\theta ))$ is analytic in $\textrm{Int}(\Theta)$ (see Theorem 2.7 in \cite{brown1986fundamentals}), Cauchy's integral formula tells us that all its higher derivatives exist and are continuous. Therefore, the derivatives $\dot\psi(t)$,~$\ddot\psi(t)$,~$\dddot{\psi}(t)$ are continuous in $t$, and $\dot\psi(t)$,$\ddot\psi(t)$ are bounded on the compact set, which follows by a well-known property that every real-valued continuous function on a compact set is necessarily bounded.
\end{remark}

\begin{lemma}\label{lem:lem1}
If Assumptions (H.1)--(H.5) hold, then as $n\rightarrow\infty$, conditionally on $\mathcal{F}_n$ in probability,
\begin{align}
  \check{\mathcal{J}}_X-\mathcal{J}_X&=O_{P|\mathcal{F}_n}(r^{-1/2}),\label{eq:ob-infor}\\
   \frac{1}{n}\frac{\partial L^*(\hat{\bm\beta}_{\rm MLE})}{\partial\bm\beta} &=O_{P|\mathcal{F}_n}(r^{-1/2}),\label{eq:ob-score}
\end{align}
where 
\begin{equation*}
  \check{\mathcal{J}}_X=
  \frac{1}{n}\frac{\partial^2 L^*(\hat{\bm\beta}_{\rm MLE})}
  {\partial\bm\beta\partial\bm\beta^T}
  =\frac{1}{nr}\sum^r_{i=1}\frac{\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i^*)\bm x_i^*(\bm x_i^*)^T}{\pi^*_i}.
\end{equation*}
\end{lemma}
\begin{proof}
By the definition condition expectation and towering property of flirtations, it yields that
\[E(\check{\mathcal{J}}_X|\mathcal{F}_n)=\mathcal{J}_X.\]
For any component $\tilde{\mathcal{J}}_X^{j_1j_2}$ of $\tilde{\mathcal{J}}_X$ where $1\le j_1,j_2\le p$,
\begin{align*}
  {\rm E}\left(\check{\mathcal{J}}_X^{j_1j_2}-\mathcal{J}_X^{j_1j_2}\Big|\mathcal{F}_n\right)^2
  =&\frac{1}{r}\sum^n_{i=1}\pi_i
     \left\{\frac{\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)x_{ij_1} x_{ij_2}}{n\pi_i}
     -\mathcal{J}_X^{j_1j_2}\right\}^2\\
  =&\frac{1}{r}\sum^n_{i=1}\pi_i
     \left\{\frac{\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)x_{ij_1} x_{ij_2}}{n\pi_i}\right\}^2
     -\frac{1}{r}\left(\mathcal{J}_X^{j_1j_2}\right)^2\\
  \le&\frac{1}{rn^2}\sum^n_{i=1}
     \frac{[\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)x_{ij_1} x_{ij_2}]^2}{\pi_i}
     -\frac{1}{r}(\mathcal{J}_X^{j_1j_2})^2\\
  \le&\frac{1}{r}\left(\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^2}{n\pi_i}\right)\sum_{i=1}^n\frac{(\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i))^2}{n}
      -\frac{1}{r}(\mathcal{J}_X^{j_1j_2})^2.
\end{align*}
 The last inequality holds by Assumptions (H.1) and (H.3).

Combining the facts that the parameter space is compact, $\hat{\bm\beta}_{\rm MLE}$ is the unique global maximum of log-likelihood function and Assumption (H.1), we have that $\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)$'s are almost surely bounded due to (H.1)-(H.3) and Remark~\ref{rm:analytic}. Then
it follows by Assumption (H.5) that
\begin{equation}\label{eq:J-J}
{\rm E}\left(\check{\mathcal{J}}_X^{j_1j_2}-\mathcal{J}_X^{j_1j_2}\Big|\mathcal{F}_n\right)^2=O_P(r^{-1}).
\end{equation}
 From Chebyshev's inequality, it is proved that Equation (\ref{eq:ob-infor}) holds.

To prove Equation (\ref{eq:ob-score}), let $t_i(\bm\beta)=y_i\bm\beta^T\bm x_i-\psi(\bm\beta^T\bm x_i)$, $t_i^*(\bm\beta)=y_i^*\bm\beta^T\bm x_i^*-\psi(\bm\beta^T\bm x_i^*)$, then
\begin{equation*}
L^*(\bm\beta)=\frac{1}{r}\sum^r_{i=1} \frac{t_i^*(\bm\beta)}{\tilde\pi_i^*}, \quad \text{and} \quad L(\bm\beta)=\sum^n_{i=1} t_i(\bm\beta).
\end{equation*}
Under the conditional distribution of subsample given $\mathcal{F}_n$,
\[
E\left\{\frac{L^*(\bm\beta)}{n}-\frac{L(\bm\beta)}{n}\bigg|\mathcal{F}_n\right\}^2
    = \frac{1}{n^2}\sum^n_{i=1}  \frac{t_i^2(\bm\beta)}{\pi_i}+\left(\frac{1}{n}\sum^n_{i=1} t_i(\bm\beta)\right)^2.
\]
It can be shown that
\[\begin{split}
|t_i(\bm\beta)|&\le |y_i\bm\beta^T\bm x_i-\dot{\psi}(\bm\beta^T\bm x_i)\bm\beta^T\bm x_i|+|\dot{\psi}(\bm\beta^T\bm x_i)\bm\beta^T\bm x_i-\psi(\bm\beta^T\bm x_i)|\\
&\le |y_i\bm\beta^T\bm x_i-\dot{\psi}(\bm\beta^T\bm x_i)\bm\beta^T\bm x_i|+|\dot{\psi}(\bm\beta^T\bm x_i)\bm\beta^T\bm x_i|+|\psi(\bm\beta^T\bm x_i)|\\
&\le C_1|y_i-\dot{\psi}(\bm\beta^T\bm x_i)|+C_2,
\end{split}\]
where $C_1=BL$ and $C_2$ is a positive constant. The last inequality holds by noting Assumptions (H.1) and (H.3).
  Therefore, as $n\rightarrow\infty$,
\[\begin{split}
\frac{1}{n^2}\sum^n_{i=1}\pi_i  \frac{t_i^2(\bm\beta)}{\pi_i^2}&\le \frac{1}{n^2}\sum^n_{i=1}  \left[\frac{C_1^2(y_i-\dot{\psi}(\bm\beta^T\bm x_i))^2}{\pi_i}+\frac{2C_1C_2|y_i-\dot{\psi}(\bm\beta^T\bm x_i)|}{\pi_i}+\frac{C_2^2}{\pi_i}\right],\\
\left(\frac{1}{n}\sum^n_{i=1} t_i(\bm\beta)\right)^2
&\le \left(\frac{1}{n}\sum^n_{i=1} C_1|y_i-\dot{\psi}(\bm\beta^T\bm x_i)|\right)^2 + \frac{2C_2 C_1}{n}\sum^n_{i=1}|y_i-\dot{\psi}(\bm\beta^T\bm x_i)| +C_2^2.


\end{split} \]

From Assumptions (H.1), (H.5) and Lemma \ref{lem:lem2},
\begin{equation}\label{eq:l2error}
E\left\{\frac{L^*(\bm\beta)}{n}-\frac{L(\bm\beta)}{n}\bigg|\mathcal{F}_n\right\}^2\le\frac{1}{r}\left\{\frac{1}{n^2}\sum^n_{i=1}\pi_i  \frac{t_i^2(\bm\beta)}{\pi_i^2}+\left(\frac{1}{n}\sum^n_{i=1} t_i(\bm\beta)\right)^2 \right\} =O_p(r^{-1}).
\end{equation}
Thus the desired result (\ref{eq:ob-score}) follows from Chebyshev's Inequality.
\end{proof}

\subsection{Proof of \textbf{Theorem~\ref{thm:as-general-alg}}}

Based on the lemmas in the above subsection, now we are ready to prove Theorem \ref{thm:as-general-alg}.
\begin{proof}
As $r\rightarrow\infty$, by (\ref{eq:l2error}), we have that $n^{-1}L^*(\bm\beta)-n^{-1}L(\bm\beta)\rightarrow0$ in conditional probability given $\mathcal{F}_n$. Note that the parameter space is compact and $\hat{\bm\beta}_{\rm MLE}$ is the unique global maximum of the continuous convex function $L(\bm\beta)$. Thus, from Theorem 5.9 and its remark of \cite{Vaart2000Asymptotic}, as $n\rightarrow\infty, r\rightarrow\infty$, conditionally on $\mathcal{F}_n$ in probability,
\begin{equation}\label{eq:12}
  \|\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}\|=o_{P|\mathcal{F}_n}(1).
\end{equation}

Using Taylor's theorem for random variables \citep[see][Chapter 4]{Ferguson1996A},
\begin{align}
0=\frac{\dot{L}^*_j(\tilde{\bm\beta})}{n}
=&\frac{\dot{L}^*_j(\hat{\bm\beta}_{\rm MLE})}{n}
   +\frac{1}{n} \frac{\partial\dot{L}^*_j(\hat{\bm\beta}_{\rm MLE})}{\partial \bm\beta^T}(\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE})
   +\frac{1}{n} R_j,
   \label{eq:taylor-expansion}
\end{align}
where $\dot{L}^*_j({\bm\beta})$ is the partial derivative of $L^*({\bm\beta})$ with respect to $\beta_j$, and the Lagrange remainder
\begin{equation*}
  \frac{1}{n}R_j=\frac{1}{n}(\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE})^T
   \int_0^1\int_0^1\frac{\partial^2\dot{L}^*_j\{\hat{\bm\beta}_{\rm MLE}
   +uv(\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE})\}}{\partial \bm\beta\partial
   \bm\beta^T}v du dv\
(\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}).
\end{equation*}
From (H.1), (H.3) and (H.5),
\begin{align*}
  \left\|\frac{\partial^2\dot{L}^*_j(\bm\beta)}
  {\partial\bm\beta\partial\bm\beta^T}\right\|=
  &\frac{1}{r}\left\|\sum_{i=1}^r \frac{-\dddot{\psi}(\bm\beta^T\bm x_i)}{\tilde\pi^*_i}x^*_{ij}\bm x^*_i{\bm x^*_i}^T\right\|
    \le\frac{C_3}{r}\sum_{i=1}^r\frac{\|\bm x^*_i\|^3}{\pi^*_i}
\end{align*}
for all $\bm\beta$, where $C_3$ is the upper bound of $\|\dddot{\psi}(\bm\beta^T\bm x_i)\|$, From (H.1-3) and Remark~\ref{rm:analytic}, it is known that $C_3$ is a constant. Thus
\begin{align*}
  \frac{1}{n}\left\|\int_0^1\int_0^1\frac{\partial^2\dot{L}^*_j\{\hat{\bm\beta}_{\rm MLE}
   +uv(\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE})\}}{\partial \bm\beta\partial\bm\beta^T}vdudv\ \right\|
  \le\frac{C_3}{2r}\sum_{i=1}^r\frac{\|\bm x^*_i\|^3}{n\pi^*_i}=O_{P|\mathcal{F}_n}(1),
\end{align*}
where the last equality holds by noticing Assumptions (H.1), (H.5) and the fact that as $\tau\rightarrow\infty$,
\begin{align}\label{eq:53}
  P\left(\frac{1}{nr}\sum_{i=1}^r\frac{\|\bm x^*_i\|^3}{\pi^*_i}\ge\tau
  \Bigg|\mathcal{F}_n\right) \le\frac{1}{nr\tau}\sum_{i=1}^r
  E\left(\frac{\|\bm x^*_i\|^3}{\pi^*_i}\Bigg|\mathcal{F}_n\right)
  =\frac{1}{n\tau}\sum_{i=1}^n\|\bm x_i\|^3\rightarrow0,
\end{align}
in probability.

Combining the above equations with the Taylor's expansion (\ref{eq:taylor-expansion}), we have
\begin{equation}\label{eq:13-2}
  \tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}=
  -\check{\mathcal{J}}_X^{-1}\left\{\frac{\dot {L}^*(\hat{\bm\beta}_{\rm MLE})}{n}
    +O_{P|\mathcal{F}_n}(\|\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}\|^2)\right\}.
\end{equation}
From  Lemma~\ref{lem:lem1} and Assumption (H.4), it is obvious that $\check{\mathcal{J}}_X^{-1} =O_{P|\mathcal{F}_n}(1)$. Therefore,
\begin{equation*}
  \tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}
  =O_{P|\mathcal{F}_n}(r^{-1/2})+o_{P|\mathcal{F}_n}(\|\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}\|),
\end{equation*}
which implies that
\[  \tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}=O_{P|\mathcal{F}_n}(r^{-1/2}).\]
\end{proof}

\subsection{Proof of \textbf{Theorem~\ref{thm:CLT}}}
\begin{proof}
  Note that
\begin{equation}\label{eq:linder1}
  \frac{\dot L^*(\hat{\bm\beta}_{\rm MLE})}{n}
  =\frac{1}{r}\sum^r_{i=1}\frac{\{y^*_i-\dot{\psi}^*_i(\hat{\bm\beta}_{\rm MLE}^T\bm x_i^*)\}\bm x^*_i}{n\pi^*_i}
  =:\frac{1}{r}\sum^r_{i=1}\eta_i.
\end{equation}
It can be seen that given $\mathcal{F}_n$, $\eta_1, \ldots, \eta_r$ are i.i.d  random variables with mean $\mathbf{0}$ and variance
\begin{align}\label{eq:linder2}
  &\text{var}(\eta_1|\mathcal{F}_n)=
 \frac{1}{n^2}\sum^n_{i=1}\frac{\{y_i-\dot{\psi}_i(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^2
 \bm x_i\bm x_i^T}{\pi_i}. 
\end{align}
Then from Lemma \ref{lem:lem2}, it is known that $\text{var}(\eta_i|\mathcal{F}_n)=O_P(1)$ as $n\rightarrow\infty$.

Meanwhile, for some $\delta>0$ and every $\varepsilon>0$,
\begin{equation}\label{eq:dbclt}
\begin{aligned}
\sum^r_{i=1} E\{\|r^{-1/2}\eta_i\|^2
    I(\|\eta_i\|>r^{1/2}\varepsilon)|\mathcal{F}_n\} &\le\frac{1}{r^{1+\delta/2}\varepsilon^{\delta}}
    \sum^r_{i=1}E\{\|\eta_i\|^{2+\delta}
    I(\|\eta_i\|>r^{1/2}\varepsilon)|\mathcal{F}_n\}\\
  &\le\frac{1}{r^{1+\delta/2}\varepsilon^{\delta}}
    \sum^r_{i=1} E(\|\eta_i\|^{2+\delta}|\mathcal{F}_n)\\
  &\le \frac{1}{r^{\delta/2}}\frac{1}{n^{2+\delta}}\frac{1}{\varepsilon^{\delta}}
    \sum^n_{i=1}\frac{\{y_i-\dot{\psi}_i(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^{2+\delta}
    \|\bm x_i\|^{2+\delta}}{\pi_i^{1+\delta}}\\
  &\le\frac{1}{r^{\delta/2}}\frac{1}{\varepsilon^{\delta}}
\left(\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^{2+\delta}}{(n\pi_i)^{1+\delta}}\right)\sum^n_{i=1}\frac{\{y_i-\dot{\psi}_i(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^{2+\delta}}{n}.
\end{aligned}
\end{equation}

Then  we obtain that
\begin{align*}
\sum^r_{i=1} E\{\|r^{-1/2}\eta_i\|^2
    I(\|\eta_i\|>r^{1/2}\varepsilon)|\mathcal{F}_n\} \le \frac{1}{r^{\delta/2}}\frac{1}{\varepsilon^{\delta}}
O_P(1)\cdot O_P(1)=o_P(1),
\end{align*}
where the last equality holds by using Assumptions (H.1), (H.5) and Lemma \ref{lem:lem2} with $k=2+\delta$.
This shows that the Lindeberg-Feller conditions are satisfied in probability.
From (\ref{eq:linder1}) and (\ref{eq:linder2}), by the Lindeberg-Feller central limit theorem \citep[Proposition 2.27 of][]{Vaart2000Asymptotic}, conditionally on $\mathcal{F}_n$, it can be shown that
\begin{equation*}
  \frac{1}{n}V_c^{-1/2}\dot{L}^*(\hat{\bm\beta}_{\rm MLE})=
  \frac{1}{r^{1/2}}\{\text{var}(\eta_i|\mathcal{F}_n)\}^{-1/2}\sum^r_{i=1}\eta_i
  \rightarrow N(0,{I}),
\end{equation*}
in distribution.
From Lemma~\ref{lem:lem1}, (\ref{eq:13-2}) and Theorem \ref{thm:as-general-alg}, we have
\begin{equation}\label{eq:15}
  \tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}=
  -\frac{1}{n}\check{\mathcal{J}}_X^{-1}\dot{L}^*(\hat{\bm\beta}_{\rm MLE})+O_{P|\mathcal{F}_n}(r^{-1}).
\end{equation}
From (\ref{eq:ob-infor}) in Lemma~\ref{lem:lem1}, it follows that
\begin{align}\label{eq:7}
  \check{\mathcal{J}}_X^{-1}-\mathcal{J}_X^{-1}
  &=-\mathcal{J}_X^{-1}(\check{\mathcal{J}}_X-\mathcal{J}_X)\check{\mathcal{J}}_X^{-1}
  =O_{P|\mathcal{F}_n}(r^{-1/2}).
\end{align}
Based on Assumption (H.4) and (\ref{eq:linder2}), it can be proved that
\begin{equation*}
  V={\mathcal{J}}_X^{-1}V_c{\mathcal{J}}_X^{-1}
  =\frac{1}{r}{\mathcal{J}}_X^{-1}\left(rV_c\right)
  {\mathcal{J}}_X^{-1}=O_{P}(r^{-1}).
\end{equation*}
Thus,
\begin{align*}
  V^{-1/2}(\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE})
  &=-V^{-1/2}n^{-1}\check{\mathcal{J}}_X^{-1}\dot{L}^*(\hat{\bm\beta}_{\rm MLE})
    +O_{P|\mathcal{F}_n}(r^{-1/2})\\
  &=-V^{-1/2}\mathcal{J}_X^{-1}n^{-1}\dot{L}^*(\hat{\bm\beta}_{\rm MLE})
    -V^{-1/2}(\check{\mathcal{J}}_X^{-1}-\mathcal{J}_X^{-1})n^{-1}\dot{L}^*(\hat{\bm\beta}_{\rm MLE})
    +O_{P|\mathcal{F}_n}(r^{-1/2})\\
  &=-V^{-1/2}\mathcal{J}_X^{-1}V_c^{1/2}V_c^{-1/2}n^{-1}\dot{L}^*(\hat{\bm\beta}_{\rm MLE})
    +O_{P|\mathcal{F}_n}(r^{-1/2}).
\end{align*}
So the result in \eqref{normal} of Theorem~\ref{thm:CLT} follows by applying Slutsky's Theorem \citep[Theorem 6, Section 6 of][]{Ferguson1996A} and the fact that
\begin{equation*}
  V^{-1/2}\mathcal{J}_X^{-1}V_c^{1/2}(V^{-1/2}\mathcal{J}_X^{-1}V_c^{1/2})^T
  =V^{-1/2}\mathcal{J}_X^{-1}V_c^{1/2}V_c^{1/2}\mathcal{J}_X^{-1}V^{-1/2}={I}.
\end{equation*}
\end{proof}

\subsection{Proof of \textbf{Theorem~\ref{thm:3}}}

\begin{proof}
Note that
\begin{align*}
  \text{tr}({V})=\text{tr}(\mathcal{J}_X^{-1}{V}_c\mathcal{J}_X^{-1})
  &=\frac{1}{n^2r}\sum_{i=1}^n\text{tr}\left[\frac{1}{\pi_i}
    \{y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^2\mathcal{J}_X^{-1}\bm x_i\bm x_i^T\mathcal{J}_X^{-1}\right]\\
  &=\frac{1}{n^2r}\sum_{i=1}^n\left[\frac{1}{\pi_i}\{y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^2
    \|\mathcal{J}_X^{-1}\bm x_i\|^2\right]\\
  &=\frac{1}{n^2r}\sum_{i=1}^n\pi_i\sum_{i=1}^n\left[\pi_i^{-1}\{y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^2
    \|\mathcal{J}_X^{-1}\bm x_i\|^2\right]\\
  &\ge\frac{1}{n^2r}\left[\sum_{i=1}^n|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\|\mathcal{J}_X^{-1}\bm x_i\|\right]^2,
\end{align*}
where the last inequality follows by the Cauchy-Schwarz inequality, and the equality in it holds if and only if $\pi_i\propto |y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\|\mathcal{J}_X^{-1}\bm x_i\|I\{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\|\mathcal{J}_X^{-1}\bm x_i\|>0\}$.
 Here we define $0/0=0$, and this is equivalent to removing data points with $|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|=0$ in the expression of $V_c$.

\end{proof}

\subsection{Proof of \textbf{Theorems~\ref{thm:r-size} and~\ref{thm:r-size1}}}


 Let ${\left\| A \right\|_F}: =  (\sum_{i = 1}^m {\sum_{j = 1}^n {A_{ij}^2} })^{1/2}$ denote the Frobenius norm, and $\|A\|_S:=  {\sup}_{\bm x \ne 0, \bm x\in \mathbb{R}^n}$ $ {{{{\left\| {A\bm x} \right\|}}}}/{{{{\left\| \bm x \right\|}}}} = {\sigma _{\max }}(A)$ denote the spectral norm.
For a given $m \times n$ matrix $A$ and an $n \times p$ matrix $B$, we want to get an approximation to the
product $AB$. In the following FAST MONTE CARLO Algorithm in \cite{drineas2006fast}, we do $r$ independent trials, in each trial we randomly sample an element of $\{ 1,2, \cdots ,n\} $ with given discrete distribution $P=:\{ {p_i}\} _{i = 1}^n$. Then we extract an $m \times r$ matrix $C$ from the columns of $A$, and extract an $r \times n$ matrix $R$ from the corresponding rows of $B$. If the $P$ is chosen appropriately in the sense that $CR$ is a nice approximation to $AB$,
then the F-norm matrix concentration inequality in Lemma \ref{lem-fMC} holds with high probability.

\begin{lemma}\label{lem-fMC} (Theorem 2.1 in \cite{drineas2006sampling}) Let $A^{(i)}$ be the $i$-th row of $A \in {R^{m \times n}}$ as row vector and $B_{(j)}$ be the $j$-th column of $B \in {R^{n \times p}}$ as column vector. Suppose sampling probabilities $\{ {p_i}\} _{i = 1}^n,(\sum\nolimits_{i = 1}^n {{p_i} = 1} )$ are such that
\[{p_i} \ge \beta \frac{{\left\| {{A^{(i)}}} \right\|\left\| {{B_{(j)}}} \right\|}}{{\sum\nolimits_{j = 1}^n {\left\| {{A^{(i)}}} \right\|\left\| {{B_{(j)}}} \right\|} }}\]
for some $\beta  \in (0,1]$. Construct $C$ and $R$ with Algorithm 1 in \cite{drineas2006fast}, and assume that $\varepsilon  \in (0,1/3)$. Then, with probability at least $1 - \varepsilon $, we have
\[{\left\| {AB - CR} \right\|_F} \le \frac{{4\sqrt {\log (1/\varepsilon )} }}{{\beta \sqrt c }}{\left\| A \right\|_F}{\left\| B \right\|_F}.\]
\end{lemma}

Now we start to prove Theorems ~\ref{thm:r-size} and \ref{thm:r-size1} by applying the above Lemma \ref{lem-fMC}.

\begin{proof}
Note the fact that the maximum likelihood estimate $\hat{\bm \beta}_{\rm MLE}$ of the parameter vector $\bm\beta$ satisfy the following estimation equation
\begin{equation}\label{eq:equation1}
 \bm X_d^T[\bm y-\dot{\psi}(\bm X_d^T{\bm \beta})]=\bm 0,
\end{equation}
 where $\dot{\psi}(\bm X_d^T{\bm \beta})$ denotes the vector whose the $i$-th element is $\dot{\psi}(\bm x_i^T{\bm \beta})$.

 Without of loss of generality, we only show the case with probability $\bm\pi^{\rm mV}$, since the proof for $\bm\pi^{\rm mVc}$ is quite similar.
 Let $S$ be an $n\times r$ matrix whose $i$-th column is $1/\sqrt{r{\pi}^{\mathrm{mV}}_{j_i}}\bm e_{j_i}$, where $\bm e_{j_i}\in \mathbb{R}^n$ denotes the all-zeros vector except that its $j_i$-th entry is set to one. Here $j_i$ denotes the $j_i$-th data point chosen from the $i$-th independent random subsampling with probabilities $\bm\pi^{\mathrm{mV}}$.
 Then $\tilde{\bm\beta}$ satisfies the following equation
\begin{equation}\label{eq:equation2}
 \bm X_d^TSS^T[\bm y-\dot{\psi}(\bm X_d^T{\bm \beta})]=\bm 0.
\end{equation}
 Let $\|\cdot\|_F$ denote the Frobenius norm, we have
 \begin{align*}
  &\quad \sigma_{\min}(\bm X_d^T SS^T)\|\dot{\psi}(\bm X_d^T\tilde{\bm\beta})-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})\|\\
  
  &\le\|{\bm X_d}^TSS^T[\dot{\psi}(\bm X_d^T\tilde{\bm\beta})-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})]\|_F\\
  &\le\|{\bm X_d}^TSS^T[\dot{\psi}(\bm X_d^T\tilde{\bm\beta})-\bm y]\|_F+\|X^TSS^T[\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})]\|_F\\
  &=\|{\bm X_d}^TSS^T[\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})]\|_F~~~~[by~(\ref{eq:equation2})]\\
  &\le \|\bm X_d^T[\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})]\|_F+\|\bm X_d^T[\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})]-{\bm X_d}^TSS^T[\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})]\|_F\\
  &=\|\bm X_d^T[\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})]-{\bm X_d}^TSS^T[\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})]\|_F~~~~[by~ (\ref{eq:equation1})]\\
  &\le \frac{4\kappa(\mathcal{J}_X^{-1})\sqrt{\log(1/\epsilon)}}{\sqrt{r}}\|\bm X_d\|_F\|\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})\|\\
  &\le \frac{4\kappa(\mathcal{J}_X^{-1})\sqrt{\log(1/\epsilon)}}{\sqrt{r}}\sigma_{\max}(\bm X_d)\sqrt{p}\|\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})\|,
\end{align*}
 where the second last inequality follows from Lemma \ref{lem-fMC} by putting $A=X_d^T, B=\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE}), C={\bm X_d}^TS, R=S^T(\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE}))$ and $\beta=1/\kappa(\mathcal{J}_X^{-1})$. Hence,
 \begin{equation}\label{eq:mid-r-size}
   \|\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})-\dot{\psi}(\bm X_d^T\tilde{\bm\beta})\|\le \frac{4\sigma_{\max}(\bm X_d)\sqrt{p\log(1/\epsilon)}}{\sigma_{\min}(\bm X_d^T SS^T)\kappa(\mathcal{J}_X^{-1})\sqrt{r}}\|\bm y-\dot{\psi}(\bm X_d^T\hat{\bm\beta}_{\rm MLE})\|.
 \end{equation}


 Then by following the facts that $\sigma_{\min}(\bm X_d^TSS^T\bm X_d)\le \sigma_{\max}(\bm X_d)\sigma_{\min}(\bm X_d^TSS^T)$ and $\sigma_{\min}^2(\bm{\tilde{X}_d})=\sigma_{\min}(\bm X_d^TSS^T\bm X_d)\ge 0.5\sigma^2_{\min}(\bm X_d)$, it holds that
 \begin{equation}\label{eq:Xss}
 \sigma_{\min}(\bm X_d^TSS^T)\ge 0.5\sigma^2_{\min}(\bm X_d)/\sigma_{\max}(\bm X_d).
 \end{equation}


Now, we turn to prove Theorem \ref{thm:r-size1}.

Similarly, assume that $\|\bm X_d^T\bm X_d-\bm X_d^TSS^T\bm X_d\|_{S} \le {c_d}\|\bm X_d^T\bm X_d-\bm X_d^TSS^T\bm X_d\|_F$ with ${c_d}\le 1$.
Then we have
 \begin{align*}
 |\sigma_{\min}(\bm X_d^T\bm X_d)-\sigma_{\min}(\bm X_d^TSS^T\bm X_d)|&\le \|\bm X_d^T\bm X_d-\bm X_d^TSS^T\bm X_d\|_{S}\\
 &\le {c_d} \|\bm X_d^T\bm X_d-\bm X_d^TSS^T\bm X_d\|_F\\
 &\le {c_d} \frac{4\kappa(\mathcal{J}_X^{-1})\sqrt{\log(1/\epsilon)}}{\sqrt{r}}\|\bm X_d\|_F^2\\
 &\le {c_d} \frac{4\kappa(\mathcal{J}_X^{-1})\sqrt{\log(1/\epsilon)}}{\sqrt{r}}p\sigma_{\max}^2(\bm X_d).
 \end{align*}

Using 
 the above inequality,  the following equation holds with probability at least $1-\varepsilon$:
 \[ \sigma_{\min}(\bm X_d^TSS^T\bm X_d)\ge 0.5\sigma^2_{\min}(\bm X_d),\]
if we set $r>64{c_d^2}\log(1/\epsilon)\sigma_{\max}^4(\bm X_d)\kappa^2(\mathcal{J}_X^{-1})p^2/(\sigma_{\min}^4(\bm X_d))$.

 \end{proof}

\subsection{Proof of \textbf{Theorem~\ref{thm:asy-2step-alg}}}

For the average weighted log-likelihood in Step 2 of Algorithm 2, we have
\begin{align*}
L_{{\tilde{\bm\beta}_0}}^{\rm{two-step}}(\bm\beta ): &= \frac{1}{{r + {r_0}}}\sum\limits_{i = 1}^{r + {r_0}} {\frac{{t_i^*(\bm\beta )}}{{\pi _i^*({\tilde{\bm\beta}_0})}}}  = \frac{1}{{r + {r_0}}}[\sum\limits_{i = 1}^{{r_0}} {\frac{{t_i^*(\bm\beta )}}{{\pi _i^*({\tilde{\bm\beta}_0})}} + \sum\limits_{i = {r_0} + 1}^{r + {r_0}} {\frac{{t_i^*(\bm\beta )}}{{\pi _i^*({\tilde{\bm\beta}_0})}}} } ]\\
&=\frac{{ {{r_0}} }}{{r + {r_0}}} \cdot \frac{1}{{ {{r_0}} }}\sum\limits_{i = 1}^{{r_0}} {\frac{{t_i^*(\bm\beta )}}{{\pi _i^*({\tilde{\bm\beta}_0})}} + \frac{{ r }}{{r + {r_0}}} \cdot \frac{1}{{ r }}\sum\limits_{i = {r_0} + 1}^{r + {r_0}} {\frac{{t_i^*(\bm\beta )}}{{\pi _i^*({\tilde{\bm\beta}_0})}}} },
\end{align*}
where $\pi _i^*({\tilde{\bm\beta}_0})$ in the first item stands for the initial subsampling strategy which satisfies (H.5).


For the sake of brevity, we begin with  the case with probability $\bm\pi^{\rm mV}$.
Denote the log-likelihood in the first and second steps by
\begin{equation*}
L^{*0}_{\tilde{\bm\beta}_0}(\bm\beta)=\frac{1}{{r_0}}\sum_{i=1}^{{r_0}} \frac{t_i^*(\bm\beta)}{\pi_i^*(\tilde{\bm\beta}_0)},\text{\quad and\quad}L^*_{\tilde{\bm\beta}_0}(\bm\beta)=\frac{1}{{r}}\sum_{i=1}^{{r}} \frac{t_i^*(\bm\beta)}{\pi_i^*(\tilde{\bm\beta}_0)},
\end{equation*}
respectively, where $\pi_i(\tilde{\bm\beta}_0)=\check{\pi}_i^{\mathrm{mVc}}$ in $L^*_{\tilde{\bm\beta}_0}(\bm\beta)$, $\check{\pi}_i^{\mathrm{mVc}}$ has the same expression as $\pi_i^{\mathrm{mVc}}$ except that $\hat{\bm\beta}_{\rm MLE}$ is replaced by $\tilde{\bm\beta}_0$. 
 
 
 To proof of Theorem~\ref{thm:asy-2step-alg}, we begin with the following Lemma \ref{lem:lem3}.


\begin{lemma}\label{lem:lem3}
If Assumptions~(H.1)--(H.4) holds, then as $n\rightarrow\infty$, conditionally on $\mathcal{F}_n$ in probability,
\begin{align}
  \check{\mathcal{J}}_X^{\tilde{\bm\beta}_0}-\mathcal{J}_X&=O_{P|\mathcal{F}_n}({r}^{-1/2}),\label{eq:4}\\
   \frac{1}{n}\frac{\partial L^*_{\tilde{\bm\beta}_0}(\hat{\bm\beta}_{\rm MLE})}{\partial\bm\beta} &=O_{P|\mathcal{F}_n}({r}^{-1/2}),\label{eq:28}
\end{align}
where
\begin{equation*}
  \check{\mathcal{J}}_X^{\tilde{\bm\beta}_0}=
  \frac{1}{n}\frac{\partial^2 L^*_{\tilde{\bm\beta}_0}(\hat{\bm\beta}_{\rm MLE})}
  {\partial\bm\beta\partial\bm\beta^T}
  =\frac{1}{n{r}}\sum_{i=1}^{{r}}\frac{\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i^*)\bm x_i^*(\bm x_i^*)^T}{\pi^*_i(\tilde{\bm\beta}_0)}.
\end{equation*}
\end{lemma}
\begin{proof}
Using the same arguments in Lemma \ref{lem:lem1}, we have
\begin{equation*}
  {\rm E}\left(\check{\mathcal{J}}_X^{\tilde{\bm\beta}_0,j_1j_2}-\mathcal{J}_X^{j_1j_2}\Big|\mathcal{F}_n,\tilde{\bm\beta}_0\right)^2
  \le\frac{1}{r}\left(\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^2}{n\pi_i(\tilde{\bm\beta}_0)}\right)\sum_{i=1}^n\frac{(\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i))^2}{n} -\frac{1}{r}(\mathcal{J}_X^{j_1j_2})^2.
\end{equation*}


Direct calculation yields
\begin{equation*}
\begin{split}
  E\left(\check{\mathcal{J}}_X^{\tilde{\bm\beta}_0,j_1j_2}-\mathcal{J}_X^{j_1j_2}\Big|\mathcal{F}_n\right)^2
  &=E_{\tilde{\bm\beta}_0}\left\{\left(\check{\mathcal{J}}_X^{\tilde{\bm\beta}_0,j_1j_2}-\mathcal{J}_X^{j_1j_2}\Big|\mathcal{F}_n,\tilde{\bm\beta}_0\right)^2\right\}\\
  &\le \frac{1}{r}\left(\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^2}{n\pi_i(\tilde{\bm\beta}_0)}\right)\sum_{i=1}^n\frac{(\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i))^2}{n} +\frac{1}{r}(\mathcal{J}_X^{j_1j_2})^2,\\
  
  \end{split}
\end{equation*}
where $E_{\tilde{\bm\beta}_0}$ means that the expectation is taken with respect to the distribution of $\tilde{\bm\beta}_0$ given $\mathcal{F}_n$. 

From Assumptions (H.1) and (H.5), it be shown that
\[\begin{split}
\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^2}{n\pi_i(\tilde{\bm\beta}_0)}&=\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^2\sum_{j=1}^n\max(|y_j-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_j)|,{\delta})\|\bm x_j\|}{n\max(|y_i-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)|,{\delta})\|\bm x_i\|}\\
&\le \max_{i=1,\ldots,n}\frac{\|\bm x_i\|\sum_{j=1}^n(|y_j-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_j)|+{\delta})\|\bm x_j\|}{n{\delta}}\\
&\le \frac{C^2\sum_{j=1}^n(|y_j-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_j)|+{\delta})}{n{\delta}}\\
&\le C^2+\frac{C^2\sum_{j=1}^n |y_j-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_j)|}{n{\delta}}\\
&= O_P(1),
\end{split}
\]
where $C$ in the second last inequality is the upper bound of $\|\bm x_i\|$ due to (H.1), and the last equality holds from Lemma~\ref{lem-moment} and Chebyshev's inequality.

It worths to mention that $\sum_{i=1}^n (\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i))^2/n$ is bounded by noting the facts that the parameter space is compact, $\hat{\bm\beta}_{\rm MLE}$ is the unique global maximum of log-likelihood function, and Assumption (H.1). Thus
\begin{equation*}
  E\left(\check{\mathcal{J}}_X^{\tilde{\bm\beta}_0,j_1j_2}-\mathcal{J}_X^{j_1j_2}\Big|\mathcal{F}_n\right)^2=O_P(r^{-1}).
\end{equation*}


On the other hand, following the same arguments in Lemma \ref{lem:lem1}, we can have
\[
E\left\{\frac{L^*_{\tilde{\bm\beta}_0}(\bm\beta)}{n}-\frac{L(\bm\beta)}{n}\bigg|\mathcal{F}_n,\tilde{\bm\beta}_0\right\}^2
  =O_P(r^{-1}).
\]

Then $E\left\{n^{-1}{L^*_{\tilde{\bm\beta}_0}(\bm\beta)}-n^{-1}{L(\bm\beta)}|\mathcal{F}_n\right\}^2=O_P(r^{-1})$, the desired result holds.
\end{proof}


Now we prove Theorem~\ref{thm:asy-2step-alg}.
\begin{proof}
Using the same arguments in Lemma \ref{lem:lem3} and Slutsky's theorem, we have
\[\begin{split}
E\left\{\frac{L^{\rm{two-step}}_{\tilde{\bm\beta}_0}(\bm\beta)}{n}-\frac{L(\bm\beta)}{n}\bigg|\mathcal{F}_n\right\}^2=&(\frac{{ {{r_0}} }}{{r + {r_0}}})^2E\left\{\frac{L^{*0}_{\tilde{\bm\beta}_0}(\bm\beta)}{n}-\frac{L(\bm\beta)}{n}\bigg|\mathcal{F}_n\right\}^2\\
&+(\frac{{ {{r}} }}{{r + {r_0}}})^2E\left\{\frac{L^{*}_{\tilde{\bm\beta}_0}(\bm\beta)}{n}-\frac{L(\bm\beta)}{n}\bigg|\mathcal{F}_n\right\}^2
  =O_P(r^{-1}).
  \end{split}
\]

Therefore $E\{n^{-1}{L^{\rm{two-step}}_{\tilde{\bm\beta}_0}(\bm\beta)}-n^{-1}{L(\bm\beta)}|\mathcal{F}_n\}^2\rightarrow 0$ as $r_0/r\rightarrow 0,r\rightarrow\infty$ and $n^{-1}{L^{\rm{two-step}}_{\tilde{\bm\beta}_0}(\bm\beta)}-n^{-1}{L(\bm\beta)}\rightarrow 0$ in conditional probability given $\mathcal{F}_n$. Also note that the parameter space is compact and $\hat{\bm\beta}_{\rm MLE}$ is the unique global maximum of the continuous convex function $L(\bm\beta)$. Thus, from Theorem 5.9 and its remark of \cite{Vaart2000Asymptotic}, we have
\[
  \|{\breve{\bm\beta}}-\hat{\bm\beta}_{\rm MLE}\|=o_{P|\mathcal{F}_n}(1).
\]
 Using Taylor's theorem,
 \begin{align*}
0&=\frac{\dot{L}^{\rm{two-step}}_{\tilde{\bm\beta_0},j}(\breve{\bm\beta})}{n}=\frac{{ {{r_0}} }}{{r + {r_0}}}\frac{\dot{L}^{*0}_{\tilde{\bm\beta_0},j}(\breve{\bm\beta})}{n}+\frac{{ {{r}} }}{{r + {r_0}}}\frac{\dot{L}^{*}_{\tilde{\bm\beta_0},j}(\breve{\bm\beta})}{n}\\
&=\frac{{ {{r}} }}{{r + {r_0}}}\left\{\frac{\dot{L}^*_{\tilde{\bm\beta_0},j}(\hat{\bm\beta}_{\rm MLE})}{n}
   +\frac{1}{n} \frac{\partial\dot{L}^*_{\tilde{\bm\beta_0},j}(\hat{\bm\beta}_{\rm MLE})}{\partial \bm\beta^T}(\breve{\bm\beta}-\hat{\bm\beta}_{\rm MLE})
   +\frac{1}{n} R_{\tilde{\bm\beta_0},j}\right\}+\frac{{ {{r_0}} }}{{r + {r_0}}}\frac{\dot{L}^{*0}_{\tilde{\bm\beta_0},j}(\breve{\bm\beta})}{n},
\end{align*}
where $\dot{L}^*_{\tilde{\bm\beta_0},j}({\bm\beta})$ is the partial derivative of $L^*_{\tilde{\bm\beta_0},j}({\bm\beta})$ with respect to $\beta_j$, and the Lagrange remainder
\begin{equation*}
  \frac{1}{n}R_{\tilde{\bm\beta_0},j}=\frac{1}{n}(\breve{\bm\beta}-\hat{\bm\beta}_{\rm MLE})^T
   \int_0^1\int_0^1\frac{\partial^2\dot{L}^*_j\{\hat{\bm\beta}_{\rm MLE}
   +uv(\breve{\bm\beta}-\hat{\bm\beta}_{\rm MLE})\}}{\partial \bm\beta\partial
   \bm\beta^T}v du dv\
(\breve{\bm\beta}-\hat{\bm\beta}_{\rm MLE}).
\end{equation*}
From Assumptions (H.1), (H.3) and (H.5),
\begin{align*}
  \left\|\frac{\partial^2\dot{L}^*_j(\bm\beta)}
  {\partial\bm\beta\partial\bm\beta^T}\right\|=
  &\frac{1}{r}\left\|\sum_{i=1}^r \frac{-\dddot{\psi}(\bm\beta^T\bm x_i)}{\pi^*_i(\tilde{\bm\beta}_0)}x^*_{ij}\bm x^*_i{\bm x^*_i}^T\right\|
    \le\frac{C_3}{r}\sum_{i=1}^r\frac{\|\bm x^*_i\|^3}{\pi^*_i(\tilde{\bm\beta}_0)}
\end{align*}
for all $\bm\beta$, where $C_3$ is the upper bound of $\|\dddot{\psi}(\bm\beta^T\bm x_i)\|$. From Remark~\ref{rm:analytic} and (H1) -(H.3), it is known that $C_3$ is a constant. Thus
\begin{align*}
  \frac{1}{n}\left\|\int_0^1\int_0^1\frac{\partial^2\dot{L}^*_j\{\hat{\bm\beta}_{\rm MLE}
   +uv(\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE})\}}{\partial \bm\beta\partial\bm\beta^T}vdudv\ \right\|
  \le\frac{C_3}{2r}\sum_{i=1}^r\frac{\|\bm x^*_i\|^3}{n\pi^*_i(\tilde{\bm\beta}_0)}=O_{P|\mathcal{F}_n}(1),
\end{align*}
where the last equality is from (H.1), (H.5) and the fact that
\[P\left(\frac{1}{nr}\sum_{i=1}^r \frac{\|\bm x^*_i\|^3}{\pi^*_i(\tilde{\bm\beta}_0)}\ge \tau\Bigg{|}\mathcal{F}_n\right)\le\frac{1}{nr\tau}\sum_{i=1}^r E\left(\frac{\|\bm x^*_i\|^3}{\pi^*_i(\tilde{\bm\beta}_0)}\Bigg|\mathcal{F}_n\right)=\frac{1}{n\tau}\sum_{i=1}^n\|\bm x_i\|^3\rightarrow 0,\]
as $\tau\rightarrow\infty$.

Since $r_0/r\rightarrow 0$, it is easy to see that
\[\frac{{ {{r_0}} }}{{r + {r_0}}}\frac{\dot{L}^{*0}_{\tilde{\bm\beta_0},j}(\breve{\bm\beta})}{n}=o_{P|\mathcal{F}_n}(r^{-1/2}),\]
and $r/(r+r_0)\rightarrow1$.
Hence,
\begin{equation}\label{eq:13}
  \breve{\bm\beta}-\hat{\bm\beta}_{\rm MLE}=
  -(\check{\mathcal{J}}^{\tilde{\bm\beta}_0}_X)^{-1}\left\{\frac{\dot {L}^*_{\tilde{\bm\beta}_0}(\hat{\bm\beta}_{\rm MLE})}{n}
    +O_{P|\mathcal{F}_n}(\|\breve{\bm\beta}-\hat{\bm\beta}_{\rm MLE}\|^2)\right\},
\end{equation}
as $r_0/r\rightarrow0$.
From  Lemma~\ref{lem:lem3} and (H.4), $(\tilde{\mathcal{J}}^{\tilde{\bm\beta}_0}_X)^{-1} =O_{P|\mathcal{F}_n}(1)$. Then it follows from Slutsky's theorem that
\begin{equation*}
  \breve{\bm\beta}-\hat{\bm\beta}_{\rm MLE}
  =O_{P|\mathcal{F}_n}(r^{-1/2})+o_{P|\mathcal{F}_n}(\|\tilde{\bm\beta}-\hat{\bm\beta}_{\rm MLE}\|),
\end{equation*}
which implies that
\[  \breve{\bm\beta}-\hat{\bm\beta}_{\rm MLE}=O_{P|\mathcal{F}_n}(r^{-1/2}).\]

 For the case $\pi_i(\tilde{\bm\beta}_0)=\check{\pi}_i^{\mathrm{mV}}$ in $L^*_{\tilde{\bm\beta}_0}(\bm\beta)$ with $\check{\pi}_i^{\mathrm{mV}}$ has the same expression as $\pi_i^{\mathrm{mV}}$ except that $\hat{\bm\beta}_{\rm MLE}, \mathcal{J}_X$ is replaced by $\tilde{\bm\beta}_0$ and $\tilde{\mathcal{J}}_X$, respectively, we have $\pi_i(\tilde{\bm\beta}_0)\ge \kappa(\tilde{\mathcal{J}}_X)^{-1}\check{\pi}_i^{\mathrm{mVc}}$.
Under Assumptions (H.1) and (H.4), we know $\kappa({\mathcal{J}}_X)$ is bounded. Thus, as $r_0\rightarrow\infty$,  $\kappa(\tilde{\mathcal{J}}_X)$ is also bounded.
Therefore, the rest of the proof is the same as that of $\check{\pi}_i^{\mathrm{mVc}}$ with minor modifications.
\end{proof}

\subsection{Proof of \textbf{Theorem~\ref{thm:CLT-2step}}}
\begin{proof}
For the sake of brevity, we begin with  the case with probability $\check{\bm\pi}^{\rm mV}$.
Denote
\begin{equation}\label{eq:56}
  \frac{\dot L^*_{\tilde{\bm\beta}_0}(\hat{\bm\beta}_{\rm MLE})}{n}
  =\frac{1}{{r}}\sum_{i=1}^{{r}}
  \frac{\{y^*_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i^*)\}\bm x^*_i}
  {n{\pi}^*_i(\tilde{\bm\beta}_0)}
  =:\frac{1}{{r}}\sum_{i=1}^{{r}}\eta_i^{\tilde{\bm\beta}_0}.
\end{equation}

It can be shown that given $\mathcal{F}_n$ and ${\tilde{\bm\beta}_0}$, $\eta_1^{\tilde{\bm\beta}_0}, \ldots, \eta_{{r}}^{\tilde{\bm\beta}_0}$ are i.i.d
random variables with zero mean and variance
\begin{align*}
  &\text{var}(\eta_i|\mathcal{F}_n,\tilde{\bm\beta}_0)
    ={r}V_c^{\tilde{\bm\beta}_0}
    =\frac{1}{n^2}\sum_{i=1}^n \pi_i(\tilde{\bm\beta}_0)\frac{\{y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^2
    \bm x_i\bm x_i^T}{{\pi}^2_i(\tilde{\bm\beta}_0)}.
    
\end{align*}
Meanwhile, for every $\varepsilon>0$,
\begin{align*}
   \sum_{i=1}^{{r}} E\{\|{r}^{-1/2}\eta_i^{\tilde{\bm\beta}_0}\|^2
    I(\|\eta_i^{\tilde{\bm\beta}_0}\|>{r}^{1/2}\varepsilon)
    |\mathcal{F}_n,\tilde{\bm\beta}_0\} & \le\frac{1}{{r}^{3/2}\varepsilon}
    \sum_{i=1}^{{r}} E\{\|\eta_i^{\tilde{\bm\beta}_0}\|^{3}
    I(\|\eta_i^{\tilde{\bm\beta}_0}\|>{r}^{1/2}\varepsilon)|
    \mathcal{F}_n,\tilde{\bm\beta}_0\} \\
   &
    \le\frac{1}{{r}^{3/2}\varepsilon}
    \sum_{i=1}^{{r}} E(\|\eta_i^{\tilde{\bm\beta}_0}\|^{3}|
    \mathcal{F}_n,\tilde{\bm\beta}_0)\\
  &\le\frac{1}{{r}^{1/2}}\frac{1}{n^{3}}
    \sum_{i=1}^n\frac{\{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\}^{3}
    \|\bm x_i\|^{3}}{\pi_i^{2}(\tilde{\bm\beta}_0)} \\
   &
    \le\frac{1}{{r}^{1/2}}\left(\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^{3}}{n^{2}\pi_i^{2}(\tilde{\bm\beta}_0)}\right)
    \sum_{i=1}^n\frac{\{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\}^{3}}
    {n}.
\end{align*}

Combining (H.1), (H.5), Lemma~\ref{lem-moment} and Chebyshev's inequality, it can be proved that
\[0\le \sum_{i=1}^n\frac{\{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\}^{3}}{n}=O_P(1),\]
and
\[\begin{split}
0\le \max_{i=1,\ldots,n}\frac{\|\bm x_i\|^3}{n^2\pi_i^2(\tilde{\bm\beta}_0)}&=\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^3\left\{\sum_{j=1}^n\max(|y_j-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_j)|,{\delta})\|\bm x_j\|\right\}^2}{n^2\left\{\max(|y_i-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_i)|,{\delta})\|\bm x_i\|\right\}^2}\\
&\le \max_{i=1,\ldots,n}\frac{\|\bm x_i\|\left\{\sum_{j=1}^n\max(|y_j-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_j)|,{\delta})\|\bm x_j\|\right\}^2}{n^2{\delta}}\\
&\le \frac{C^3[\sum_{j=1}^n(|y_j-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_j)|+{\delta})]^2}{n^2{\delta}}\\
&\le \frac{C^3}{\delta}\left(\frac{\sum_{j=1}^n |y_j-\dot{\psi}(\tilde{\bm\beta}_0^T\bm x_j)|}{n}+{\delta}\right)^2\\
&= O_P(1),
\end{split}
\]
where $C$ in second last inequality is the upper bound of $\|\bm x_i\|$ due to (H.1), and last equality is from Lemma~\ref{lem-moment} and Chebyshev's inequality. Hence
\[\sum_{i=1}^{{r}} E\{\|{r}^{-1/2}\eta_i^{\tilde{\bm\beta}_0}\|^2
    I(\|\eta_i^{\tilde{\bm\beta}_0}\|>{r}^{1/2}\varepsilon)
    |\mathcal{F}_n,\tilde{\bm\beta}_0\}=o_P(1).\]
This shows that the Lindeberg-Feller conditions are satisfied in probability.
By the Lindeberg-Feller central limit theorem \citep[Proposition 2.27 of][]{Vaart2000Asymptotic}, conditionally on $\mathcal{F}_n$ and $\tilde{\bm\beta}_0$,
\begin{equation*}
  \frac{1}{n}(V_c^{\tilde{\bm\beta}_0})^{-1/2}\dot L^*(\hat{\bm\beta}_{\rm MLE})=
  \frac{1}{{r}^{1/2}}\{\text{var}(\eta_i|\mathcal{F}_n,\tilde{\bm\beta}_0)\}^{-1/2}\sum_{i=1}^{{r}}\eta_i
  \rightarrow N(0,I),
\end{equation*}
in distribution.


By the same approach of the proof in Lemma \ref{lem:lem3}, we have
$$
\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^2}{n\pi_i(\tilde{\bm\beta}_0)}= O_P(1),
$$

Similarly,
\[\max_{i=1,\ldots,n}\frac{\|\bm x_i\|^2}{n\pi_i}= O_P(1).\]
From (H.1) and (H.5), the distance between $V_c^{\tilde{\bm\beta}_0}$ and $V_c$ is
\[\|V_c-V_c^{\tilde{\bm\beta}_0}\|\le \frac{1}{r}\left\{\max_{i=1,\ldots,n}\left(\frac{\|\bm x_i\|^2}{n\pi_i}\right)+\max_{i=1,\ldots,n}\left(\frac{\|\bm x_i\|^2}{n\pi_i(\tilde{\bm\beta}_0)}\right)\right\}\sum_{i=1}^n\frac{\{y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^2}{n}
  =O_{P|\mathcal{F}_n}({r}^{-1}),\]
  where the last equation follows from  the facts that
  \[0\le \max_{i=1,\ldots,n}\left(\frac{\|\bm x_i\|^2}{n\pi_i}\right)+\max_{i=1,\ldots,n}\left(\frac{\|\bm x_i\|^2}{n\pi_i(\tilde{\bm\beta}_0)}\right)=O_P(1),\]
  and
  \[0\le \sum_{i=1}^n\frac{\{y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)\}^2}{n}=O_P(1).\]
Here the last equality holds from Lemma~\ref{lem-moment} and Chebyshev's inequality.


From Lemma~\ref{lem:lem3} and~Theorem \ref{thm:asy-2step-alg},
\begin{equation}\label{eq:57}
  {\breve{\bm\beta}}-\hat{\bm\beta}_{\rm MLE}=
  -\frac{1}{n}(\check{\mathcal{J}}_X^{\tilde{\bm\beta}_0})^{-1}
  \dot L^*_{\tilde{\bm\beta}_0}(\hat{\bm\beta}_{\rm MLE})+O_{P|\mathcal{F}_n}({r}^{-1})
\end{equation}
Based on Equation (\ref{eq:4}), we further have
\begin{align*}
  (\check{\mathcal{J}}_X^{\tilde{\bm\beta}_0})^{-1}-\mathcal{J}_X^{-1}
  &=-\mathcal{J}_X^{-1}(\check{\mathcal{J}}_X^{\tilde{\bm\beta}_0}-\mathcal{J}_X)
    (\check{\mathcal{J}}_X^{\tilde{\bm\beta}_0})^{-1}
  =O_{P|\mathcal{F}_n}({r}^{-1/2}).
\end{align*}
Then
\begin{align*}
  &V^{-1/2}({\breve{\bm\beta}}-\hat{\bm\beta}_{\rm MLE})\\
  &=-V^{-1/2}n^{-1}(\check{\mathcal{J}}_X^{\tilde{\bm\beta}_0})^{-1}
    \dot L^*(\hat{\bm\beta}_{\rm MLE})
    +O_{P|\mathcal{F}_n}({r}^{-1/2})\\
  &=-V^{-1/2}\mathcal{J}_X^{-1}n^{-1}\dot L^*(\hat{\bm\beta}_{\rm MLE})
    -V^{-1/2}\{(\check{\mathcal{J}}_X^{\tilde{\bm\beta}_0})^{-1}-\mathcal{J}_X^{-1}\}
    n^{-1}\dot L^*(\hat{\bm\beta}_{\rm MLE})
    +O_{P|\mathcal{F}_n}({r}^{-1/2})\\
  &=-V^{-1/2}\mathcal{J}_X^{-1}(V_c^{\tilde{\bm\beta}_0})^{1/2}
    (V_c^{\tilde{\bm\beta}_0})^{-1/2}n^{-1}\dot L^*(\hat{\bm\beta}_{\rm MLE})
    +O_{P|\mathcal{F}_n}({r}^{-1/2}).
\end{align*}
 It can be shown that
 \begin{equation*}
 \begin{split}
   V^{-1/2}\mathcal{J}_X^{-1}(V_c^{\tilde{\bm\beta}_0})^{1/2}(V^{-1/2}\mathcal{J}_X^{-1}(V_c^{\tilde{\bm\beta}_0})^{1/2})^T
   &=V^{-1/2}\mathcal{J}_X^{-1}(V_c^{\tilde{\bm\beta}_0})\mathcal{J}_X^{-1}V^{-1/2}\\
   &=V^{-1/2}\mathcal{J}_X^{-1}(V_c)\mathcal{J}_X^{-1}V^{-1/2}+O_{P|\mathcal{F}_n}(r^{-1/2}) \\
   &=I+O_{P|\mathcal{F}_n}(r^{-1/2}).
   \end{split}
 \end{equation*}
So the desired result follows by Slutsky's theorem.

 As for the case $\pi_i(\tilde{\bm\beta}_0)=\check{\pi}_i^{\mathrm{mV}}$ in $L^*_{\tilde{\bm\beta}_0}(\bm\beta)$.  $\check{\pi}_i^{\mathrm{mV}}$ has the same expression as $\pi_i^{\mathrm{mV}}$ except that $\hat{\bm\beta}_{\rm MLE}, \mathcal{J}_X$ is replaced by $\tilde{\bm\beta}_0$ and $\tilde{\mathcal{J}}_X$, respectively, we have $\tilde{\mathcal{J}}_X\rightarrow\mathcal{J}_X$ as $r_0\rightarrow\infty$. Thus we can replace $\tilde{\mathcal{J}}_X$ in $\pi_i(\tilde{\bm\beta}_0)$ by $\mathcal{J}_X$ through continue mapping theory and also note that
$\pi_i(\tilde{\bm\beta}_0)\ge \kappa({\mathcal{J}}_X)^{-1}\check{\pi}_i^{\mathrm{mVc}}$. Thus the rest of the proof is the same as that of $\check{\pi}_i^{\mathrm{mVc}}$ with  minor modifications.
\end{proof}


\subsection{Discuss the case of unbound covariates}
\label{sec:discuss-case-unbound}
For the purpose of obtaining the asymptotic results in terms of unbound covariates situation, here we pose some alternative assumptions to replace (H.1) and (H.2). It needs to show that Theorems~\ref{thm:as-general-alg}, \ref{thm:CLT}, \ref{thm:asy-2step-alg} and \ref{thm:CLT-2step} still hold.

 \begin{itemize}
\item [\textbullet] (A.1): $\frac{1}{n^2}\sum^n_{i=1}\frac{\{|y_i-\dot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)|\}^k\bm x_i\bm x_i^T}{\pi_i}=O_P(1)$ for $k=1,2,3$.
\item [\textbullet] (A.2): Restrict our analysis of regression coefficient vector in the following compact domain:
        $$J_B=\left\lbrace  \bm\beta \in \mathbb{R}^{p}:\: \bm\beta^{T}\bm x \in \Theta, \|\bm\beta\|_1 \le B \right\rbrace, $$
        where $B$ is a constant and $\Theta$ is a compact set.
\item [\textbullet] (A.3): If $\bm\beta \in J_B$, then $\frac{1}{n}\sum\limits_{i = 1}^n {{{|{\varphi ^{(l)}}({\bm\beta ^T}{\bm x_i})|}^k}}  = {O_P}(1)$ for $k=1,2,3$ and $l=0,1,2,3$.
\item [\textbullet] (A.4): $\frac{1}{n}\sum\limits_{i = 1}^n {\left\|{\bm x_i}\right\|^3}  = {O_P}(1)$.
\item [\textbullet] (A.5): $\frac{1}{n}\sum\limits_{i = 1}^n {{{\left\| {{y_i}{\bm x}_i} \right\|}^s}}  = {O_P}(1)$ for $s=1,2$.
\item [\textbullet] (A.6): $\frac{1}{n^2}\sum^n_{i=1}
     \frac{[\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)x_{ij_1} x_{ij_2}]^2}{\pi_i}={O_P}(1)$ for $1\le j_1,j_2\le p$.
\end{itemize}


The aforementioned three assumptions are essentially moment conditions and are very general. It is clearly to see that  (H.1) and (H.2) are the sufficient conditions for (A.1)--(A.10).

I. \textbf{Theorem~\ref{thm:as-general-alg} for unbound covariates}:

Go on with the line of proving {Theorem~\ref{thm:as-general-alg}}, we only need to check (\ref{eq:J-J}), (\ref{eq:l2error}) and (\ref{eq:53}).

For (\ref{eq:J-J}), (A.6) gives
\begin{align*}
  {\rm E}\left(\tilde{\mathcal{J}}_X^{j_1j_2}-\mathcal{J}_X^{j_1j_2}\Big|\mathcal{F}_n\right)^2
\le &\frac{1}{rn^2}\sum^n_{i=1}
     \frac{[\ddot{\psi}(\hat{\bm\beta}_{\rm MLE}^T\bm x_i)x_{ij_1} x_{ij_2}]^2}{\pi_i}
     -\frac{1}{r}(\mathcal{J}_X^{j_1j_2})^2\\
      &={O_P}({r^{ - 1}}).
\end{align*}

For (\ref{eq:l2error}), observe that
\[t_i^{}(\bm\beta ) \le |{y_i}{\bm\beta ^T}{\bm x_i} - \psi ({\bm\beta ^T}{\bm x_i})| \le B\left\| {{y_i}{\bm x_i}} \right\| + |\psi ({\bm\beta ^T}{\bm x_i})|\]
which implies by (A.5) with $s=1$ and (A.3) with $k=1$
\[\frac{1}{n}\sum\limits_{i = 1}^n {t_i^{}(\bm\beta )}  \le \frac{B}{n}\sum\limits_{i = 1}^n {\left\| {{y_i}{\bm x_i}} \right\|}  + \frac{1}{n}\sum\limits_{i = 1}^n {|\psi ({\bm\beta ^T}{\bm x_i})|}  = B{O_P}(1) + {O_P}(1) = {O_P}(1)\]
Moreover,
\[t_i^2(\bm \beta ) = |{y_i}{\bm\beta ^T}{\bm x_i} - \psi ({\bm \beta ^T}{\bm x_i}){|^2} \le {B^2}{\left\| {{y_i}{\bm x_i}} \right\|^2} + 2B\left\| {{y_i}{\bm x_i}} \right\||\psi ({\bm\beta ^T}{\bm x_i})| + |\psi ({\bm\beta ^T}{\bm x_i}){|^2}\]
By Cauchy inequality, (A.5) with $s=2$ and (A.3)  with $k=2$ show that
\begin{align*}
\frac{1}{n}\sum\limits_{i = 1}^n {t_i^2(\bm\beta )}  & \le \frac{{{B^2}}}{n}\sum\limits_{i = 1}^n {{{\left\| {{y_i}{\bm x_i}} \right\|}^2}}  + 2B{(\frac{1}{n}\sum\limits_{i = 1}^n {{{\left\| {{y_i}{\bm x_i}} \right\|}^2}} )^{1/2}}{(\frac{1}{n}\sum\limits_{i = 1}^n {|\psi ({\bm\beta ^T}{\bm x_i}){|^2}} )^{1/2}} + \frac{1}{n}\sum\limits_{i = 1}^n {|\psi ({\bm\beta ^T}{\bm x_i}){|^2}} \\
& \le {B^2}{O_P}(1) + 2B\sqrt {{O_P}(1){O_P}(1)}  + {O_P}(1) = {O_P}(1).
\end{align*}
Thus (\ref{eq:l2error}) is verified.

For (\ref{eq:53}), it is easy to see by (A.4)
\begin{align}\label{eq:53ub}
  P\left(\frac{1}{nr}\sum_{i=1}^r\frac{\|\bm x^*_i\|^3}{\pi^*_i}\ge\tau
  \Bigg|\mathcal{F}_n\right) \le \frac{1}{n\tau}\sum_{i=1}^n\|\bm x_i\|^3 =\frac{1}{\tau }{O_P}(1) \rightarrow0,
\end{align}
as $\tau\rightarrow\infty$.

II. \textbf{\textbf{Theorems~\ref{thm:CLT}}, \ref{thm:asy-2step-alg} and \ref{thm:CLT-2step} for unbounded covariates}:

To prove the other required theorem with unbounded covariates, we deal with more assumptions (A.7)-(A.9) below
 \begin{itemize}
\item [\textbullet] (A.7): If $\bm\beta \in J_B$, then $\frac{1}{n^{2+\delta}}\sum^n_{i=1}\frac{\{y_i-\dot{\psi}_i({\bm\beta}^T\bm x_i)\}^{2+\delta}
    \|\bm x_i\|^{2+\delta}}{\pi_i^{1+\delta}}= {O_P}(1)$ for some $\delta >0$.

\item [\textbullet] (A.8): If $\bm\beta \in J_B$, then $\max_{i=1,\ldots,n}\frac{\|\bm x_i\|\sum_{j=1}^n(|y_j-\dot{\psi}(\bm\beta ^T\bm x_j)|+{\delta})\|\bm x_j\|}{n{\delta}}= {O_P}(1)$ for any positive $\delta $.

\item [\textbullet] (A.9): If $\bm\beta \in J_B$, then $\max_{i=1,\ldots,n}\frac{\|\bm x_i\|\left\{\sum_{j=1}^n\max(|y_j-\dot{\psi}(\bm\beta^T\bm x_j)|,{\delta})\|\bm x_j\|\right\}^2}{n^2{\delta^2}}= {O_P}(1)$.
 \end{itemize}
The derivations are similar to the bounded case, so we omit the proof.
\label{eq:clt}
{"config": "arxiv", "file": "1806.06761/draft.tex"}