Divide-and-Conquer Information-Based Optimal Subdata Selection Algorithm

TL;DR

This paper introduces a divide-and-conquer version of the IBOSS algorithm for selecting informative subdata from large datasets that cannot fit into memory, maintaining efficiency and statistical properties.

Contribution

It develops a divide-and-conquer IBOSS method for large data, providing theoretical guarantees and practical evaluation for memory-limited environments.

Findings

Estimator retains efficiency with a small number of partitions

The method is computationally feasible for very large datasets

Numerical experiments confirm theoretical properties

Abstract

The information-based optimal subdata selection (IBOSS) is a computationally efficient method to select informative data points from large data sets through processing full data by columns. However, when the volume of a data set is too large to be processed in the available memory of a machine, it is infeasible to implement the IBOSS procedure. This paper develops a divide-and-conquer IBOSS approach to solving this problem, in which the full data set is divided into smaller partitions to be loaded into the memory and then subsets of data are selected from each partitions using the IBOSS algorithm. We derive both finite sample properties and asymptotic properties of the resulting estimator. Asymptotic results show that if the full data set is partitioned randomly and the number of partitions is not very large, then the resultant estimator has the same estimation efficiency as the original IBOSS estimator. We also carry out numerical experiments to evaluate the empirical performance of the proposed method.

Figures

Tables

Table 1. Table 1 : Computing times for different number of partitions B 𝐵 B with N = 10 7 𝑁 superscript 10 7 N=10^{7} , p = 10 2 𝑝 superscript 10 2 p=10^{2} , and k = 10 4 𝑘 superscript 10 4 k=10^{4} . “FULL” uses all the observations in each data block; “IBOSS” uses the IBOSS algorithm to select a subdata set of k / B 𝑘 𝐵 k/B observations; and “UNI uses Poisson subsampling to select a subdata set of k / B 𝑘 𝐵 k/B observations.

Methods	$B=2$	$B=5$	$B=10$	$B=25$	$B=50$
FULL	98.80	95.76	90.82	86.41	85.28
IBOSS	29.74	24.90	20.79	18.81	18.42
UNI	6.81	5.67	4.32	6.75	14.08

Equations

$$y_i={\beta }_0+\sum _{j=1}^p{z}_{ij}{\beta }_j+{\epsilon }_i,i=1,...,N,$$

$${\hat{\mathbf{\beta }}}_{\mathrm{full}}={\left(\sum _{i=1}^N{\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{T}}\right)}^{-1}\sum _{i=1}^N{\mathbf{x}}_i{y}_i,$$

$$\mathbb{V}(\hat{\bm{\beta}}_{\mathrm{full}}|\mathcal{Z})=\bm{M}_{\mathrm{full}}^{-1}=\bigg{(}\sum_{i=1}^{N}\bm{x}_{i}\bm{x}_{i}^{\rm T}\bigg{)}^{-1},$$

$${\mathbf{M}}_{\mathbf{\delta }}=\sum _{i=1}^N{\delta }_i{\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{T}}.$$

$${\mathbf{\delta }}_{\mathrm{D}}^{opt}=\underset{\mathbf{\delta }:|\mathbf{\delta }|=k}{\mathrm{argmax}}\left|\sum _{i=1}^N{\delta }_i{\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{T}}\right|.$$

$$|{\mathbf{M}}_{\mathbf{\delta }}|\le {\zeta }_N:=\frac{k^{p+1}}{4^p}\prod _{j=1}^p(z_{(N)j}-z_{(1)j}{)}^2,$$

$$\frac{|({\mathcal{X}}_{{\mathcal{D}}_{\mathrm{BS}}}^{\ast }{)}^{\mathrm{T}}{\mathcal{X}}_{{\mathcal{D}}_{\mathrm{BS}}}^{\ast }|}{{\zeta }_N}\ge \max (C_R,C_E),$$

$$C_R=$$

$$C_E=$$

$$\frac{|({\mathcal{X}}_{{\mathcal{D}}_S}^{\ast }{)}^{\mathrm{T}}{\mathcal{X}}_{{\mathcal{D}}_S}^{\ast }|}{{\zeta }_N}\ge \frac{{\lambda }_{\min }^p({\mathbf{R}}_{{\mathcal{D}}_S})}{p^p}\prod _{j=1}^p{\left(\frac{z_{(N-r+1)j}-z_{(r)j}}{z_{(N)j}-z_{(1)j}}\right)}^2,$$

$${\hat{\mathbf{\beta }}}^{{\mathcal{D}}_{\mathrm{BS}}}={\left(\sum _{b=1}^B\sum _{i=1}^{n_B}{\mathbf{x}}_{bi}^{\ast }{{\mathbf{x}}_{bi}^{\ast }}^{\mathrm{T}}\right)}^{-1}\left(\sum _{b=1}^B\sum _{i=1}^{n_B}{\mathbf{x}}_{bi}^{\ast }{y}_{bi}^{\ast }\right),$$

$$\mathbb{V}({\hat{\beta }}_0^{{\mathcal{D}}_{\mathrm{BS}}}|\mathcal{Z})\ge \frac{{\sigma }^2}{k},\text{ and}$$

$$\frac{4{\sigma }^2}{k(z_{(N)j}-z_{(1)j}{)}^2}\le$$

$$V_{aj}$$

$$V_{ej}$$

$$V_{oj}=\frac{4p{\sigma }^2}{k{\lambda }_{\min }({\mathbf{R}}_{{\mathcal{D}}_S})(z_{(N-r+1)j}-z_{(r)j}{)}^2},j=1,...,p.$$

$$\mathbb{V}({\hat{\beta }}_j^{{\mathcal{D}}_{\mathrm{BS}}}|\mathcal{Z})=O_P\left\{\frac{p}{k(z_{(N)j}-z_{(1)j}{)}^2}\right\},j=1,...,p,$$

$$\frac{r}{N[1-F_j\{(1-ϵ)F_j^{-1}(1-N^{-1})\}]}\to 0\text{ and }\frac{F_j(1-1/N)}{F_j(1-B/N)}\to 1,$$

$$\frac{r}{N{F}_j\{(1-ϵ)F_j^{-1}(N^{-1})\}}\to 0,\text{ and }\frac{F_j(1/N)}{F_j(B/N)}\to 1,$$

$$F_j(z)=1-\exp \{-z^{\gamma }h(z)\},$$

$$\displaystyle\frac{p}{k\big{\{}F_{j}^{-1}(1-N^{-1})-F_{j}^{-1}(N^{-1})\big{\}}^{2}},\quad j=1,...,p.$$

$$|({\mathcal{X}}_{{\mathcal{D}}_{\mathrm{BS}}}^{\ast }{)}^{\mathrm{T}}{\mathcal{X}}_{{\mathcal{D}}_{\mathrm{BS}}}^{\ast }|\ge k(k-1{)}^p{\lambda }_{\min }^p({\mathbf{R}}_{{\mathcal{D}}_{\mathrm{BS}}})\prod _{j=1}^p\mathrm{var}(z_{j{\mathcal{D}}_{\mathrm{BS}}}^{\ast }).$$

$$(k-1)\mathrm{var}(z_{j{\mathcal{D}}_{\mathrm{BS}}}^{\ast })=\sum _{i=1}^k{\left(z_{ij}^{\ast }-{\bar{z}}_{j{\mathcal{D}}_{\mathrm{BS}}}^{\ast }\right)}^2=\sum _{b=1}^B\sum _{i=1}^{k_B}{\left(z_{bij}^{\ast }-{\bar{z}}_{j{\mathcal{D}}_{\mathrm{BS}}}^{\ast }\right)}^2$$

$$\ge$$

$$=$$

$$\ge$$

$$\mathrm{var}(z_{j{\mathcal{D}}_{\mathrm{BS}}}^{\ast })\ge$$

$$|({\mathcal{X}}_{{\mathcal{D}}_{\mathrm{BS}}}^{\ast }{)}^{\mathrm{T}}{\mathcal{X}}_{{\mathcal{D}}_{\mathrm{BS}}}^{\ast }|\ge$$

$$\times \prod _{j=1}^p\sum _{b=1}^B{\left(\frac{z_{b(n_B-r_B+1)j}-z_{b(r_B)j}}{z_{(N)j}-z_{(1)j}}\right)}^2.$$

$$\frac{|({\mathcal{X}}_{{\mathcal{D}}_{\mathrm{BS}}}^{\ast }{)}^{\mathrm{T}}{\mathcal{X}}_{{\mathcal{D}}_{\mathrm{BS}}}^{\ast }|}{M_N}\ge$$

$$(k-1)\mathrm{var}(z_{j{\mathcal{D}}_{\mathrm{BS}}}^{\ast })\ge$$

$$\ge$$

$$\mathrm{var}(z_{j{\mathcal{D}}_{\mathrm{BS}}}^{\ast })\ge$$

$$\frac{|({\mathcal{X}}_{{\mathcal{D}}_{\mathrm{BS}}}^{\ast }{)}^{\mathrm{T}}{\mathcal{X}}_{{\mathcal{D}}_{\mathrm{BS}}}^{\ast }|}{M_N}\ge$$

Full text

1 Introduction

In big data analysis, a common challenge is that available computing facilities are inadequate to meet the computational needs. To overcome this challenge, there are two fundamental approaches: one is the divide-and-conquer approach in which a large data set is divided into smaller partitions and results are calculated from each partition and then combined; the other approach is to use a subset of the full data as a surrogate, and use calculation results from the subdata to approximate the full data calculation results.

The divide-and-conquer approach naturally fits into parallel and distributed computing systems, and it is effective to deal with the challenge that the data volume is too large to be loaded into the memory although it may not always reduce the computing time with a single processor. Most distributed computational platforms are expensive to access and they are not suitable for daily routine usages such as exploratory data analysis and pilot model prototyping. The divide-and-conquer method has attracted many researchers in machine learning and statistics, leading to key advancements in Lin and Xie (2011); Jordan (2012); Chen and Xie (2014); Shang and Cheng (2017); Battey et al. (2018), among others. In the streaming setting where data blocks are accessible sequentially for only one time, the online updating method has been developed (Schifano et al., 2016; Xue et al., 2018). In order to produce an estimator that preserves the same convergence rate as the full data estimator, a key requirement for the divide-and-conquer method is that the number of partitions cannot be too large (Shang and Cheng, 2017).

The approach of using a subset of the full data is an effective method to reduce the computational burden, and it is often the only way to extract useful information from massive data sets when available computing power is limited. An advantage of this approach is that once a subset of data is taken, thorough analysis can often be performed on a regular computer such as a laptop. This is important in statistics because data analysis is not just computing. Note that the primary goal of selecting an informative subset from the full data agrees with the basic motivation of optimal design of experiments in which one wants to get as much information as possible with a fixed design budget. The major challenge of subdata selection from big data is how to identify informative data points computationally fast in order to maintain the computational advantage in the whole procedure of data analysis. For this purpose, statistical leverage scores are often used to define subsampling probabilities in linear regression (Drineas et al., 2012; Ma et al., 2015). In the context of logistic regression, Wang et al. (2018) derived the optimal subsampling probabilities that minimize the asymptotic variance of the subsampling estimator under the A- and L-optimality criteria, and Wang (2018) further developed a more efficient estimation approach based on optimal subsamples. Yao and Wang (2019) extended this technique to include the softmax regression.

The aforementioned work uses random subsampling to take subsets of full data, and Wang et al. (2019) showed that an estimator obtained from this approach for linear regression has a variance that is bounded from below by a term that is proportional to the inverse of the subset sample size, i.e., the estimator will not converge to the true parameter if the subset sample size is fixed no matter how fast the full data sample size goes to infinity. From the characterization of the optimal subdata under the D-optimality criterion, they further develop a computationally efficient IBOSS algorithm. The resulting estimator from this algorithm converges to the true parameter even if the subdata sample size is fixed as long as the full data sample size gets large and the supports of the covariate distributions are unbounded. The IBOSS algorithm selects data points by examining each column of the design matrix, so it is infeasible to apply if the data volume is too large to be loaded into the random-access memory (RAM). This paper combines the divide-and-conquer method and the IBOSS algorithm to solve this issue. The basic approach is to divide a large data set into smaller partitions so that each partition can be loaded into the RAM and the IBOSS algorithm is applied on the smaller partitions. In practice, big data are seldom stored in a single file; they are often stored in multiple data blocks. The proposed method fits this scenario naturally. We will investigate both finite sample properties and asymptotic properties of this procedure, and use numerical experiments to evaluate its empirical performance.

The rest of the paper is organized as the following. We present notations for the model setup and introduce the IBOSS method in Section 2. The main theoretical results on the divide-and-conquer IBOSS algorithm will be presented in Section 3. Section 4 evaluates the empirical performance of the divide-and-conquer IBOSS algorithm using numerical examples. Section 5 concludes the paper and all proofs are provided in the Appendix.

2 IBOSS framework and detailed algorithm

Assume that the full data $\mathcal{D}_{N}=\{(\mathbf{z}_{1},y_{1}),...,(\mathbf{z}_{N},y_{N})\}$ satisfy the following linear regression model

[TABLE]

where $y_{i}$ are the responses, $\mathbf{z}_{i}=(z_{i1}$, …, $z_{ip})^{\rm T}$ are the covariates, $\beta_{0}$ and $\beta_{1},...,\beta_{p}$ are the intercept and slope parameters, respectively, and $\varepsilon_{i}$ are model errors. For easy of presentation, we define some notations. Let $\bm{x}_{i}=(1,\mathbf{z}_{i})$ and $\mathcal{Z}=(\mathbf{z}_{1}^{\rm T},...,\mathbf{z}_{N}^{\rm T})^{\rm T}$. Denote $\bm{\beta}_{1}=(\beta_{1},\beta_{2},...,\beta_{p})^{\rm T}$ and $\bm{\beta}=(\beta_{0},\bm{\beta}_{1}^{\rm T})^{\rm T}$. Here we assume that the error terms $\varepsilon_{i}$ are uncorrelated and satisfy $\mathbb{E}(\varepsilon_{i})=0$ and $\mathbb{V}(\varepsilon_{i})=\sigma^{2}$. Without loss of generality, we assume that $\sigma^{2}=1$ when discussing the selection of subdata sets.

To estimate the unknown $\bm{\beta}$, the best linear unbiased estimator (BLUE) is the least squares estimator. With the full data $\mathcal{D}_{N}$, it has an expression of

[TABLE]

with variance-covariance matrix

[TABLE]

where $\bm{M}_{\mathrm{full}}$ is the observed Fisher information matrix of $\bm{\beta}$ for the full data if $\varepsilon_{i}$ are normally distributed.

When the sample size $N$ of the full data set is too large, full data analysis may take too long to afford. To extract some useful information from the data in time with limited computing resources, a subset of the full data can be selected and thoroughly analyzed. For this purpose, Wang et al. (2019) proposed the IBOSS method. We summarize the motivation and procedure in the following.

Use $\bm{\delta}=\{\delta_{1},\delta_{2},...,\delta_{N}\}$ to indicate a subset of data where $\delta_{i}=1$ if $(\mathbf{z}_{i},y_{i})$ is included in the subset and $\delta_{i}=0$ otherwise. Here $\bm{\delta}$ may depend on $\mathcal{Z}$ but it does not depend on $\bm{y}$, so the subdata selection rule is ancillary to the regression parameter, and the least squares estimator is still the BLUE based on the subset of the full data. With this notation, the information matrix for $\bm{\beta}$ based on a subset of the full data indexed by $\bm{\delta}$ is

[TABLE]

To extract the maximum amount of information from the full data with a fixed subset sample size $k$ so that $|\bm{\delta}|=\sum_{i=1}^{N}\delta_{i}=k$, Wang et al. (2019) proposed to select the subset of data that maximizes $\bm{M}_{\bm{\delta}}$ with respect to $\bm{\delta}$ under the the D-optimality criterion in optimal design of experiments, namely, to find

[TABLE]

In general, there is no analytical solution to $\bm{\delta}^{opt}_{\mathrm{D}}$, and numerical search for an exact solution is computationally NP-hard. Wang et al. (2019) derived an upper bound for $|\bm{M}_{\bm{\delta}}|$, and then proposed an algorithm to approximate it. To be specific, they showed that for any $\bm{\delta}$ with $|\bm{\delta}|=k$,

[TABLE]

where $z_{(N)j}$ and $z_{(1)j}$ are the sample maximum and sample minimum, respectively, for the $j$-th covariate. Although the upper bound $\zeta_{N}$ is typically unachievable, it indicates that the more informative data points are in the tail regions of the covariates. This motivated the IBOSS algorithm to approximate the upper bound, and here is a summary of the IBOSS procedure. Assume that $r=k/(2p)$ is an integer. For each covariate, using a partition-based selection algorithm (Martínez, 2004), select $r$ data points with the smallest $z_{ij}$ values and $r$ data points with the largest $z_{ij}$ values among the data points that are not yet included in the subsample. The data points obtained according to the $p$ covariates are then combined to form a subdata of size $k$ for statistical analysis.

The original IBOSS paper by Wang et al. (2019) did not present all practical details of implementing the IBOSS procedure. For example, it is mentioned that if some data points have been included in the subdata by some covariate, the IBOSS algorithm requires to exclude them from consideration when using other covariates to select data points, but the implementation details were not provided. This step may add significant CPU time if not implemented appropriately. Thus, for completeness, we present a detailed IBOSS algorithm taking into account necessary practical considerations in Algorithm 1.

3 Divide-and-conquer IBOSS algorithm

The IBOSS approach in Algorithm 1 needs to look at the data column by column. Although for most programming languages, such as R (R Core Team, 2018) and Julia (Bezanson et al., 2017), data matrices are stored by columns in the memory, data files are often stored on hard drive by rows. Thus, if the data volume exceeds the size of the available memory, it is difficult to implement the IBOSS algorithm. To deal with the problem of limited memory, a natural approach is the divide-and-conquer procedure, in which the full data set is divided into partitions by rows and then each partition can be loaded into the memory and processed individually. To combine this procedure with the IBOSS procedure, we present the divide-and-conquer IBOSS method in Algorithm 2.

Remark* 1**.*

Algorithm 2 takes a subdata $\mathcal{D}_{S}^{(b)}$ from each partition $\mathcal{D}_{N}^{(b)}$ and then combines them to have a final subdata $\mathcal{D}_{{\tiny{\mathrm{BS}}}}$. This is recommended if only one machine is used, and the combined subdata should be small enough so that thorough analysis can be performed on the machine. If a distributed computing system is used and data partitions are processed in parallel by multiple machines, we can calculate the least squares estimators and their estimated variance-covariance matrices on each subdata $\mathcal{D}_{S}^{(b)}$ in all machines and then aggregate these estimates using a linear combination with the inverses of the variance-covariance matrices as weights. The final estimators from these two approaches are identical (Lin and Xie, 2011; Schifano et al., 2016). The latter procedure takes advantage of parallel computing facilities, but it has to carry out additional calculations on each machine if further regression diagnostics are to be performed.

Remark* 2**.*

Using one machine to implement Algorithm 2, the step to divide the full data $\mathcal{D}_{N}$ into blocks $\mathcal{D}_{N}^{(1)}$, …, $\mathcal{D}_{N}^{(B)}$ can be done by using the UNIX command split. This command is also available for Windows through Cygwin.

3.1 Theoretical properties

In this section, we investigate the performance of Algorithm 2 theoretically. With out loss of generality, we assume that $r_{B}=r/B=k/(2pB)$ and $n_{B}=N/B$ are both integers, where $B$ is the number of partitions on the full data.

Let $z_{(i)j}$ ($i=1,...,N;j=1,...,p$) be the $i$-th order statistic on the $j$-th covariate in the full data and $z_{b(i)j}$ ($b=1,...,B;i=1,...,n_{B};j=1,...,p$) be the $i$-th order statistics on the $j$-th covariate in the $b$-th block. Denote the covariate matrix for the subdata $\mathcal{D}_{{\tiny{\mathrm{BS}}}}$ by $\mathcal{Z}^{*}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$ and denote the sample correlation matrix of $\mathcal{Z}^{*}_{\mathcal{D}_{S}}$ by $\mathbf{R}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$. Let $\mathcal{X}^{*}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}=(\mathbf{1},\mathcal{Z}^{*}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}})$.

For the divide-and-conquer IBOSS method, we first present a result showing the quality of using Algorithm 2 to approximate $\zeta_{N}$, a typically unachievable upper bound of $|\bm{M}_{\bm{\delta}}|$ defined in (2).

Theorem 1**.**

The subdata $\mathcal{D}_{{\tiny{\mathrm{BS}}}}$ selected using Algorithm 2 satisfies that

[TABLE]

where

[TABLE]

and $\lambda_{\min}(\mathbf{R}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}})$ is the smallest eigenvalue of $\mathbf{R}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$.

Remark* 3**.*

The lower bound of the approximation ratio in (3) has two terms: $C_{R}$ in (4) and $C_{E}$ in (5). Heuristically, $C_{R}$ is a bound corresponding to the scenario that the full data is divided randomly so that the covariate distributions for different partitions are the same. In this scenario $C_{R}$ is a sharper bound compared with $C_{E}$, i.e., $C_{R}>C_{E}$. On the other hand, $C_{E}$ is a bound corresponding to some extreme cases, e.g., ranges of covariate values for all blocks do not overlap. In this case, $C_{E}$ may be larger than $C_{R}$.

Theorem 1 is aligned with Theorem 3 of Wang et al. (2019) for the original IBOSS algorithm, which shows that for the subdata $\mathcal{D}_{S}$ obtained from Algorithm 1, the following inequality holds.

[TABLE]

where $\mathcal{X}^{*}_{\mathcal{D}_{S}}=(\mathbf{1},\mathcal{Z}^{*}_{\mathcal{D}_{S}})$, $\mathcal{Z}^{*}_{\mathcal{D}_{S}}$ is the covariate matrix for the subdata $\mathcal{D}_{S}$, and $\lambda_{\min}^{p}(\mathbf{R}_{\mathcal{D}_{S}})$ is the minimum eigenvalue of $\mathbf{R}_{\mathcal{D}_{S}}$, the correlation matrix of $\mathcal{Z}^{*}_{\mathcal{D}_{S}}$. Comparing (3) and (6), we see that there may be some information loss due to the divide-and-conquer procedure because $(z_{b(n_{B}-r_{B}+1)j}-z_{b(r_{B})j})^{2}$ and $(z_{(N-r_{B}+1)j}-z_{(r_{B})j})^{2}/B$ are typically smaller than $(z_{(N-r+1)j}-z_{(r)j})^{2}$. However, if $B$ is not too large, $z_{b(n_{B}-r_{B}+1)j}-z_{b(r_{B})j}$ and $z_{(N-r+1)j}-z_{(r)j}$ are on the same order under reasonable assumptions, so the lower bounds of the approximation ratios are at the same order. Theorem 1 also indicates that $|(\mathcal{X}^{*}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}})^{\rm T}\mathcal{X}^{*}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}|$ may achieve the same order of the unachievable upper bound $\zeta_{N}$ under some reasonable assumptions if $B$ and $p$ do not go to infinity.

Now we investigate the properties of the resulting estimator form Algorithm 2. Let $\hat{\bm{\beta}}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}=(\hat{\beta}_{0}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}},\hat{\beta}_{1}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}},...,\hat{\beta}_{p}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}})^{\rm T}$ be the least squares estimator calculated from the divide-and-conquer IBOSS subdata, namely,

[TABLE]

where $\bm{x}_{bi}^{*}=(1,\mathbf{z}_{bi}^{*})^{\rm T}$ and $(\mathbf{z}_{bi}^{*},y_{bi}^{*})$ is the $i$-th observation in the subdata $\mathcal{D}_{S}^{(b)}$ from the $b$-th block. To investigate the theoretical properties of $\hat{\bm{\beta}}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$, we focus on the variances because $\hat{\bm{\beta}}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$ is unbiased for $\bm{\beta}$.

The following theorem provides finite sample bounds on variances of $\hat{\beta}_{j}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$, $j=0,1,...,p$. No quantity is required to go to infinity here.

Theorem 2**.**

If the sample correlation matrix for covariates $\mathcal{Z}^{*}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$ in the subdata $\mathcal{D}_{{\tiny{\mathrm{BS}}}}$ is column full rank, then the following results hold for the estimator, $\hat{\bm{\beta}}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$, obtained from Algorithm 2:

[TABLE]

with probability one for $j=1,...,p$, where

[TABLE]

and $\lambda_{\min}(\bm{R}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}})$ is the minimum eigenvalue of the sample correlation matrix of $\mathcal{Z}^{*}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$.

Theorem 2 shows that for the intercept, the variance of the estimator $\hat{\beta}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}_{0}$ is bounded from below by $\sigma^{2}/k$, which is the same as the bound for the original IBOSS algorithm. It indicates that this variance cannot go to zero for a fixed $k$. However, for slope parameters, the variances of the estimators $\hat{\beta}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}_{j}$ may converge to zero under some assumptions even $k$ is fixed.

Theorem 2 is aligned with Theorem 4 of Wang et al. (2019), which shows that the variances of the slope estimators from the original IBOSS procedure are bounded from above by

[TABLE]

Comparing the upper bounds $V_{aj}$ and $V_{ej}$ with the upper bounds $V_{oj}$ for the original IBOSS procedure, we also see that the divide-and-conquer IBOSS algorithm may subject to some information loss. However, the information loss can be ignored asymptotically if the full data is partitioned randomly and the number of partitions does not go to infinity too fast. We derive the following theorem showing that the orders of $\mathbb{V}(\hat{\beta}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}_{j}|\mathcal{Z})$ can be the same as those of the variances for the original IBOSS estimators.

Theorem 3**.**

Assume that covariate distributions are in the domain of attraction of the generalized extreme value distribution, and $\underset{N\rightarrow\infty}{\lim\inf}\lambda_{\min}(\bm{R}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}})>0$. For large enough $N$, when the full data is divided randomly in Algorithm 2, the following results hold for the estimator, $\hat{\beta}_{j}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$, $j=1,...,p$.

[TABLE]

if one of the following conditions holds: 1) $r$ and $B$ are fixed; 2) the support of $F_{j}$ is bounded and $r/N\rightarrow 0$, where $F_{j}$ is the marginal distribution function of the $j$-th component of $\mathbf{z}$; 3) the upper endpoint for the support of $F_{j}$ is $\infty$ and the lower endpoint for the support of $F_{j}$ is finite, and

[TABLE]

for all $\epsilon>0$; 4) the upper endpoint for the support of $F_{j}$ is finite and the lower endpoint for the support of $F_{j}$ is $-\infty$, and

[TABLE]

for all $\epsilon>0$; and 5) the upper endpoint and the lower endpoint for the support of $F_{j}$ are $\infty$ and $-\infty$, respectively, and both (9) and (10) hold.

The convergence rates of $\mathbb{V}(\hat{\beta}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}_{j}|\mathcal{Z})$ described in Theorem 3 are the same as those given in Theorem 5 of Wang et al. (2019) for the original IBOSS estimators of slope parameters.

For case 1), if $p$ is fixed, then the total subdata size $k$ is also fixed. As long as the supports of covariate distributions are not bounded, $z_{(N-r+1)j}-z_{(r)j}\rightarrow\infty$ in probability as $n\rightarrow\infty$, and thus the variances $\mathbb{V}(\hat{\beta}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}_{j}|\mathcal{Z})$ converge to zero in probability. This is not the case for random subsampling-based methods (Ma et al., 2015; Wang et al., 2018). In this case, the asymptotic expression for the variance-covariance matrix of $\hat{\bm{\beta}}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}$ when $\mathbf{z}_{i}$ follow normal and lognormal distributions are identical to those in Theorem 6 of Wang et al. (2019). For case 2), $z_{(N-r+1)j}-z_{(r)j}$ goes to some finite constant, so it is necessary that $k\rightarrow\infty$ in order for $\mathbb{V}(\hat{\beta}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}_{j}|\mathcal{Z})$ to converge to zero. For case 3), the condition in (9) impose constrains on $r$, $B$, $N$, and the tail behavior of the covariate distributions. For example, if

[TABLE]

where $\gamma>0$ and $h(z)$ is a slowly varying function, then the first constrain in (9) is to require that $\log r/\log N\rightarrow 0$ (Hall, 1979). The distribution in (11) includes many commonly seen distributions, such as exponential distribution, gamma distribution, Gumbel distribution, Laplace distribution, and normal distribution. For these distributions, if $\log B/\log N\rightarrow 0$, then the second constrain in (9) also holds. For case 4), the condition in (10) is essentially the same as that in (9) if one takes transformation $\mathbf{z}=-\mathbf{z}$. For case 5, it is a combination of cases 3 and 4. Under the given conditions, $F_{j}^{-1}(1-N^{-1})/z_{(N)j}\rightarrow 1$ and $F_{j}^{-1}(N^{-1})/z_{(1)j}\rightarrow 1$, so the rates for $\mathbb{V}(\hat{\beta}^{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}_{j}|\mathcal{Z})$ to converge to zero are

[TABLE]

4 Numerical experiments

In this section, we use numerical experiments to evaluate the empirical performance of the divide-and-conquer IBOSS algorithm.

We generate data from model (1) with $\beta_{0}=1$ and $\bm{\beta}_{1}$ being a 50 dimensional vector of ones. We assume that $\varepsilon_{i}$ are independent following $N(0,1)$. The following five cases are considered to generate the covariate matrix $\mathcal{Z}$.

1. C1se1:

Normal. We generate $\mathcal{Z}$ from a multivariate normal distribution, $N(\bm{0},\bm{\Sigma})$, where the $(i,j)$-th element of $\bm{\Sigma}$ is 0.5 if $i\neq j$ and 1 otherwise. 2. C2se2:

LogNormal. We generate $\mathcal{Z}$ from a multivariate lognormal distribution, namely, generate $\mathcal{V}$ from $N(\bm{0},\bm{\Sigma})$ as defined in Case 1 and then set $\mathcal{Z}=e^{\mathcal{V}}$, where the exponentiation is element-wise. 3. C3se3:

$\bm{T_{2}}$. We generate $\mathcal{Z}$ from a multivariate $t$ distribution with two degrees of freedom $t_{2}(\bm{0},\bm{\Sigma})$ with $\bm{\Sigma}$ being the same as in Case 1. 4. C4se4:

Mix with order. We generate $\mathcal{Z}$ from five different distributions. Specifically, generate $\mathcal{Z}_{1}$ from $N(\bm{0},\bm{\Sigma})$; generate $\mathcal{Z}_{2}$ from $t_{2}(\bm{0},\bm{\Sigma})$; generate $\mathcal{Z}_{3}$ from $t_{3}(\bm{0},\bm{\Sigma})$; generate $\mathcal{Z}_{4}$ with its elements independently following the same uniform distribution between zero and two; and generate $\mathcal{Z}_{5}$ from the multivariate lognormal distribution defined in Case 2. Here, $\mathcal{Z}_{i}$ ($i=1,...,5$) all contain $n/5$ rows, and $\bm{\Sigma}$ is the same as defined in Case 1. Set $\mathcal{Z}=(\mathcal{Z}_{1}^{\rm T},\mathcal{Z}_{2}^{\rm T},\mathcal{Z}_{3}^{\rm T},\mathcal{Z}_{4}^{\rm T},\mathcal{Z}_{5}^{\rm T})^{\rm T}$. 5. C5se5:

Mix random order. We generate $\mathcal{Z}$ using the same procedure as in Case 4, and then we randomize the row orderings.

Here, Cases 1-4 are the same as those in Wang et al. (2019). Cases 4 and 5 would produce identical results for the full data analysis approach but produce different results for the divide-and-conquer IBOSS algorithm as different row orderings will affect the data partitions. The purpose of considering both Cases 4 and 5 is to investigate the effect of different data partitions on the the divide-and-conquer IBOSS algorithm. We divide the full data according to the sequence of row numbers in the numerical experiments.

We repeat the simulation for $T=1000$ times to calculate the empirical mean squared error (MSE) as MSE=$T^{-1}\sum_{t=1}^{T}\|\bm{\beta}_{1}^{(t)}-\bm{\beta}_{1}\|^{2}$ for the slope parameter, where $\bm{\beta}_{1}^{(t)}$ is the estimate at the $t$-th repetition.

Figure 1 plots log of empirical MSEs against full data sample size $N$ of $N=5\times 10^{3},10^{4},10^{5}$ and $10^{6}$, with total subdata size $k=1000$. Here the natural logarithm is taken to have better presentations. For the number of partitions $B$, we considered four values: $B=1,2,5$, and $10$, and the corresponding values of $r_{B}$ are $r_{B}=10,5,2$, and $1$, respectively. Note that when $B=1$, the method is the original IBOSS approach. For comparison, the uniform Poisson subsampling method is also implemented. From Figure 1, the empirical MSE decreases as $N$ gets large for the divide-and-conquer IBOSS approach, while it stays constant for the random sampling approach. Different values of $B$ have some effect on the performance of divide-and-conquer IBOSS approach, the effect is not very significant, especially for Cases 2, 3, and 5. Comparing results in Cases 4 and 5, we see that different row orderings affect the performance of the divide-and-conquer IBOSS algorithm, and the random row ordering (random partitions) is preferable.

We also consider the scenario that $r_{B}$ is fixed while $B$ changes. This is to mimic the scenario that one may have $B$ machines to use while the computational capacity for each machine is limited. This is often the case for parallel distributed computing systems. We set $r_{B}=5$ and choose $B$ to be $B=1$, 2, 5, 10, and 20, which gives the values of $k$ as $k=500$, 1000, 2500, 5000, and 10000, respectively. We fixed $N=10^{6}$. Figure 2 presents the results. It is seen that both the divide-and-conquer method and the Poisson subsampling method improve as $B$ increases because the total subdat size $k$ increases. The divide-and-conquer method dominates the Poisson subsampling method and the difference in their performance depends on the covariate distribution.

To further explore the effect of larger number of partitions $B$ with a fixed $k=2000$, we reduce the value of $p$ to be $p=2$ so that the maximum value of $B$ is $B=500$ corresponding to $r_{B}=1$. We consider $B=1$, 2, 5, 10, 20, 50, 100, 125, 250, and 500, which gives values of $r_{B}$ as $r_{B}=500$, 250, 100, 50, 25, 10, 5, 4, 2, and 1, respectively. Figure 3 gives the results. We see that the effect of $B$ on the performance of the divide-and-conquer IBOSS algorithm is small compared with its significant improvement relative to the Poisson subsampling method. Additionally, there is no clear increasing pattern for the empirical MSE as $B$ gets large. This is different from the observations for regular cases in the existing literature where the estimation efficiency of the divide-and-conquer estimators often decreases as the number of partitions increases (e.g. Lin and Xie, 2011; Schifano et al., 2016; Battey et al., 2018).

We investigate the computing times of the divide-and-conquer IBOSS algorithm for the case that data blocks are stored on hard drive. We set $N=10^{7}$, $p=100$, $k=10^{4}$, and choose values of $B$ to be $B=2,5,10,25$, and $50$, giving values of $r_{B}=25,10,5,2$, and $1$, respectively. For comparisons, we implement the full data divide-and-conquer approach in which all observations in each block of data are used to calculate least squares estimators and these estimators are then aggregated to form a full data estimator. We also implement the Poisson subsampling method on each data block to obtain subdata sets and then combine these subdata sets to calculate an estimate. Table 1 presents the results on computing times, where the time to load each data block into the memory is also counted. The calculation was performed on a computer with Intel Core i7 processors, a 64 GB RAM, and SSD hard drives, running Ubuntu 18.04 Linux system. Only one CPU was used for fair comparisons. It is seen that the divide-and-conquer IBOSS algorithm is faster than the full data approach, but is not as fast as the the uniform subsampling approach. For the divide-and-conquer IBOSS approach and the full data approach, as the number of blocks increases, the computing times are reduced. However, this is not the case for the uniform subsampling approach, for which the computing time decreases first but then increases as $B$ increases. This is because as $B$ gets large, the data loading time becomes the dominating term for the time required by the divide-and-conquer uniform Poisson subsampling method.

5 Concluding remarks

In this paper, we have developed a divide-and-conquer IBOSS method in the context of linear regression to solve the issue that the data volume is too large to apply the D-optimality motivated IBOSS method in memory. We have derived theoretical guarantees for the resultant estimator and drawn comparisons with the original IBOSS approaches. We have also carried out numerical experiments to evaluate and demonstrate the proposed method.

We found that the estimation efficiency of the divide-and-conquer IBOSS algorithm may not be a monotonic function of the number of partitions as in other divide-and-conquer methods, which indicates the possibility of a nontrivial optimal number of partitions. This is an interesting and challenging question warrants for further investigations.

6 Proofs

6.1 Proof of Theorem 1

Proof.

For $l\neq j$, let $z_{j}^{(i)l}$ be the concomitant of $z_{(i)l}$ for $z_{j}$, i.e., if $z_{(i)l}=z_{sl}$ then $z_{j}^{(i)l}=z_{sj}$, $i=1,...,N$. Let $\bar{z}_{j\mathcal{D}_{{\tiny{\mathrm{BS}}}}}^{*}$ and $\mathrm{var}(z_{j\mathcal{D}_{{\tiny{\mathrm{BS}}}}}^{*})$ be the sample mean and sample variance for covariate $z_{j}$ of subdata $\mathcal{D}_{{\tiny{\mathrm{BS}}}}$. For the proof of Theorem 3 in Wang et al. (2019), we know that

[TABLE]

Let $k_{B}=k/B$ be the number of data points to take from each partition, and $z_{bij}^{*}$ be the $i$-th observation on the $j$-th covariate in the subdata $\mathcal{D}_{S}^{(b)}$ from the $b$-th partition. The sample variance for each $j$ satisfies,

[TABLE]

where $\bar{z}_{bj}^{**}=\left(\sum_{i=1}^{r_{B}}+\sum_{i=n_{B}-r_{B}+1}^{n_{B}}\right)z_{b(i)j}/(2r_{B})$, $\bar{z}_{bj}^{*l}=\sum_{i=1}^{r_{B}}z_{b(i)j}/{r_{B}}$,

and $\bar{z}_{bj}^{*u}=\sum_{i=n_{B}-r_{B}+1}^{n_{B}}z_{b(i)j}/(r_{B})$. Thus,

[TABLE]

which, combined with (12), shows that

[TABLE]

This shows that

[TABLE]

Each numerator in the summation of the bound in (14) relay on the covariate range of each data partition. If the full data is not divided randomly, this may not produce a sharp bound and using the full data covariate ranges may produce a better bound. We use this idea to derive the bound $C_{E}$ in the following. From Algorithm 2, for each $j=1,...,p$, the $r_{B}$ data points with the smallest value of $z_{j}$ and the $r_{B}$ data points with the largest value of $z_{j}$ must be included in $\mathcal{D}_{{\tiny{\mathrm{BS}}}}$. Thus, for each sample variance,

[TABLE]

In this case, $\bar{z}_{j}^{**}=(\sum_{i=1}^{r_{B}}+\sum_{i=N-r_{B}+1}^{N})z_{(i)j}/(2r_{B})$, $\bar{z}_{j}^{*l}=\sum_{i=1}^{r}z_{(i)j}/(r_{B})$,

and $\bar{z}_{j}^{*u}=\sum_{i=n-r+1}^{n}z_{(i)j}/(r_{B})$. Thus,

[TABLE]

This, combined with (12), shows that

[TABLE]

The proof finishes if we combine (14) and (16). ∎

6.2 Proof of Theorem 2

Proof.

The proof for (7) is similar to the proof of inequality (19) in Wang et al. (2019). For (8), from the proof of Theorem 4 in Wang et al. (2019), we know that

[TABLE]

and

[TABLE]

From the (17), (18), and the fact that $(\bm{R}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}^{-1})_{jj}\geq 1$, we have

[TABLE]

From (13), (15) and the fact that $(\bm{R}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}}^{-1})_{jj}\leq\lambda_{\min}^{-1}(\bm{R}_{\mathcal{D}_{{\tiny{\mathrm{BS}}}}})$, we have

[TABLE]

The proof finishes by combining (19), (20), and (21). ∎

6.3 Proof of Theorem 3

Proof.

For the first case that $B$ and $r$ are fixed, from (8) with $V_{ej}$, we only need to verify that

[TABLE]

which is true according to Theorems 2.8.1 and 2.8.2 in Galambos (1987).

For the second case, $z_{b(n_{B}-r_{B}+1)j}-z_{b(r_{B})j}$ and $z_{(N)j}-z_{(1)j}$ converge to the same fixed constant in probability, and $z_{b(n_{B}-r_{B}+1)j}-z_{b(r_{B})j}$ are bounded by this constant for all $b$. Thus,

[TABLE]

converges to the same finite constant. Thus, (22) can be easily verified.

For the third case, let $g_{N,j}=F_{j}^{-1}(1-1/N)$ and $g_{n_{B},j}=F_{j}^{-1}(1-1/n_{B})$. When (9) holds, from the proof of Theorem 5 in Wang et al. (2019), we have $z_{(N)j}/g_{N,j}=1+o_{P}(1)$ and $z_{(n_{B}-r_{B}+1)j}/g_{n_{B},j}=1+o_{P}(1)$. Thus, noting that $z_{b(r)j}$ and $z_{(1)j}$ are bounded in probability when the lower endpoint for the support of $F_{j}$ is finite, we have

[TABLE]

Note that $\left|\frac{z_{b(n_{B}-r_{B}+1)j}-z_{b(r)j}}{z_{(N)j}-z_{(1)j}}\right|$ are bounded by the same constant for all $b$, so

[TABLE]

Combining this and (8) with $V_{aj}$, the result follows.

For the forth case, the result follows from reversing the signs of the covariates in the proof for the third case.

The proof for the fifth case is obtained by combining the proof for the third case and the proof for the fourth case. ∎