# Speaker embeddings by modeling channel-wise correlations

###### Abstract

Speaker embeddings extracted with deep 2D convolutional neural networks are
typically modeled as projections of first and second order statistics of
channel-frequency pairs onto a linear layer, using either average or attentive
pooling along the time axis. In this paper we examine an alternative pooling
method, where pairwise correlations between channels for given frequencies are
used as statistics. The method is inspired by style-transfer methods in
computer vision, where the style of an image, modeled by the matrix of
channel-wise correlations, is transferred to another image, in order to
produce a new image having the style of the first and the content of the
second. By drawing analogies between image style and speaker characteristics,
and between image content and phonetic sequence, we explore the use of such
channel-wise correlations features to train a ResNet architecture in an end-
to-end fashion. Our experiments on VoxCeleb demonstrate the effectiveness of
the proposed pooling method in speaker recognition111Submitted to
Interspeech-2021.

Index Terms: speaker recognition, style-transfer, deep learning

##  1 Introduction

Extracting low-dimensional representations from speech utterances that
characterize the speaker is a challenging tasks with numerous application.
Apart from speaker recognition, speaker representations are used in automatic
speech recognition (e.g. for speaker adaptation), in speaker diarization and
separation, and in speech synthesis, for mimicking a target speaker using
text-to-speech or voice conversion.

Until recently, i-vectors were the most commonly used means for representing
speaker characteristics [1, 2]. With the advent of deep learning, neural
architectures such as time-delay neural networks (TDNNs [3]), Long Short-Term
Memory networks (LSTMs), 2D convolutional neural networks (such as ResNets and
its extensions [4]) and Transformers have emerged, surpassing the performance
of i-vectors in most applications and speaker recognition settings.

In order to obtain fixed-size representation vectors from utterances of
variable duration and arbitrary word sequences (i.e. text-independent), a
pooling method is required. The most common pooling method is a simple
aggregation over the time axis, where the mean vector is concatenated with
standard deviation (known as statistical pooling layer). The speaker embedding
(a vector of few hundred dimensions) is obtained by projecting this pooled
representation (a vector of few thousand dimensions) onto a linear layer or a
shallow feed-forward neural network. Other pooling approaches have been
proposed, such as single and multi-head attention, as well as the use of
higher-order statistics [5, 6, 7].

In this paper we present an alternative way of modeling speaker information,
that is largely inspired by style-transfer [8]. This seminal algorithm in
computer vision is capable of producing new images that combine the content of
an image with the style of other images or artworks. The style of an image is
modeled by the correlations between the channels of the convolutional neural
network. By considering styles as broadly analogous to speakers, we experiment
with modeling speakers by channel-wise correlation. This is in contrast to
averaging along the time dimension of each channel frequency pair of the
output tensor that is used in statistical pooling. The other sources of
variability that can also be considered as style of the recording, such as
noise and session effects222Session effects are related to the recording
device and they are also known as channel effects in speaker recognition
literature. are suppressed by training the model end-to-end using supervised
learning with speakers as targets. The resulting embeddings are highly
speaker-discriminative, yielding state-of-the-art results on VoxCeleb using
plain cosine similarity.

The rest of the paper is organized as follows. In Sect. 2 we briefly review
some of the related work, such as pooling methods and style-transfer. In Sect.
3 we introduce the proposed method, and provide all the implementation and
architectural details. In Sect. 4, we present our experiments on a standard
speaker recognition benchmark (VoxCeleb). Finally, in Sect. 5 we provide
conclusions and future work directions.

##  2 Related work

###  2.1 Style-transfer in computer vision and voice conversion

A few years ago, the seminal work of L.A. Gatys et al. [8] in neural style-
transfer gained the attention of the computer vision community and the
industry, and received wide media coverage. Since then, several works extended
their method, some of which can be found in [9]. A theoretical interpretation
of the way style and image texture is encoded in the channel-wise correlations
(which can be considered as a Gram-Matrix) is given in [10]. The authors treat
style-transfer as a domain adaptation problem, and prove the equivalence
between trying to match the Gram-Matrices of feature maps and minimizing the
Maximum Mean Discrepancy (MMD) with the second-order polynomial kernel.

One of the few attempts to apply neural-style transfer in speech synthesis is
described in [11]. The authors employ a pretrained CNN and apply style-
transfer to train a speech generation and voice conversion system. Moreover,
they demonstrate the capacity of the intermediate tensors to encode speaker
characteristics using the VCTK dataset. However, the scope of their work is
different to ours, the network is not trained end-to-end while the
representations are high-dimensional. Furthermore, their method is not tested
on a speaker recognition benchmark containing noise and high intrinsic and
extrinsic within-speaker session variability (such as VoxCeleb or NIST-SRE).
Hence, the ability of their utterance-level representation to suppress session
variability has not been examined.

###  2.2 Pooling layers in speaker embeddings

The vast majority of speaker embedding networks use the statistical pooling
layer introduced by x-vectors, i.e. mean and standard deviation (std) for each
frequency and channel combination [12]. The main family of alternative pooling
methods enhance these statistics with an attention mechanism. Examples are
single and multi-head attention, as well as Net and Ghost-VLAD [6, 13, 14].
Both attention and VLAD methods employ a set of heads (i.e. trainable network
parameters which are broadly analogous to the means of a mixture model) to
cluster the output features into broad phonetic categories in an unsupervised
and end-to-end fashion. The resulting statistical pooling is defined as the
concatenation of the per-head statistics. Their main difference lies in the
normalization; the attention normalizes only over the time axis while the VLAD
over the attentive heads. Second order statistics (i.e. std features) may also
be added as in the case of average pooling. More recent approaches exist, such
as vector-based attentive pooling [15], phonetically-aware attention [16], and
cross-attention [17].

A second alternative to standard statistical pooling is the use of higher-
order statistics (skewness and kurtosis) [7]. Skewness is a measure of the
asymmetry of a distribution with respect to its mode, while the kurtosis
measures its “tailedness”. The authors did not obtain any consistent
improvement compared to using merely first and second statistics, as in
standard statistical pooling. Recently (and concurrently to our work) a set of
alternative pooling methods were examined in [18]. The authors show that the
use of second order statistics is more effective compared to first order,
while they also perform experiments with using the full covariance matrix,
without showing improvements over the standard statistical pooling. Their work
is close to ours, although we differ in several ways, such as in the suggested
pooling method and in the several implementation details (e.g. instead of
covariances we propose frequency-dependent channel-wise correlation matrices).

##  3 The proposed method

###  3.1 Notation

We denote by ${\mathbf{X}}=\\{{\mathbf{x}}_{t}\\}_{t=1}^{T_{i}}$ the sequence
of acoustic features of an utterance with $T_{i}$ frames, which in our method
are $F_{i}=80$ dimensional filterbank features. The architecture is a
ResNet-34 and the output of the last ResNet block is a 3D tensor (we omit the
batch axis) denoted by ${\mathbf{Y}}\in{\mathbf{R}}^{T\times F\times C}$,
where $T$ denotes the size of the temporal axis, $F$ the size of the frequency
axis, and $C$ the number of channels. Typical values for these quantities are
$(T,F,C)=(50,10,256)$; note though that $T$ is variable and can be derived by
$T_{i}$ divided by the cumulative temporal stride (equal to $2^{3}=8$ in our
architecture). Similarly, $F=10$ since the cumulative frequency stride is
again equal to 8, while $F_{i}=80$.

###  3.2 Statistical Pooling

The statistical pooling layer in speaker embeddings networks with 2D CNN
architectures is a concatenation of the mean and std of each of the $F\times
C$ frequency-channel pairs $[{\mathbf{P}}^{m};{\mathbf{P}}^{s}]$, where

${\mathbf{P}}^{m}_{f,c}=\frac{1}{T}\sum_{t}{\mathbf{Y}}_{t,f,c}$ (1)

and

${\mathbf{P}}^{s}_{f,c}=\sqrt{\frac{1}{T}\sum_{t}({\mathbf{Y}}_{t,f,c}-{\mathbf{P}}^{m}_{f,c})^{2}}$
(2)

The speaker embeddings are extracted by projecting the (vectorized) pooling
layer onto a lower dimensional space (e.g. $d_{e}=256$) using a single linear
layer.

###  3.3 Style modeling in images

The style of an image in computer vision for a given convolutional layer is
modeled by the symmetric positive-semidefinite matrix ${\mathbf{S}}$ with
elements defined as

${\mathbf{S}}_{c,c^{\prime}}=\frac{1}{TF}\sum_{t,f}{\mathbf{Y}}_{t,f,c}{\mathbf{Y}}_{t,f,c^{\prime}}$
(3)

Note though that images differ from spectral representations of acoustic
signals. Image characteristics are considered invariant along both spatial
axes, while spectral representations are invariant only along the time axis.
In fact, this property explains the use of frequency-channel pairs in
statistical pooling instead of channels alone, as shown in (1)-(2). Due to
this property of spectral representations, pooling along both time and
frequency axes as in (3) is expected to yield less discriminative speaker
representations.

###  3.4 Pooling via frequency-dependent channel-wise correlations

A natural modification of (3) that addresses the above issue is to pool only
along the time axis, i.e.

${\mathbf{S}}_{f,c,c^{\prime}}=\frac{1}{T}\sum_{t}{\mathbf{Y}}_{t,f,c}{\mathbf{Y}}_{t,f,c^{\prime}}$
(4)

The 3D tensor ${\mathbf{S}}$ can be viewed as a list of $F$ symmetric and
positive semidefinite matrices ${\mathbf{S}}_{f}={\mathbf{S}}_{f,:,:}$ and has
therefore $FC(C+1)/2$ free variables. As its size can be too large we propose
to reduce it in a frequency-dependent manner as shown later.

A further proposed modification is the use of mean and variance normalization
along the time axis of ${\mathbf{Y}}$ for each frequency-channel pair, i.e.
${\mathbf{Z}}_{:,f,c}\xleftarrow{norm}\ {\mathbf{Y}}_{:,f,c}$. When variance
normalization is applied, the variables in the diagonal of ${\mathbf{S}}_{f}$
(i.e. the variance of each channel-frequency combination) become equal to 1,
and therefore they can be removed from the pooling layer resulting in
$C(C-1)/2$ free variables for each frequency $f$ and $FC(C-1)/2$ overall.
Moreover, ${\mathbf{S}}_{f,c,c^{\prime}}\in[-1,1]$ for all $f,c,c^{\prime}$.
The pooling layer is followed by a linear layer which reduces the size of the
representation to $d_{e}=256$ (the speaker embedding space), as in the
standard statistical pooling.

Note that when mean and variance normalization are applied, the pooled
representation encodes complementary information to that of the standard
pooling; the means of each channel-frequency are being subtracted, while their
variance is also set equal to one.

###  3.5 Additional extensions and modifications

In the pooling method described above, we add the following operations in
other to decrease the size of the pooling layer and regularize the network.
Recall that regularization is needed in order to extract embeddings that
generalize well to speakers unseen during training.

####  3.5.1 Channel dropout

Dropout is a standard choice for regularizing a network. We employ it in a
channel-wise fashion, that is, we sample a 1D mask of size $C$ (with dropout
probability $p_{d}=0.25$) and we broadcast it along the time and frequency
axes. Channel-wise dropout is chosen because of its Bayesian interpretation as
model averaging between models with variable number of convolutional filters
[19].

####  3.5.2 Frequency ranges

As neighboring frequency bins in ${\mathbf{Y}}$ have similar statistics, we
consider merging them into $F_{r}$ frequency ranges. After experimentation, we
concluded that merging $f_{r}=2$ consecutive frequency bins in a non-
overlapping manner yields more stable performance. Merging is implemented by
reshaping ${\mathbf{Y}}$, so that $F_{r}\leftarrow Ff_{r}^{-1}$ and
$T_{r}\leftarrow Tf_{r}$. For example, when experimenting with $f_{r}=2$, we
obtain $F_{r}=5$ and $T_{r}=100$ given that $F=10$ and $T=50$. Note that by
setting $f_{r}=F$ we recover the original style-transfer used in computer
vision, where pooling is performed along both spatial axes.

####  3.5.3 Frequency-dependent channel reduction

The number of channels in ${\mathbf{Y}}$ is $C=256$ which would result in
$C(C-1)/2\approx 33K$ variables per frequency bin. So we transform
${\mathbf{Y}}$ as follows

${\mathbf{Y}}_{t,f,c^{\prime}}\leftarrow\sum_{c}{\mathbf{L}}_{f,c,c^{\prime}}{\mathbf{Y}}_{t,f,c}.$
(5)

where ${\mathbf{L}}$ is a learnable tensor of shape $[F_{r},C,C^{\prime}]$,
acting on the channel axis and reducing its size from $C$ to $C^{\prime}$ (we
choose $C^{\prime}=64$). The rationale for using a 3D kernel instead of a 2D
is that different linear combinations of channels should be better suited for
obtaining speaker discriminative information from each frequency range.

###  3.6 Summary and order of operations

We summarise here the set of operations and their order we use in order to
extract speaker embeddings. Recall that the output tensor of the last ResNet
block is ${\mathbf{Y}}$ with shape $[T,F,C]$.

  1. 1.

Dropout We apply channel-wise dropout with $p_{d}$ by sampling a binary mask
of shape $C$ and broadcasting it along the other two axes.

  2. 2.

Frequency ranges We reshape the tensor ${\mathbf{Y}}$ with new shape
$[T_{r},F_{r},C]$, by merging $f_{r}$ consecutive frequency bins into $F_{r}$
frequency ranges.

  3. 3.

Channel Reduction We reduce the channel size by applying ${\mathbf{L}}$ as
shown in (5), i.e. ${\mathbf{Y}}\xleftarrow{{\mathbf{L}}}{\mathbf{Y}}$, where
the shape of ${\mathbf{Y}}$ becomes $[T_{r},F_{r},C^{\prime}]$ .

  4. 4.

Normalization We apply mean and variance normalization along the time axis,
i.e. ${\mathbf{Z}}\xleftarrow{norm}{\mathbf{Y}}$.

  5. 5.

Pooling We apply pooling as shown in (4), but using ${\mathbf{Z}}$ instead of
${\mathbf{Y}}$, by which we obtain ${\mathbf{S}}$, a 3D tensor with shape
$[F_{r},C^{\prime},C^{\prime}]$.

  6. 6.

Flattening We flatten the pooled representation to a vector of size
$d_{p}=F_{r}C^{\prime}(C^{\prime}-1)/2$, which corresponds to the number of
free variables in ${\mathbf{S}}$.

  7. 7.

Embedding We reduce the size of the pooled representation vector to $d_{e}$
using a linear layer. The resulting vector is the speaker embedding.

##  4 Experiments

###  4.1 VoxCeleb 1 & 2 datasets

We evaluate the systems on the three VoxCeleb test sets, namely the original
VoxCeleb 1 _test_ set (composed of only 40 speakers) [20] and the (much
larger) _Extended_ and _Hard_ sets defined in [14]. As training set, we use a
subset of the VoxCeleb 2 _development_ set which we augmented with babble,
noise, reverberation, and music using Kaldi [21, 3]. The training set then
contains 5750 speakers and approximately 4.9M utterances including the
augmented versions. As a _held-out_ set used to control the training, we used
455 utterances including augmentations which are not part of the training set,
from 90 speakers that are part of the training set.

###  4.2 Architecture, training, and loss function

The backbone of the network is a 34-layer ResNet. All convolutional kernels
are $3\times 3$, the strides are $[1,2,2,2]$ for the 4 blocks in both time and
frequency axes, while the corresponding number of channels are
$[64,128,256,256]$ (the first convolutional layer also outputs 64 channels).
Finally, the number of convolutional layers per block is [3, 4, 6, 3], which
is the standard configuration of the 34-layer ResNet.

The network is trained using multi-speaker classification. As optimization
criterion we employ the Additive Angular Margin (AAM) loss, which has shown
notable improvements over the plain Cross Entropy loss [22]. Moreover, the
$l_{2}$-norm that applies to both embeddings and speaker representations of
the classification head (i.e. the final linear layer before the Softmax) makes
cosine-similarity the natural choice for measuring the similarity between
enrolment and test utterances. We set the scale coefficient of the loss equal
to 30, while for the margin coefficient we adopt a curriculum learning
approach and progressively increase it as [0.1,0.2,0.3] during training [23].

As optimizer we use stochastic gradient descent with momentum equal to 0.9.
The minibatch size is 256, however in order to train the model with a single
GPU (of 12GB RAM) we split the minibatch into 16 “microbatches” of 16 examples
each and use gradient accumulation. The initial learning rate is equal to 0.2
which we divide by 2 when the loss does not improve for more than 3000 model
updates in the held-out set. The experiments are conducted using TensorFlow
[24].

We emphasize that the model architecture, the training method, and its
hyperparameters are optimized with the baseline, where mean and std features
are used for pooling. Moreover, we did not observe any differences in time
required per epoch between the baseline and the proposed architectures.

Table 1:  Results on the Original, Hard, and Extended VoxCeleb test sets.  | System |  |  |  |  | VoxCeleb-O | VoxCeleb-H | VoxCeleb-E
---|---|---|---|---|---|---|---|---
|  |  |  |  |  | minDCF | EER(%) | minDCF | EER(%) | minDCF | EER(%)
B1 | Baseline | (mean & std pooling) | 0.100 | 1.61 | 0.157 | 2.72 | 0.100 | 1.59
B2 | Baseline | (mean pooling) | 0.103 | 1.63 | 0.160 | 2.81 | 0.103 | 1.63
|  | $C^{\prime}$ | $F_{r}$ | ${\mathbf{L}}$ | Normalization |  |  |  |  |  |
P1 | Proposed | 128 | 1 | 2D | mean | 0.125 | 1.91 | 0.186 | 3.06 | 0.119 | 1.84
P2 | Proposed | 64 | 5 | 2D | mean | 0.115 | 1.71 | 0.169 | 2.90 | 0.110 | 1.76
P3 | Proposed | 64 | 5 | 3D | mean | 0.095 | 1.40 | 0.147 | 2.52 | 0.094 | 1.46
P4 | Proposed | 64 | 10 | 3D | mean | 0.104 | 1.70 | 0.161 | 2.73 | 0.104 | 1.63
P5 | Proposed | 128 | 1 | 2D | mean & var. | 0.095 | 1.52 | 0.161 | 2.74 | 0.104 | 1.63
P6 | Proposed | 64 | 5 | 2D | mean & var. | 0.083 | 1.44 | 0.149 | 2.50 | 0.096 | 1.47
P7 | Proposed | 64 | 5 | 3D | mean & var. | 0.084 | 1.39 | 0.139 | 2.38 | 0.090 | 1.38
P8 | Proposed | 64 | 10 | 3D | mean & var. | 0.094 | 1.48 | 0.145 | 2.50 | 0.092 | 1.49

###  4.3 Experimental Results

To evaluate our method we experiment we examine several configurations and
ablations. The results we obtained are given in Table 1, and are obtained
using cosine-similarity, without score or any other normalization. For the two
baseline experiments we train and evaluate a ResNet-34 with the standard
pooling, one with and one without including std features in the pooling layer.
We observe that std features indeed yields slightly increased performance. We
should emphasize that for the baseline experiments the tensor operations are
not included, i.e. mean and std pooling is performed directly on the output of
the last ResNet block.

The first experiment using a variant of the proposed pooling (P1) is
essentially a direct implementation of the pooling as suggested in image
style-transfer (recall that pooling with $F_{r}=1$ is equivalent to (3), i.e.
to pooling along both time and frequency axes.). Two such networks are
examined, one with mean and a second one with mean and variance normalization,
denoted by P1 and P5, respectively. Note that when mean-only normalization is
applied, the number of free-parameters is $F_{r}C_{r}(C_{r}+1)/2$, i.e. the
size of the pooled representation vector is increased by $F_{r}C_{r}$. By
comparing P1 and P5 we observe that variance normalization helps significantly
the model to attain results that are very close to the baseline B1. We should
also note that for the experiments with $F_{r}=1$, we set $C^{\prime}=128$ in
order to have pooling layers of same order of magnitude, since their size
grows linearly with $F_{r}$.

In the next set of experiments we use either $F_{r}=5$ or $F_{r}=10$, where
the latter indicates no reshaping of the tensor ${\mathbf{Y}}$ (as it
corresponds to $f_{r}=1$). The results confirm that the frequency axis should
be treated differently compared to the time axis, and pooling should result in
a 3D tensor, as shown in (4). They moreover demonstrate that reshaping the
tensor by a factor of $f_{r}=2$ (i.e. $F_{r}=5$) yields the best performance.

We may also observe that the suggested frequency-dependent channel reduction
(i.e. using a 3D tensor ${\mathbf{L}}$ instead of a 2D) yields consistently
improved performance (by comparing P2 with P3, and P6 with P7). Allowing each
frequency range having its own weighted combination of channels seems to be an
effective way to increase the expressiveness of the pooled representation.

We also mention that we attempted to combine our best system (P7) with B1 in
the pooling layer. The improvement in performance attained was small. For
example the EER dropped by 0.03% and the minDCF by 0.002 in both in VoxCeleb-E
and VoxCeleb-H. Although further experimentation is required, it suggests that
combining the correlations-based with mean and std features does not yield any
notable improvement.

Finally, all experiments underline the importance of mean and variance
normalization (P5-8), as opposed to mean-only (P1-4). In other words,
transforming the Gram-Matrices to correlation matrices results in more
speaker-discriminative representations. This finding challenges the common
view of frame aggregation being the natural means of extracting speaker-
relevant information from utterances. Correlations, as inherently normalized
statistical quantities, should be more resilient as features to the diverse
within-speaker session effects, as well as to the effects related to the
variable phonetic content of utterances encountered in text-independent
speaker recognition. Temporal aggregation, plain or attentive, may not be the
most effective way to remove such undesired sources of variability. Removing
them via modeling correlations between different representations of the
utterances (such as time-channel maps of given frequency ranges) may attain
this goal more effectively and naturally.

##  5 Conclusions and research directions

In this paper, we introduced a new way of pooling information extracted from a
CNN for modeling speaker information. The approach is inspired by neural
style-transfer in computer vision, where style is modeled by a Gram-Matrix of
channel-wise dot-products. We proposed a way to adapt this method to speech
signals, e.g. by calculating several such Gram-Matrices, one per frequency
range. We further demonstrated that transforming them to correlation matrices
increases the speaker-discriminability of the extracted representations. Our
experiments on VoxCeleb indicate that the proposed method yields notable
improvements over the standard mean-std pooling.

There are several ways by which our approach can be extended. For estimating
the correlations between channels, the dot-product can be replaced by kernel
functions. The application of kernel functions in image style-transfer is an
active research direction, and some of these methods should be applicable to
our task [9]. The property of Gram-Matrix being a symmetric positive semi-
definite matrix can be explored in several ways, using e.g. statistical
divergences (e.g. Kullback-Leibler, Jensen-Shannon, a.o.). Alternative
training/finetuning schemes (e.g. prototypical, triplet, contrastive and self-
supervised losses [23, 25, 26, 27]) that do not require a classification head
are probably better suited for such kind of experimentation. Attention methods
are also worth-examining, e.g. for assigning data-driven weights in the time
axis before estimating the correlation coefficients. Finally, the whole chain
of transforms we apply on the tensor ${\mathbf{Y}}$ is not a result of
exhaustive experimentation, leaving room for further improvements.

##  6 Acknowledgements

The work was supported by Czech National Science Foundation (GACR) project
“NEUREM” No. 19-26934X, and European Union’s Horizon 2020 project No. 833635
ROXANNE. This project has received funding from the European Union’s Horizon
2020 research and innovation programme under grant agreement No 101007666 /
ESPERANTO / H2020-MSCA-RISE-2020.

## References

  * [1] P. Kenny, G. Boulianne, P. Ouellet, and P. Dumouchel, “Joint factor analysis versus eigenchannels in speaker recognition,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 15, no. 4, pp. 1435–1447, 2007\.
  * [2] N. Dehak, P. J. Kenny, R. Dehak, P. Dumouchel, and P. Ouellet, “Front-end factor analysis for speaker verification,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 19, no. 4, pp. 788–798, 2011.
  * [3] D. Snyder, D. Garcia-Romero, G. Sell, D. Povey, and S. Khudanpur, “X-vectors: Robust DNN embeddings for speaker recognition,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018.
  * [4] T. Zhou, Y. Zhao, and J. Wu, “Resnext and res2net structure for speaker verification,” _arXiv preprint arXiv:2007.02480_ , 2020.
  * [5] K. Okabe, T. Koshinaka, and K. Shinoda, “Attentive statistics pooling for deep speaker embedding,” _Proc. Interspeech 2018_ , pp. 2252–2256, 2018.
  * [6] M. À. India Massana, P. Safari, and F. J. Hernando Pericás, “Self multi-head attention for speaker recognition,” in _Interspeech 2019: the 20th Annual Conference of the International Speech Communication Association: 15-19 September 2019: Graz, Austria_. International Speech Communication Association (ISCA), 2019, pp. 4305–4309.
  * [7] L. You, W. Guo, L.-R. Dai, and J. Du, “Multi-task learning with high-order statistics for x-vector based text-independent speaker verification,” in _Proc. Interspeech 2019_ , 2019, pp. 1158–1162.
  * [8] L. A. Gatys, A. S. Ecker, and M. Bethge, “Image style transfer using convolutional neural networks,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 2414–2423.
  * [9] Y. Jing, Y. Yang, Z. Feng, J. Ye, Y. Yu, and M. Song, “Neural style transfer: A review,” _IEEE transactions on visualization and computer graphics_ , vol. 26, no. 11, pp. 3365–3385, 2019.
  * [10] Y. Li, N. Wang, J. Liu, and X. Hou, “Demystifying neural style transfer,” in _Proceedings of the 26th International Joint Conference on Artificial Intelligence_ , 2017, pp. 2230–2236.
  * [11] J. Chorowski, R. J. Weiss, R. A. Saurous, and S. Bengio, “On using backpropagation for speech texture generation and voice conversion,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 2256–2260.
  * [12] D. Snyder, D. Garcia-Romero, D. Povey, and S. Khudanpur, “Deep neural network embeddings for text-independent speaker verification,” in _INTERSPEECH_ , 2017, pp. 999–1003.
  * [13] P. Safari, M. India, and J. Hernando, “Self-attention encoding and pooling for speaker recognition,” _arXiv preprint arXiv:2008.01077_ , 2020.
  * [14] A. Nagrani, J. S. Chung, W. Xie, and A. Zisserman, “Voxceleb: Large-scale speaker verification in the wild,” _Computer Speech & Language_, vol. 60, p. 101027, 2020.
  * [15] Y. Wu, C. Guo, H. Gao, X. Hou, and J. Xu, “Vector-based attentive pooling for text-independent speaker verification,” _Proc. Interspeech 2020_ , pp. 936–940, 2020.
  * [16] T. Zhou, Y. Zhao, J. Li, Y. Gong, and J. Wu, “CNN with phonetic attention for text-independent speaker verification,” in _2019 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2019, pp. 718–725.
  * [17] S. M. Kye, Y. Kwon, and J. S. Chung, “Cross attentive pooling for speaker verification,” _arXiv preprint arXiv:2008.05983_ , 2020.
  * [18] S. Wang, Y. Yang, Y. Qian, and K. Yu, “Revisiting the statistics pooling layer in deep speaker embedding learning,” in _2021 12th International Symposium on Chinese Spoken Language Processing (ISCSLP)_. IEEE, 2021, pp. 1–5.
  * [19] Y. Gal and Z. Ghahramani, “A theoretically grounded application of dropout in recurrent neural networks,” _Advances in Neural Information Processing Systems_ , pp. 1027–1035, 2016.
  * [20] A. Nagrani, J. S. Chung, and A. Zisserman, “Voxceleb: a large-scale speaker identification dataset,” in _INTERSPEECH_ , 2017, pp. 2616–2620.
  * [21] D. Povey, A. Ghoshal, G. Boulianne, L. Burget, O. Glembek, N. Goel, M. Hannemann, P. Motlicek, Y. Qian, P. Schwarz _et al._ , “The kaldi speech recognition toolkit,” in _IEEE 2011 workshop on automatic speech recognition and understanding_ , no. EPFL-CONF-192584. IEEE Signal Processing Society, 2011.
  * [22] J. Deng, J. Guo, N. Xue, and S. Zafeiriou, “Arcface: Additive angular margin loss for deep face recognition,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 4690–4699.
  * [23] J. S. Chung, J. Huh, S. Mun, M. Lee, H. S. Heo, S. Choe, C. Ham, S. Jung, B.-J. Lee, and I. Han, “In defence of metric learning for speaker recognition,” _arXiv preprint arXiv:2003.11982_ , 2020.
  * [24] H. Zeinali, L. Burget, J. Rohdin, T. Stafylakis, and J. Cernocky, “How to improve your speaker embeddings extractor in generic toolkits,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 6141–6145.
  * [25] C. Zhang and K. Koishida, “End-to-end text-independent speaker verification with triplet loss on short utterances,” in _Proc. Interspeech 2017_ , 2017, pp. 1487–1491.
  * [26] M. Ravanelli and Y. Bengio, “Learning speaker representations with mutual information,” in _Proc. Interspeech 2019_ , 2019, pp. 1153–1157.
  * [27] T. Stafylakis, J. Rohdin, O. Plchot, P. Mizera, and L. Burget, “Self-supervised speaker embeddings,” in _Proc. Interspeech 2019_ , 2019, pp. 2863–2867.


# Scalable subsampling:
computation, aggregation and inference

Dimitris N. Politis
Department of Mathematics
and Halicioglu Data Science Institute
University of California, San Diego
La Jolla, CA 92093-0112, USA
<EMAIL_ADDRESS>
###### Abstract

Subsampling is a general statistical method developed in the 1990s aimed at
estimating the sampling distribution of a statistic $\hat{\theta}_{n}$ in
order to conduct nonparametric inference such as the construction of
confidence intervals and hypothesis tests. Subsampling has seen a resurgence
in the Big Data era where the standard, full-resample size bootstrap can be
infeasible to compute. Nevertheless, even choosing a single random subsample
of size $b$ can be computationally challenging with both $b$ and the sample
size $n$ being very large. In the paper at hand, we show how a set of
appropriately chosen, non-random subsamples can be used to conduct
effective—and computationally feasible—distribution estimation via
subsampling. Further, we show how the same set of subsamples can be used to
yield a procedure for subsampling aggregation—also known as subagging—that is
scalable with big data. Interestingly, the scalable subagging estimator can be
tuned to have the same (or better) rate of convergence as compared to
$\hat{\theta}_{n}$. The paper is concluded by showing how to conduct
inference, e.g., confidence intervals, based on the scalable subagging
estimator instead of the original $\hat{\theta}_{n}$.

Keywords: Bagging, Big Data, bootstrap, confidence intervals, distributed
inference, subagging.

##  1 Introduction

Assume data $X_{1},\ldots,X_{n}$ that are independent, identically distributed
(i.i.d.) taking values in an arbitrary space. Often, this space will be ${\bf
R}^{d}$ but other choices exist, e.g., it can be a function space, a space of
networks, etc. A statistic $\hat{\theta}_{n}=T_{n}(X_{1},\ldots,X_{n})$ is
employed to estimate a parameter $\theta$ associated with the common
distribution of the data. Assume that $\theta$ takes values on a normed linear
space ${\bf\Theta}$ with norm denoted by $||\cdot||$.

Let $J_{n}(x)=P\\{\tau_{n}g(\hat{\theta}_{n}-\theta)\leq x\\}$ where the rate
of convergence $\tau_{n}$ diverges to infinity as $n$ increases. In general,
the function $g(\cdot)$ may be taken to be the norm $||\cdot||$ but other
choices are possible. For example, if ${\bf\Theta}={\bf R}^{p},$ then
$g(\cdot)$ can taken to be the sup-norm on ${\bf R}^{p}$; knowing (or
estimating) the quantiles of $J_{n}$ would then lead to simultaneous
confidence intervals and/or simultaneous hypothesis tests for all $p$
coordinates of $\theta$. Note that if ${\bf\Theta}={\bf R},$ then $g(\cdot)$
can be taken to be the identity function, leading to one-sided inference. We
will also assume:

Assumption A.  There exists a non-degenerate probability distribution $J$,
such that $J_{n}(x)\to J(x)$ as $n\to\infty$ for all $x$ points at which
$J(x)$ is continuous.

Subsampling is a general statistical method developed in the 1990s aimed at
estimating the sampling distribution $J_{n}$ in order to conduct nonparametric
inference such as the construction of confidence intervals and hypothesis
tests; see Politis, Romano and Wolf (1999) and the references therein. To
describe it, let the subsample size $b$ be a positive integer less than $n$,
and consider all the subsets of size $b$ of the sample $X_{1},\ldots,X_{n}$.
There are $Q=$$n\choose b$ such subsets that can be ordered in an arbitrary
fashion and denoted by ${\cal B}_{j}$ for $j=1,\ldots,Q$.

Compute the subsample statistics $\hat{\theta}_{b,j}=T_{b}({\cal B}_{j})$ for
$j=1,\ldots,Q$, and the subsampling distribution

${}L_{n,b}(x)=\frac{1}{Q}\sum_{i=1}^{Q}1\\{\tau_{b}g(\hat{\theta}_{b,i}-\hat{\theta}_{n})\leq
x\\}.$

Under Assumption A and the additional conditions

$b\to\infty\ \ \mbox{as}\ \ n\to\infty\ \ \mbox{but with}\ \ b/n\to 0$ (1)

and

$\tau_{b}/\tau_{n}\to 0$ (2)

it was shown that $L_{n,b}(x)\stackrel{{\scriptstyle
P}}{{\longrightarrow}}J(x)$ for all $x$ points of continuity of $J$, where
$\stackrel{{\scriptstyle P}}{{\longrightarrow}}$ denotes convergence in
probability; see Theorems 3.1 and 3.2 of Politis and Romano (1994). Note that
${}L_{n,b}(x)\stackrel{{\scriptstyle P}}{{\longrightarrow}}J(x)$ implies that
${}L_{n,b}(x)-J_{n}(x)\stackrel{{\scriptstyle P}}{{\longrightarrow}}0$, i.e.,
${}L_{n,b}(x)$ can be used as a proxy for the unknown $J_{n}(x)$ so as to
conduct statistical inference based on $\hat{\theta}_{n}$. Also note that eq.
(2) follows immediately if the rate of convergence satisfies

$\tau_{n}=n^{\alpha}{\cal L}(n)\ \ \mbox{for some}\ \ \alpha>0,$ (3)

and some slowly varying function ${\cal L}(\cdot)$.

It was also recognized early on that if $n$ is large, it is not realistic to
compute $\hat{\theta}_{b,j}$ for all $j=1,\ldots,Q\ $ since $Q$ can be
astronomically large. For that reason, Corollary 2.2 of Politis and Romano
(1994) showed that a stochastic approximation to $L_{n,b}(x)$ can be used
instead. The stochastic approximation relies on $B$ randomly chosen subsamples
from the set $\\{{\cal B}_{j}$, $j=1,\ldots,Q\\}$. Here, $B$ is the
practitioner’s choice; it should be large, tending to infinity as
$n\to\infty$, but it can be chosen small enough for the stochastic
approximation to be computable.

Note that choosing a random subsample from the set $\\{{\cal B}_{j}$,
$j=1,\ldots,Q\\}$ amounts to sampling $b$ values without replacement from the
set $\\{X_{1},\ldots,X_{n}\\}$. This brings analogies to the so-called
bootstrap with smaller resample size where the sampling is with
replacement—cf. Bickel and Freedman (1981), Bretagnolle (1983), and Swanepoel
(1986); see also Politis and Romano (1993), and Bickel, Götze and van Zwet
(1997), and Corollary 2.3.1 of Politis, Romano, and Wolf (1999) for the
comparison of sampling with vs. without replacement.

Subsampling has seen a resurgence in the Big Data era of the 21st century
where Efron’s (1979) full-resample size bootstrap can be infeasible to
compute; see e.g., Jordan (2013), Kleiner et al. (2014), and Sengupta et al.
(2016). Nevertheless, even choosing a single random subsample of size $b$ can
be computationally challenging with both $b$ and $n$ very large; approximate
solutions such as Poisson sampling may be needed as described in Bertail,
Chautru, and Clémencon (2017).

In the next section we show how a set of appropriately chosen, non-random
subsamples can be used to conduct effective—and computationally
feasible—distribution estimation via subsampling. Further, in Section 4 we
show how the same set of subsamples can be used to yield a procedure for
subsampling aggregation—also known as subagging—that is scalable with big data
in an attempt to remedy computability issues discussed in Section 3.
Interestingly, the scalable subagging estimator can be tuned to have the same
(or better) rate of convergence as compared to $\hat{\theta}_{n}$. The paper
is concluded by providing details on how to conduct inference, e.g. confidence
intervals, based on the scalable subagging estimator instead of the original
$\hat{\theta}_{n}$.

##  2 Scalable subsampling distribution estimation

Recall the set of all size $b$ subsamples $\\{{\cal B}_{j}$,
$j=1,\ldots,Q\\}$, and re-arrange it so that the first subsamples are obtained
as blocks of consecutive data points, i.e., ${\cal
B}_{j}=(X_{(j-1)h+1},X_{(j-1)h+2},\ldots,X_{(j-1)h+b})$. Recall that the block
size $b$ is an integer in $[1,n$], and so is $h$; in particular, $h$ controls
the amount of overlap (or separation) between ${\cal B}_{j}$ and ${\cal
B}_{j+1}$. If $h=1$, then the overlap is the maximum possible; if $h\sim 0.2\
b$, then there is an approximate 80% overlap between ${\cal B}_{j}$ and ${\cal
B}_{j+1}$; if $h=b$, then there is $no$ overlap between ${\cal B}_{j}$ and
${\cal B}_{j+1}$ but these two blocks are adjacent; finally, if $h\sim 1.2\
b$, then there is a block of about $0.2\ b$ data points from the data sequence
$X_{1},\ldots,X_{n}$ that separate the blocks ${\cal B}_{j}$ and ${\cal
B}_{j+1}$. Note that, in general, $b$ and $h$ are functions of $n$, but these
dependences will not be explicitly denoted.

The collection of all available block–subsample of size $b$, is then $\\{{\cal
B}_{j},\ j=1,\ldots,q\\}$ where $q=\lfloor(n-b)/h\rfloor+1$ and
$\lfloor\cdot\rfloor$ denotes the integer part function. We claim that this
non-random collection is sufficient for effective and computationally feasible
subsampling distribution estimation. To see why, note that subsampling using
the aforementioned block–subsamples has been found to be consistent in the
setting where the data sequence $X_{1},\ldots,X_{n}$ is a finite stretch of a
strong-mixing, stationary time series; see e.g. Politis and Romano (1994,
Section 3.2). Since the i.i.d. case can be considered as a special case of a
stationary time series, the claim follows.

To elaborate, we define the subsample statistics
$\hat{\theta}_{b,j}=T_{b}({\cal B}_{j})$ for $j=1,\ldots,q$ as before, and
construct the subsampling distribution as

${}L_{n,b,h}(x)=\frac{1}{q}\sum_{i=1}^{q}1\\{\tau_{b}g(\hat{\theta}_{b,i}-\hat{\theta}_{n})\leq
x\\}.$ (4)

###### Proposition 2.1

Assume Assumption A, and conditions (1) and (3). Also assume that either
$h=1$, or that $h$ satisfies

$h\sim c_{1}\ b\ \ \mbox{for some constant}\ \ c_{1}>0.$ (5)

Then, $L_{n,b,h}(x)\stackrel{{\scriptstyle P}}{{\longrightarrow}}J(x)$ for all
$x$ points of continuity of $J$.

The Proposition follows from Corollary 3.2 of Politis and Romano (1994) who
worked under the assumption that $1\leq h\leq b$; the case where $h>b$ —but
still with $h=O(b)$— can be proven in a similar way. The essence of the
argument is that

$EL_{n,b,h}(x)\approx J_{b}(x)\to J(x)\ \ \mbox{as}\ \ b\to\infty$ (6)

where $x$ a point of continuity of $J$. In addition,

$Var(L_{n,b,h}(x))=O(b/n).$ (7)

Eq. (6) and (7) together with Chebyshev’s inequality imply
$L_{n,b,h}(x)\stackrel{{\scriptstyle P}}{{\longrightarrow}}J(x)$.

Note that the bound (7) holds true regardless of the choice of $h$, i.e.,
whether $h=1$ or $h$ satisfies condition (5); it is just the proportionality
constant in $O(b/n)$ that becomes smaller (but is bounded below) as $h$
decreases. Therefore, for reasons of parsimony and computational tractability,
we will not propose using full-overlap block-subsamples, i.e., the case $h=1$,
in what follows. Instead we will work under condition (5), in which case
$q=O(n/b$). Hence, assuming $\hat{\theta}_{n}$ can be computed in
$O(n^{\zeta}$) operations (for some constant $\zeta>0$), the construction of
$L_{n,b,h}$ and its quantiles has computational complexity
$O(n^{\zeta})+O(nb^{\zeta-1})=O(n^{\zeta})$ which is the same as computing the
statistic $\hat{\theta}_{n}$ itself; therefore, we may call the construction
of $L_{n,b,h}$ under condition (5) as being scalable.

In terms of practicalities, the recommendation is: (a) choose an appropriate
block size $b$, e.g. use estimation or calibration as described in Chapter 9
of Politis et al. (1999), or adapt the method of Götze and Rackauskas (2001),
and Bickel and Sakov (2008) to the scalable computation of $L_{n,b,h}(x)$ that
is based on block-subsamples; and (b) choose $h$ small, i.e., choose $c_{1}$
small in condition (5), to make $q$ large (so that the the bound in (7) is
small), but also lying within feasibily/computability constraints.

To give an example, suppose that $n=200,000$, and $b$ is chosen to be 500.
Recall that $L_{n,b,h}(x)$ is a sum of $q$ summands, where
$q=\lfloor(n-b)/h\rfloor+1$. If one chooses $c_{1}=2$, then $q\approx 200$
which is feasible. If a larger $q$ is deemed computationally feasible, one
could take $c_{1}=1$ or even $c_{1}=1/2$ leading to $q$ equalling 400 or 800
respectively.

Implementation problems might ensue only if $\zeta>1$ and $n$ is huge, but in
this case even the computation of the original statistic $\hat{\theta}_{n}$
may be problematic; we will deal with this issue in the next section.

##  3 Computability issues

As discussed in the last section, the statistic $\hat{\theta}_{n}$ will
generally be computable in $O(n^{\zeta}$) operations. If $\zeta$ is small,
then no issues incur. Unfortunately, examples abound with $\zeta>1$, making
the computability of $\hat{\theta}_{n}$ questionable in the Big Data era.

Some examples are as follows:

  1. 1.

The $X_{i}$ are univariate, and $\hat{\theta}_{n}$ is the sample mean (or
median) of $X_{1},\ldots,X_{n}$. Then, $\zeta=1$.

  2. 2.

The $X_{i}$ take values in ${\bf R}^{d}$, and $\hat{\theta}_{n}$ is the sample
mean of $X_{1},\ldots,X_{n}$. Then, $\hat{\theta}_{n}$ is computable in
$O(dn$) operations. If $d$ is a constant, then $\zeta=1$ as above. However, it
may be that $d$ grows with $n$; if $d$ grows linearly with $n$, then
$\zeta=2.$

  3. 3.

Suppose that $X_{i}=(Y_{i},W_{i})$ where $Y_{i}$ is the univariate response
associated with regressor $W_{i}$ that takes values in ${\bf R}^{d}$; this is
the standard regression situation. If $d$ is large, then LASSO regression can
be employed; see Tibshirani (1996). A popular method to compute the LASSO has
computational complexity $O(d^{3}+d^{2}n)$ as long as $d<n$; see Efron et al.
(2004). If $d$ grows linearly with $n$, e.g., when $d\sim n/2$, then
$\zeta=3.$

There is a growing body of work dealing with the possibility that the sample
size $n$ is so large that it may not be feasible to compute
$\hat{\theta}_{n}$. One branch of this literature is devoted to ‘optimal
subsampling’ whose meaning is a little different than the subsampling-based
inference discussed so far. In a nutshell, if $\hat{\theta}_{n}$ is not
computable, one can use just one of the subsample statistics
$\hat{\theta}_{b,i}$ to estimate $\theta$. The question ‘which one to use’ is
addressed by the literature on ‘optimal subsampling’; see Yao and Wang (2021)
for a review.

The problem with this approach is that the practitioner is effectivelly
throwing out the majority of the data. A Divide-and-Conquer alternative is
proposed in the next section.

##  4 Scalable subsampling aggregation

###  4.1 Computation

Ensemble methods are a family of techniques that combine several estimators
(or predictors) with the objective of producing a new estimator or predictor
with better statistical properties. One of the most popular and widely used
ensemble methods is bootstrap aggregating (Breiman, 1996), also known as
bagging; it is a resampling technique whose main purpose is to reduce the
variability of a given estimator/predictor rendering it more stable.

Subsample aggregation, also known as subagging, was proposed by Bühlmann and
Yu (2002). In the context of the present paper, the subagging estimator is
defined as

$\bar{\theta}_{b,SA}=\frac{1}{Q}\sum_{i=1}^{Q}\hat{\theta}_{b,i}$

for an appropriate choice of $b$. Here, and for the remainder of the paper, we
will assume that $\theta$ is univariate, i.e., ${\bf\Theta}={\bf R}.$

Under some regularity conditions, Bühlmann and Yu (2002) showed that

$E\bar{\theta}_{b,SA}=E\hat{\theta}_{b,1},\ \ \mbox{and}\ \
Var(\bar{\theta}_{b,SA})\leq(b/n)Var(\hat{\theta}_{b,1}).$

Hence, if the Bias of $\hat{\theta}_{b,1}$ is tolerable, subagging yields a
welcome variance reduction.

With Big Data it is, of course, infeasible to compute (and average) all
$Q=$$n\choose b$ subsample statistics. As already mentioned, even choosing a
random subset of these $Q$ subsamples presents computational difficulties.
Note, however, that although $\bar{\theta}_{b,SA}$ is an average of $Q$
values, the variance is reduced by dividing by $n/b$ not $Q$; the reason is
that the $Q$ subsamples have typically high overlap, and hence their
associated subsample statistics are highly dependent.

We can achieve the same variance reduction effect by using just the first $q$
subsamples in the ordering described in Section 2, i.e., ${\cal
B}_{j}=(X_{(j-1)h+1},X_{(j-1)h+2},\ldots,X_{(j-1)h+b})$ for $j=1,\ldots,q$.
For better variance reduction —as well as computational savings— we will
choose $h\geq b$ in what follows. We therefore define the scalable subagging
estimator as

$\bar{\theta}_{b,n,SS}=\frac{1}{q}\sum_{i=1}^{q}\hat{\theta}_{b,i}.$

###### Proposition 4.1

Assume $h\geq b$. Then,

$E\bar{\theta}_{b,n,SS}=E\hat{\theta}_{b,1},\ \ \mbox{and}\ \
Var(\bar{\theta}_{b,n,SS})=q^{-1}Var(\hat{\theta}_{b,1})$

where $q=\lfloor(n-b)/h\rfloor+1$.

The proof is immediate noting that
$\hat{\theta}_{b,1},\ldots,\hat{\theta}_{b,q}$ are i.i.d., their independence
being ensured by the assumption $h\geq b$.

###### Remark 4.1

Note that $q=O(n/b)$ since $h\geq b$; hence, $\bar{\theta}_{b,n,SS}$ can be
computed with $O(qb^{\zeta})=O(nb^{\zeta-1})$ operations, a significant
savings over the $O(n^{\zeta})$ operations needed for $\hat{\theta}_{n}$.
Recall the discussion of Section 3 on (i) identifying the ‘optimal subsample’
denoted by ${\cal B}_{j^{*}}$, and then (ii) using as your final estimator the
subsample statistic $\hat{\theta}_{b,j^{*}}$. The computational cost of part
(ii) is of course $O(b^{\zeta})$; if you add to this the nontrivial
computational cost of part (i) —that may require numerical optimization— the
overall complexity may well exceed the $O(nb^{\zeta-1})$ needed to compute the
scalable subagging estimator $\bar{\theta}_{b,n,SS}$.

We can further consider whether $\bar{\theta}_{b,n,SS}$ can have a rate of
convergence that is comparable to that of $\hat{\theta}_{n}$; this is the
subject of the next subsection.

###  4.2 Rate of convergence and choice of $b$ for scalable subagging

If the computational complexity of $\bar{\theta}_{b,n,SS}$ is
$O(nb^{\zeta-1})$, what is to stop us from taking $b$ very small, even $b=1$,
to make it $O(n)$? The answer is the bias of $\hat{\theta}_{b,1}$ that figures
in Proposition 4.1. If the bias of $\hat{\theta}_{b,1}$ is zero (or
negligible), then indeed $b$ could be chosen very small.

Consider the following three Bias Conditions (BC):

  * (I)

Estimator $\hat{\theta}_{n}$ is exactly unbiased, as is the case with linear
statistics; see Example 1.4.1 of Politis et al. (1999). Then, $b$ can be taken
to equal one but, of course, scalable subagging is not needed here as the
computational complexity of $\hat{\theta}_{n}$ is $O(n)$ already.

  * (II)

Estimator $\hat{\theta}_{n}$ is asymptotically unbiased, and its bias is
asymptotically negligible even after multiplication by $\tau_{n}$; in other
words, the limit law $J$ of Assumption A is centered at zero.

  * (III)

Estimator $\hat{\theta}_{n}$ is asymptotically unbiased but its bias does not
vanish after multiplication by $\tau_{n}$; in other words, the limit law $J$
of Assumption A is centered at a nonzero value.

Although the subsampling distribution estimator $L_{n,b,h}$ can work under all
above eventualities —including BC case (III)—, to investigate the rate of
convergence of scalable subagging we will work under the assumption of BC case
(II). Note that an estimator falling under BC case (III) could be analytically
debiased —by subtracting from it a consistent estimate of its bias—, allowing
the debiased estimator to be handled under BC case (II).

Hence, we formulate the following assumption:

Assumption B.  Assume that $E\hat{\theta}_{n}^{2}<\infty$ for all $n$, and
that there exist constants $\gamma>\alpha>0$, $C\in{\bf R}-\\{0\\}$, and
$\sigma^{2}>0$ such that

$n^{\gamma}(E\hat{\theta}_{n}-\theta)\to C,\ \ \mbox{and}\ \
Var(\tau_{n}\hat{\theta}_{n})\to\sigma^{2}\ \ \mbox{as}\ \ n\to\infty$

where $\tau_{n}=n^{\alpha}$.

In the above, we have simplified eq. (3) by omitting the slowly varying
function ${\cal L}(n)$; however, ${\cal L}(n)$ could be incorporated in the
technical arguments that follow if/when so desired. Note that the assumption
$\gamma>\alpha$, implies that the bias of $\hat{\theta}_{n}$ is negligible
even after multiplication by $\tau_{n}$, as in BC case (II) above.

Recall that in Section 2 we operated under the assumption that $h$ is of the
same order of magnitude as $b$. So far, this was not required in the
discussion of scalable subagging; indeed, Proposition 4.1 only needed $h\geq
b$. For the sake of argument, in what follows we will entertain the
possibility that $h$ grows at least as fast as $b$, and maybe faster.
Therefore, we will assume:

$b\sim c_{2}n^{\beta}\ \ \mbox{and}\ \ h\sim c_{3}n^{\delta}\ \ \mbox{as}\ \
n\to\infty\ \ \mbox{for positive constants}\ \ c_{2},c_{3},\ \ \mbox{and
constants}\ \ 0<\beta\leq\delta<1.$ (8)

The following proposition shows that $\bar{\theta}_{b,n,SS}$ can be tuned to
have the same (or better) rate of convergence as compared to
$\hat{\theta}_{n}$. We will use the notation $a_{n}={\Theta}(d_{n})$ to denote
‘exact order’, i.e., that there exist constants $\underline{c},\bar{c}$
satisfying $\underline{c}\cdot\bar{c}>0$, and such that
$\underline{c}d_{n}\leq a_{n}\leq\bar{c}d_{n}.$

###### Proposition 4.2

Assume Assumption B, eq. (8) and $h\geq b$. Also assume $\alpha\leq 1/2$.
Choose a value of $\beta$ satisfying $\beta\geq\alpha/\gamma$, and choose a
value of $\delta$ satisfying

$1-2\beta(\gamma-\alpha)\leq\delta\leq 1+2\alpha(\beta-1).$ (9)

Then:

(i) $MSE(\bar{\theta}_{b,n,SS})=O(\tau_{n}^{-2})$ where MSE is short for Mean
Squared Error.

(ii) If $\delta$ is chosen to equal the upper bound of eq. (8), i.e., if
$\delta=1+2\alpha(\beta-1)$, then
$MSE(\bar{\theta}_{b,n,SS})={\Theta}(\tau_{n}^{-2})$; in this case,
$\bar{\theta}_{b,n,SS}$ has the same rate of convergence as
$\hat{\theta}_{n}$. However, if $\delta<1+2\alpha(\beta-1)$, then
$MSE(\bar{\theta}_{b,n,SS})=o(\tau_{n}^{-2}),$ and $\bar{\theta}_{b,n,SS}$ has
faster rate of convergence than $\hat{\theta}_{n}$.

(iii) If $\delta$ can be chosen to equal the lower bound of eq. (8), i.e., if
$\delta=1-2\beta(\gamma-\alpha)$, then both
$[Bias(\bar{\theta}_{b,n,SS})]^{2}$ and $\ {Var(\bar{\theta}_{b,n,SS})}$ have
the same order of magnitude, namely they are both (asymptotically)
proportional to $n^{-2\beta\gamma}$ with generally different proportionallity
constants; in this case, $\bar{\theta}_{b,n,SS}$ falls under BC case (III).
However, if $\delta>1-2\beta(\gamma-\alpha)$, then
$[Bias(\bar{\theta}_{b,n,SS})]^{2}=o\left({Var(\bar{\theta}_{b,n,SS})}\right)$,
i.e., $\bar{\theta}_{b,n,SS}$ falls under BC case (II).

Proof. First note that $\bar{\theta}_{b,n,SS}$ is the average of
$\hat{\theta}_{b,1},\ldots,\hat{\theta}_{b,q}$ that are i.i.d. due to the fact
that $h\geq b$. Hence, $E\bar{\theta}_{b,n,SS}=E\hat{\theta}_{b,1}$, and
$Bias(\bar{\theta}_{b,n,SS})\sim Cb^{-\gamma})=Cn^{-\beta\gamma}$ by
Assumption B which further implies

$Var(\bar{\theta}_{b,n,SS})=q^{-1}Var(\hat{\theta}_{b,1})\sim\frac{\sigma^{2}}{n^{1-\delta+2\alpha\beta}},$
(10)

since $q=\lfloor(n-b)/h\rfloor+1={\Theta}(n^{1-\delta})$. Therefore,

$MSE(\bar{\theta}_{b,n,SS})\sim
C^{2}n^{-2\beta\gamma}+\frac{\sigma^{2}}{n^{1-\delta+2\alpha\beta}}.$

Choose $\delta\geq 1-2\beta(\gamma-\alpha)$, to ensure $2\beta\gamma\geq
1-\delta+2\alpha\beta$, yielding

$MSE(\bar{\theta}_{b,n,SS})={\Theta}(\frac{1}{n^{1-\delta+2\alpha\beta}}).$

Finally, to have
$\frac{1}{n^{1-\delta+2\alpha\beta}}=O(\frac{1}{n^{2\alpha}})$, we need
$\delta\leq 1+2\alpha(\beta-1)$; with equality, the two rates are the same,
i.e., if $\delta=1+2\alpha(\beta-1)$, then
$\frac{1}{n^{1-\delta+2\alpha\beta}}={\Theta}(\frac{1}{n^{2\alpha}})$.
$\diamond$

###### Remark 4.2

(Choice of $b$.) Elaborating on part (ii) of Proposition 4.2, if we wanted to
maximize the rate of convergence of $\bar{\theta}_{b,n,SS}$ we would operate
at the minimum value of $\delta$, i.e., $\delta=1-2\beta(\gamma-\alpha)$, if
the latter is a possible value for $\delta$. Recall that eq. (8) postulates
$\delta>0$; so in order for $\delta$ to be chosen to equal
$1-2\beta(\gamma-\alpha)$, the latter must be positive which can be ensured
only when

$\beta<\frac{1}{2(\gamma-\alpha)}.$ (11)

We could then choose $\beta$ to further minimize $MSE(\bar{\theta}_{b,n,SS})$
which would entail choosing $\beta$ to be as large as possible; since
$\beta\leq\delta$ from eq. (8), it follows that we would need to take
$\beta=\delta$. Furthermore, solving the equation
$\beta=1-2\beta(\gamma-\alpha)$ implies that the optimal value of $\beta$ is

$\beta=\frac{1}{1+2(\gamma-\alpha)}$ (12)

which is a value satisfying eq. (11), thus allowing $\delta$ to operate at the
minimum value of eq. (9).

###### Corollary 4.1

Under the assumptions of Proposition 4.2, choose a value of $\beta$ satisfying

$\frac{1}{1+2(\gamma-\alpha)}\leq\beta<\frac{1}{2(\gamma-\alpha)},$ (13)

and let $\delta=\beta$. Then,

$MSE(\bar{\theta}_{b,n,SS})={\Theta}(\frac{1}{n^{1-\beta+2\alpha\beta}}).$

Furthermore:

(i) If $\beta>\frac{1}{1+2(\gamma-\alpha)}$, then
$[Bias(\bar{\theta}_{b,n,SS})]^{2}=o\left({Var(\bar{\theta}_{b,n,SS})}\right)$,
i.e., $\bar{\theta}_{b,n,SS}$ falls under BC case (II).

(ii) If $\beta=\frac{1}{1+2(\gamma-\alpha)}$, then
$[Bias(\bar{\theta}_{b,n,SS})]^{2}={\Theta}\left({Var(\bar{\theta}_{b,n,SS})}\right)$,
i.e., $\bar{\theta}_{b,n,SS}$ falls under BC case (III). In addition, the
choice $\beta=\frac{1}{1+2(\gamma-\alpha)}$ minimizes the
$MSE(\bar{\theta}_{b,n,SS})$, yielding

$MSE(\bar{\theta}_{b,n,SS})={\Theta}(n^{-2\gamma/[1+2(\gamma-\alpha)]})$ (14)

in which case the optimized rate of convergence of $\bar{\theta}_{b,n,SS}$ is
$n^{\gamma/[1+2(\gamma-\alpha)]}$.

Note that letting $\beta$ equal the lower bound $\frac{1}{1+2(\gamma-\alpha)}$
kills two birds with one stone: (a) optimizes the rate of convergence of
$\bar{\theta}_{b,n,SS}$, and (b) minimizes the computational complexity in
computing $\bar{\theta}_{b,n,SS}$ making it
$O(nb^{\zeta-1})=O(n^{1+\beta(\zeta-1)})=O\left(n^{1+\frac{(\zeta-1)}{1+2(\gamma-\alpha)}}\right).$

###### Remark 4.3

(Choice of $h$.) It was already establishes that $\delta$ should be taken
equal to its allowed minimum value to optimize $\bar{\theta}_{b,n,SS}$.
Furthermore, looking at eq. (10), it is apparent that to optimize
$Var(\bar{\theta}_{b,n,SS})$ we would need to maximize $q$, i.e., minimize
$h$. It is therefore recommended —when computationally feasible— to choose $h$
at its allowed minimum, i.e., choose $h=b$.

We now discuss some examples:

  1. 1.

Linear statistics. Consider a linear statistic $\hat{\theta}_{n}$ that can be
represented as $\hat{\theta}_{n}=n^{-1}\sum_{i=1}^{n}G(X_{i})$ for some
appropriate function $G$. Note that $\hat{\theta}_{n}$ estimates
$\theta=EG(X_{1})$, and is exactly unbiased for that. As mentioned under the
description of BC case (I), subagging is not needed here because
$\hat{\theta}_{n}$ can be easily computed. Furthermore, this case can not
really fit under the premises of Proposition 4.2 since Assumption B does not
hold; recall that Assumption B implies that $\hat{\theta}_{n}$ has nonzero
bias. However, we can intuite what would happen in this case by pretending
that Assumption B holds (approximately) with a huge value of $\gamma$. Letting
$\gamma\to\infty$, the above discussion implies that the optimal $\beta$ and
$\delta$ should tend to zero. Hence, one should take $b=h=1$, in which case
$\bar{\theta}_{b,n,SS}$ reduces to the original statistic $\hat{\theta}_{n}$.

  2. 2.

Approximately linear statistics. A statistic $\hat{\theta}_{n}$ is called
approximately linear if it can be represented as
$\hat{\theta}_{n}=n^{-1}\sum_{i=1}^{n}G(X_{i})+o_{P}(n^{-1/2})$. Then,
$\hat{\theta}_{n}$ is $\sqrt{n}$–consistent for $\theta=EG(X_{1})$, i.e.,
$\alpha=1/2$. In many interesting examples, e.g., sample quantiles,
$M$-estimators, etc., it may be verified that the bias of $\hat{\theta}_{n}$
is of order $1/n$, i.e., $\gamma=1$. In such a case, the two bounds in eq. (9)
collapse, yielding a single choice for $\delta$. Further optimizing MSE as in
Remark 4.2 implies $\bar{\theta}_{b,n,SS}$ is $\sqrt{n}$–consistent as well,
with the choice $\beta=\delta=1/2$, i.e., block size increasing proportionally
to $\sqrt{n}$.

  3. 3.

Nonparametric function estimators. Consider the case where $\theta$ represents
the value of function $f$ at a point of interest; the function $f$ can be a
probability density, a spectral density, regression function, etc. that should
be estimated in a nonparametric setting; see Rosenblatt (1991) for a unified
description of nonparametric function estimation. Let $\hat{\theta}_{n}$
denote a kernel-smoothed estimator of $\theta$, and suppose that a nonnegative
kernel is used. In this case, the MSE–optimal bandwidth is
${\Theta}(n^{-1/5})$. However, this bandwidth choice brings $\hat{\theta}_{n}$
under the realm of BC case (III), and the premises of Proposition 4.2 do not
apply. As an experiment, consider a degree of undersmoothing in constructing
$\hat{\theta}_{n}$. To fix ideas, suppose the bandwidth is chosen to be
${\Theta}(n^{-1/4})$ instead, yielding

$Bias(\hat{\theta}_{n})=O(\frac{1}{n^{1/2}})\ \ \mbox{and}\ \
Var(\hat{\theta}_{n})={\Theta}(\frac{1}{n^{3/4}});$

in this case, $\alpha=3/8$ and $\gamma=1/2$. According to Remark 4.2 and
Corollary 4.1, $\beta$ and $\delta$ should be optimally chosen to both equal
$0.8$; hence, the rate of convergence of $\bar{\theta}_{b,n,SS}$ becomes
$n^{2/5}$. Interestingly, this rate is not only faster than the rate of
$\hat{\theta}_{n}$ that used the suboptimal bandwidth ${\Theta}(n^{-1/4})$;
the $n^{2/5}$ rate is actually the fastest rate achievable by any kernel-
smoothed estimator that uses a nonnegative kernel with its associated
MSE–optimal bandwidth. Nevertheless, $\bar{\theta}_{b,n,SS}$ can be computed
faster than $\hat{\theta}_{n}$, and is thus preferable.

###  4.3 Inference beyond point estimation

Having established that $\bar{\theta}_{b,n,SS}$ is a consistent estimator
whose rate of convergence towards $\theta$ is fast (and sometimes faster than
that of $\hat{\theta}_{n}$), the question now is how to conduct inference,
e.g., confidence intervals, hypothesis tests, etc. based on
$\bar{\theta}_{b,n,SS}$.

Recall that the assumption $h\geq b$ implies that
$\hat{\theta}_{b,1},\ldots,\hat{\theta}_{b,q}$ are i.i.d. but these values are
not an increasing sequence; to see that, note that the value of
$\hat{\theta}_{b,1}$ changes with $b$ (which increases with $n$). Rather,
$\hat{\theta}_{b,1},\ldots,\hat{\theta}_{b,q}$ can be thought as the $n$th row
of a triangular array with i.i.d. entries, and common distribution given
(approximately, and after centering and standardizing) by $J_{b}$.

Since $\bar{\theta}_{b,n,SS}$ is the sample mean of
$\hat{\theta}_{b,1},\ldots,\hat{\theta}_{b,q}$, a Central Limit Theorem for
triangular arrays —such as Theorem B.O.1 of Politis et al. (1999)— implies the
following useful result.

###### Corollary 4.2

Assume the premises of Corollary 4.1. In addition, assume that there exist
positive numbers $\epsilon$ and $\Delta$ such that
$E|\hat{\theta}_{n}|^{2+\epsilon}\leq\Delta<\infty$ for any $n$. Then,

$n^{\frac{-1+\beta-2\alpha\beta}{2}}\left(\bar{\theta}_{b,n,SS}-\theta\right)\stackrel{{\scriptstyle{\cal
L}}}{{\Longrightarrow}}N(C_{\beta},\sigma^{2})$ (15)

where $\stackrel{{\scriptstyle{\cal L}}}{{\Longrightarrow}}$ denotes
convergence in law as $n\to\infty$. Furthermore:

(i) If $\beta>\frac{1}{1+2(\gamma-\alpha)}$, then $C_{\beta}=0$.

(ii) If $\beta=\frac{1}{1+2(\gamma-\alpha)}$, then $C_{\beta}=C$.

In the above, the constants $C$ and $\sigma^{2}$ are as defined in Assumption
B.

Hence, under case (i) of the above, to construct an estimate of the
distribution of $\bar{\theta}_{b,n,SS}$ all that is needed is a consistent
estimator of $\sigma^{2}$; this is easy to obtain as the sample variance of
the $n$th row of the triangular array, i.e., letting

$\hat{\sigma}^{2}=\frac{b^{2\alpha}}{q}\sum_{i=1}^{q}\left(\hat{\theta}_{b,i}-\bar{\theta}_{b,n,SS}\right)^{2}.$
(16)

For example, an approximate 95% confidence interval for $\theta$ under case
(i) would be:

$\bar{\theta}_{b,n,SS}\pm 1.96\ \hat{\sigma}\cdot
n^{\frac{1-\beta+2\alpha\beta}{2}}.$

Actually, case (ii) of the Corollary 4.2 is more interesting since it ensures
the fastest rate of convergence of $\bar{\theta}_{b,n,SS}$. In this case, the
nontrivial asymptotic bias of the distribution of $\bar{\theta}_{b,n,SS}$
presents some difficulties at first inspection. Nevertheless, subsampling
comes again at the rescue since eq. (15) shows that $\bar{\theta}_{b,n,SS}$
has a well-defined asymptotic distribution; the fact that the latter is not
centered at zero is immaterial. In other words, $\bar{\theta}_{b,n,SS}$
satisfies Assumption A with $\bar{\theta}_{b,n,SS}$ in place of
$\hat{\theta}_{n}$. Hence, the scalable subsampling construction of Section 2
can be applied to yield a consistent estimate of the distribution of
$\bar{\theta}_{b,n,SS}$.

###### Remark 4.4

(Connections with distributed inference.) The case $h=b$, i.e., splitting the
sample into $q$ non-overlapping parts, is closely related to the classical
notion of $q$-fold cross-validation, as well as the more recent notion of
Divide-and-Conquer (DaC) methods; see Jordan (2013). To elaborate, the
scalable subsampling estimator $\bar{\theta}_{b,n,SS}$ has been studied before
in the following DaC distributed inference contexts (all with $h=b$):
U-statistics by Lin and Xi (2010); generalized linear models by Chen and Xie
(2014); M-estimators by Zhang, Duchi and Wainwright (2013); and a certain
class of symmetric statistics (that includes L-statistics and smoothed
functions of the sample mean) by Chen and Peng (2021). Note also that Bradic
(2016) employed subagging using non-overlapping blocks of data, and applied it
to variable selection in large-scale regression. The current section is meant
to serve many purposes. One is to show that these ideas are universally
applicable under minimal assumptions, such as Assumption B; for example, the
asymptotic normality results of Chen and Peng (2021) actually follow from our
Corollary 4.2 simply by checking its premises. Furthermore, it is important to
note that the usefulness of scalable subagging and DaC distributed inference
extends well beyond the realm of asymptotically linear, $\sqrt{n}$–consistent
statistics that have been considered so far; see our Section 4.2 including the
example on nonparametric function estimation. Finally, the interplay of the
tuning parameters $h$ and $b$ opens up interesting possibilities, e.g., the
possibility that $\bar{\theta}_{b,n,SS}$ has a faster rate of convergence than
$\hat{\theta}_{n}$ itself.

###  4.4 Possibly dependent data

All subsampling constructions in this paper, including the scalable
subsampling distribution $L_{n,b,h}(x)$, and scalable subagging estimator
$\bar{\theta}_{b,n,SS}$, remain valid if there is (weak) dependence in the
data, i.e., if $X_{1},\ldots,X_{n}$ are a stretch of a strictly stationary,
strong mixing time series. The reason is that the choice of block-subsamples
described in Section 2 and used throughout the paper is actually the choice
that is recommended in order to subsample time series; see e.g. Politis and
Romano (1994). Hence, all results in this paper remain true as stated in the
case where $X_{1},\ldots,X_{n}$ are weakly dependent with the exception of the
results in Section 4.3 that may require a little tweak.

To see why, note that if $h=b$, then
$\hat{\theta}_{b,1},\ldots,\hat{\theta}_{b,q}$ will be exactly independent
only if the data $X_{1},\ldots,X_{n}$ are independent. If $X_{1},\ldots,X_{n}$
are stationary and strong mixing, we can ensure that
$\hat{\theta}_{b,1},\ldots,\hat{\theta}_{b,q}$ are approximately independent
if we require that $h-b\to\infty$, e.g., letting $h=b+\lfloor\sqrt{b}\rfloor$;
this would ensure that blocks ${\cal B}_{j}$ and ${\cal B}_{j+1}$ are
separated by about $\sqrt{b}\ $ data points, rendering them approximately
independent as $b$ increases with $n$. We summarize this discussion in a
corollary that covers the possibility of dependent data with exponentially
decreasing strong mixing coefficients.

###### Corollary 4.3

Assume the premises of Corollary 4.2. If $h-b\to\infty$, then the results of
Corollary 4.2 remain true even when the data $X_{1},\ldots,X_{n}$ are a
stretch of a strictly stationary, strong mixing time series, with
exponentially decreasing strong mixing coefficients.

Note that condition $h-b\to\infty$ ensures that the large sample variance of
$\bar{\theta}_{b,n,SS}$ is $\sigma^{2}$ as given in eq. (15); the asymptotic
normality of $\bar{\theta}_{b,n,SS}$ (with different variance) would remain
true for any $h\geq b$.

An analog of Corollary 4.3 can also be formulated when the strong mixing
coefficients only decay polynomially fast. In that case, however, the
$\epsilon$ appearing in the moment assumption of Corollary 4.2 can not be any
positive number, i.e., it can not be taken arbitrarily close to zero. Rather,
the minimum value of $\epsilon$ allowed would be dictated by the polynomial
rate of decay of strong mixing; see the assumptions of Theorem B.O.1 of
Politis et al. (1999).

Acknowledgement. Many thanks are due to Ery Arias-Castro, Jelena Bradic and
Yiren Wang for several helpful suggestions. Research partially supported by
NSF grant DMS 19-14556.

## References

  * [1] Bertail, P., Chautru, E. and Clémencon, S. (2017). Empirical Processes in Survey Sampling with (Conditional) Poisson Designs, Scandinavian Journal of Statistics, 44: 97-111.
  * [2] Bickel, P. J. and Freedman, D. A. (1981). Some asymptotic theory for the bootstrap. Ann. Statist., 9: 1196-1217.
  * [3] Bickel, P.J., Götze F., and van Zwet W. R. (1997). Resampling Fewer Than n Observations: Gains, Losses, and Remedies for Losses, Statistica Sinica, 7: 1-31.
  * [4] Bickel, P. J. and Sakov, A. (2008). On the choice of m in the m out of n bootstrap and confidence bounds for extrema, Statistica Sinica, 18(3): 967-985.
  * [5] Bradic, J. (2016). Randomized maximum-contrast selection: Subagging for large-scale regression. Elect. J. Statist., 10: 121-170.
  * [6] Breiman, L. (1996). Bagging predictors. Machine Learning, 24: 123-140.
  * [7] Bretagnolle, J. (1983). Lois limites du bootstrap de certaines fonctionelles, Ann. Inst. H. Poincaré, Sec. B (N.S.), 3: 281-296.
  * [8] Bühlmann, P. and Yu, B. (2002). Analyzing bagging. Ann. Statist., 30(4): 927-961.
  * [9] Chen, S.X. and Peng, L. (2021). Distributed statistical inference for massive data. Ann. Statist. 49(5): 2851-2869.
  * [10] Chen, X. and Xie, M. (2014). A split-and-conquer approach for analysis of extraordinarily large data. Statist. Sinica, 24 1655–1684.
  * [11] Efron, B. (1979). Bootstrap methods: another look at the jackknife. Ann. Statist. 7: 1-26.
  * [12] Efron, B., Hastie, T., Johnstone, J. and Tibshirani, R. (2004). Least angle regression. Ann. Statist., 32(2): 407-499.
  * [13] Götze, F. and Rackauskas, A. (2001). Adaptive choice of bootstrap sample sizes. In State of the Art in Probability and Statistics, pp. 286-309, IMS publications, Hayward, California.
  * [14] Jordan, M.J. (2013). On statistics, computation and scalability, Bernoulli, 19(4): 1378-1390.
  * [15] Kleiner, A., Talwalkar, A., Sarkar, P. and Jordan, M.J. (2014). A scalable bootstrap for massive data. J. Roy. Soc., Ser. B, 76(4): 795-816.
  * [16] Lin, N. and Xi, R. (2010). Fast surrogates of U-statistics. Comput. Statist. Data Anal. 54 16-24.
  * [17] McElroy, T.S. and Politis, D.N. (2020). Time Series: A First Course with Bootstrap Starter, Chapman and Hall/CRC Press, Boca Raton.
  * [18] Politis, D.N., and Romano, J.P. (1993). Estimating the Distribution of a Studentized Statistic by Subsampling, In Bulletin of the International Statistical Institute, 49th Session, Firenze, August 25 - September 2, 1993, Book 2, pp. 315-316.
  * [19] Politis, D.N., and Romano, J.P. (1994), Large sample confidence regions based on subsamples under minimal assumptions, Ann. Statist., 22: 2031-2050.
  * [20] Politis, D.N., Romano, J.P. and Wolf, M. (1999), Subsampling, Springer, New York.
  * [21] Rosenblatt, M. (1991). Stochastic curve estimation, NSF-CMPS Regional Conference Series in Probability and Statistics (vol. 3), Institute of Mathematical Statistics, Hayward, California.
  * [22] Sengupta, S., Volgushev, S. and Shao, X. (2016). A subsampled double bootstrap for massive data, Journal of the American Statistical Association, 111(515): 1222-1232.
  * [23] Swanepoel, J.W.H. (1986). A note on proving that the (modified) bootstrap works. Comm. Statist. Theory Methods, 15: 3193-3203.
  * [24] Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. J. Roy. Statist. Soc. Ser. B., 58: 267–288.
  * [25] Yao, Y. and Wang, H. (2021). A review on optimal subsampling methods for massive datasets. Journal of Data Science, 19(1): 1–22.
  * [26] Zhang, Y., Duchi, J. C. and Wainwright, M. J. (2013). Communication-efficient algorithms for statistical optimization. J. Mach. Learn. Res. 14 3321–3363.