Year: 2021

Abstract:
Experiments with pre-trained models such as BERT are often based on a single checkpoint. While the conclusions drawn apply to the artifact tested in the experiment (i.e., the particular instance of the model), it is not always clear whether they hold for the more general procedure which includes the architecture, training data, initialization scheme, and loss function. Recent work has shown that repeating the pre-training process can lead to substantially different performance, suggesting that an alternative strategy is needed to make principled statements about procedures. To enable researchers to draw more robust conclusions, we introduce MultiBERTs, a set of 25 BERT-Base checkpoints, trained with similar hyper-parameters as the original BERT model but differing in random weight initialization and shuffling of training data. We also define the Multi-Bootstrap, a non-parametric bootstrap method for statistical inference designed for settings where there are multiple pre-trained models and limited test data. To illustrate our approach, we present a case study of gender bias in coreference resolution, in which the Multi-Bootstrap lets us measure effects that may not be detected with a single checkpoint. The models and statistical library are available online, along with an additional set of 140 intermediate checkpoints captured during pretraining to facilitate research on learning dynamics.

text: Experiments with pre-trained models such as BERT are often based on a single checkpoint. While the conclusions drawn apply to the artifact tested in the experiment (i.e., the particular instance of the model), it is not always clear whether they hold for the more general procedure which includes the architecture, training data, initialization scheme, and loss function. Recent work has shown that repeating the pre-training process can lead to substantially different performance, suggesting that an alternative strategy is needed to make principled statements about procedures. To enable researchers to draw more robust conclusions, we introduce MultiBERTs, a set of 25 BERT-Base checkpoints, trained with similar hyper-parameters as the original BERT model but differing in random weight initialization and shuffling of training data. We also define the Multi-Bootstrap, a non-parametric bootstrap method for statistical inference designed for settings where there are multiple pre-trained models and limited test data. To illustrate our approach, we present a case study of gender bias in coreference resolution, in which the Multi-Bootstrap lets us measure effects that may not be detected with a single checkpoint. The models and statistical library are available online,1 along with an additional set of 140 intermediate checkpoints captured during pretraining to facilitate research on learning dynamics.

heading: 1 INTRODUCTION
text: Contemporary natural language processing (NLP) relies heavily on pretrained language models, which are trained using large-scale unlabeled data (Bommasani et al., 2021). BERT (Devlin et al., 2019) is a particularly popular choice: it has been widely adopted in academia and industry, and aspects of its performance have been reported on in thousands of research papers (see, e.g., Rogers et al., 2020, for an overview). Because pre-training large language models is computationally expensive (Strubell et al., 2019), researchers often rely on the release of model checkpoints through libraries such as HuggingFace Transformers (Wolf et al., 2020), which enable them to use large-scale language models without repeating the pre-training work. Consequently, most published results are based on a small number of publicly released model checkpoints.
While this reuse of model checkpoints has lowered the cost of research and facilitated head-to-head comparisons, it limits our ability to draw general scientific conclusions about the performance of a particular class of models (Dror et al., 2019; D’Amour et al., 2020; Zhong et al., 2021). The key issue is that reusing model checkpoints makes it hard to generalize observations about the behavior of a single model artifact to statements about the underlying pre-training procedure which created it. Pre-training such models is an inherently stochastic process which depends on the initialization of the model’s parameters and the ordering of training examples; for example, D’Amour et al. (2020) report substantial quantitative differences across multiple checkpoints of the same model architecture on several “stress tests” (Naik et al., 2018; McCoy et al., 2019). It is therefore difficult to know how much of the success of a model based on the original BERT checkpoint is due to BERT’s design, and how much is due to idiosyncracies of a particular artifact. Understanding this difference is critical if we are to generate reusable insights about deep learning for NLP, and improve the state-of-the-art going forward (Zhou et al., 2020; Dodge et al., 2020; Aribandi et al., 2021).
1http://[anonymized]
This paper describes MultiBERTs, an effort to facilitate more robust research on the BERT model. Our primary contributions are:
• We release MultiBERTs, a set of 25 BERT-Base checkpoints to facilitate studies of robustness to parameter initialization and order of training examples (§2). Releasing these models preserves the benefits to the community of a single checkpoint release (i.e., low cost of experiments, apples-to-apples comparisons between studies based on these checkpoints), while enabling researchers to draw more general conclusions about the BERT pre-training procedure.
• We present the Multi-Bootstrap, a non-parametric method to quantify the uncertainty of experimental results based on multiple pre-training seeds (§3), and provide recommendations for how to use the Multi-Bootstrap and MultiBERTs in typical experimental scenarios. We implement these recommendations in a software library.
• We illustrate the approach with a practical use case: we investigate the impact of counterfactual data augmentation on gender bias, in a BERT-based coreference resolution systems (Webster et al., 2020) (§4). Additional examples are provided in Appendix E, where we document challenges with reproducing the widely-used original BERT checkpoint.
The release also includes an additional 140 intermediate checkpoints, captured during training for 5 of the runs (28 checkpoints per run), to facilitate studies of learning dynamics. Our checkpoints and statistical libraries are available at: http://[anonymized].
Additional Related Work. The MultiBERTs release builds on top of a large body of work that seeks to analyze the behavior of BERT (Rogers et al., 2020). In addition to the studies of robustness cited above, several authors have introduced methods to reduce BERT’s variability during finetuning (Zhang et al., 2021; Mosbach et al., 2021; Dodge et al., 2020; Lee et al., 2020; Phang et al., 2018). Other authors have also studied the time dimension, which motivates our release of intermediate checkpoints (Liu et al., 2021; Hao et al., 2020; Saphra & Lopez, 2019; Chiang et al., 2020; Dodge et al., 2020). Similarly to §3, authors in the NLP literature have recommended best practices for statistical testing (Koehn, 2004; Dror et al., 2018; Berg-Kirkpatrick et al., 2012; Card et al., 2020; Søgaard et al., 2014; Peyrard et al., 2021), many of which are based on existing tests to estimate the uncertainty of test sample. In concurrent work, Deutsch et al. (2021) considered bootstrapping methods similar to the Multi-Bootstrap, in the context of summarization metrics evaluation. Also in concurrent work, the Mistral (Karamcheti et al., 2021) project released a set of 10 GPT-2 models with intermediate checkpoints at different stages of pre-training. Our work is complementary, focusing on BERT, introducing a larger number of pre-training seeds, and presenting a methodology to draw robust conclusions about model performance.

heading: 2 RELEASE DESCRIPTION
text: We first describe the MultiBERTs release: how the checkpoints were trained and how their performance compares to the original BERT on two common language understanding benchmarks.

heading: 2.1 TRAINING
text: Overview. The MultiBERTs checkpoints are trained following the code and procedure of Devlin et al. (2019), with minor hyperparameter modifications necessary to obtain comparable results on GLUE (Wang et al., 2019); a detailed discussion of these differences is provided in Appendix E. We use the BERT-Base, Uncased architecture with 12 layers and embedding size 768. We trained the models on a combination of BooksCorpus (Zhu et al., 2015) and English Wikipedia. Since the exact dataset used to train the original BERT is not available, we used a more recent version that was collected by Turc et al. (2019) with the same methodology.
Checkpoints. We release 25 models trained for two million steps each, each training step involving a batch of 256 sequences. For five of these models, we release 28 additional checkpoints captured over the course of pre-training (every 20,000 training steps up to 200,000, then every 100,000 steps). In total, we release 165 checkpoints, about 68 GB of data.
Training Details. As in the original BERT paper, we used batch size 256 and the Adam optimizer (Kingma & Ba, 2014) with learning rate 1e-4 and 10,000 warm-up steps. We used the default values for all the other parameters, except the number of steps, which we set to two million, and sequence length, which we set to 512 from the beginning with up to 80 masked tokens per sequence.2 We follow the BERT code and initialize the layer parameters from a truncated Normal distribution, using mean 0 and standard deviation 0.02. We train using the same configuration as Devlin et al. (2019)3, with each run taking about 4.5 days on 16 Cloud TPU v2 chips.
Environmental Impact Statement. We estimate compute costs at around 1728 TPU-hours for each pre-training run, and around 208 GPU-hours plus 8 TPU-hours for associated fine-tuning experiments (§2.2, including hyperparameter search and 5x replication). Using the calculations of Lacoste et al. (2019)4, we estimate this as about 250 kg CO2e for each of our 25 models. Counting the 25 runs each of CDA-incr and CDA-full from §4, associated coreference models (20 GPU-hours per pretraining model), and additional experiments of Appendix E, this gives a total of about 12.0 metric tons CO2e before accounting for offsets or clean energy. Based on the report by Patterson et al. (2021) of 78% carbon-free energy in Google Iowa (us-central1), we estimate that reproducing these experiments would emit closer to 2.6 tons CO2e, or slightly more than two passengers on a round-trip flight between San Francisco and New York.
By releasing the trained checkpoints publicly, we aim to enable many research efforts on reproducibility and robustness without requiring this cost to be incurred for every subsequent study.
2.2 PERFORMANCE BENCHMARKS
GLUE Setup. We report results on the development sets of the GLUE tasks: CoLA (Warstadt et al., 2019), MNLI (matched) (Williams et al., 2018), MRPC (Dolan & Brockett, 2005), QNLI (v2) (Rajpurkar et al., 2016; Wang et al., 2019), QQP (Chen et al., 2018), RTE (Bentivogli et al., 2009), SST2 (Socher et al., 2013), and SST-B (Cer et al., 2017). In all cases we follow the same approach as Devlin et al. (2019). For each task, we fine-tune BERT for 3 epochs using a batch size of 32. We run a parameter sweep on learning rates [5e-5, 4e-5, 3e-5, 2e-5] and report the best score. We run the procedure five times for each of the 25 models and average the results.
2Specifically, we keep the sequence length constant (the paper uses 128 tokens for 90% of the training then 512 for the remaining 10%) to expose the model to more tokens and simplify the implementation. As we were not able to reproduce original BERT exactly using either 1M or 2M steps (see Appendix E for discussion), we release MultiBERTs trained with 2M steps under the assumption that higher-performing models are more interesting objects of study.
3We train with TensorFlow (Abadi et al., 2015) version 2.5 in v1 compatibility mode. 4https://mlco2.github.io/impact/
SQuAD Setup. We report results on the development sets of SQuAD versions 1.1 and 2.0 (Rajpurkar et al., 2016; 2018), using a setup similar to that of Devlin et al. (2019). For both sets of experiments, we use batch size 48, learning rate 5e-5, and train for 2 epochs.
Results. Figures 1 and 2 show the distribution of the MultiBERTs models’ performance on the development sets of GLUE and SQuAD, in comparison to the original BERT checkpoint.5 On most tasks, original BERT’s performance falls within the same range as MultiBERTs (i.e., original BERT is between the minimum and maximum of the MultiBERTs’ scores). However, original BERT outperforms all MultiBERTs models on QQP, and under-performs them on SQuAD. The discrepancies may be explained by both randomness and differences in training setups, as investigated further in Appendix E.
To further illustrate the performance variability inherent to pre-training and fine-tuning, we analyze the instance-level agreement between the models in Appendix C.

heading: 3 HYPOTHESIS TESTING USING MULTIPLE CHECKPOINTS
text: The previous section compared MultiBERTs with the original BERT, finding many similarities but also some differences (e.g., in the case of SQuAD). To what extent can these results be explained by random noise? More generally, how can we quantify the uncertainty of a set of experimental results when there are multiple sources of randomness?
In parallel to the MultiBERTs release, we propose a more principled and standardized method to compare training procedures. We recommend a non-parametric bootstrapping procedure, the “Multi-Bootstrap”, which enables us to make inference about model performance in the face of multiple sources of uncertainty: the randomness due to the pre-training seed, the fine-tuning seed, and the finite test data. The main idea is to use the average behavior over seeds as a means of summarizing expected behavior in an ideal world with infinite samples.
Although we present Multi-Bootstrap in the context of BERT analysis, the method could be applied in all setups that involve a set of checkpoints pre-trained with the same method, a finite test set, and (possibly) multiple rounds of fine-tuning. The Multi-Bootstrap is implemented as a library alongside the MultiBERTs release.

heading: 3.1 INTERPRETING STATISTICAL RESULTS
text: The Multi-Bootstrap provides an estimate of the amount of remaining uncertainty when summarizing the performance over multiple seeds. The following notation will help us state this precisely. We assume access to model predictions f(x) for each instance x in the evaluation set. We consider randomness arising from:
1. The choice of pre-training seed S ∼ M 2. The choice of fine-tuning seed T ∼ N 3. The choice of test sample (X,Y ) ∼ D
The Multi-Bootstrap procedure allows us to account for all of the above. Specifically, MultiBERTs enables us to estimate the variance due to the choice of pre-training seed (1), which would not be possible with a single artifact. Note that multiple fine-tuning runs are not required in order to use the procedure.
For each pre-training seed s, let fs(x) denote the learned model’s prediction on input features x and let L(s) denote the expected performance metric of fs on a test distribution D over features X and labels Y . For example, the accuracy would be L(s) = E[1{Y = fs(X)}]. We can use the test sample (which we will assume has nx examples) to estimate the performance for each of the seeds in MultiBERTs, which we denote as L̂(s).
5We used https://storage.googleapis.com/bert_models/2020_02_20/uncased_ L-12_H-768_A-12.zip, as linked from https://github.com/google-research/bert.
The performance L(s) depends on the seed, but we are interested in summarizing the model over all seeds. A natural summary is the average over seeds, ES∼M [L(S)], which we will denote by θ. Then, using ns independently sampled seeds, we can compute an estimate θ̂ as
θ̂ = 1
ns
ns∑
j=1
L̂(Sj) .
Because θ̂ is computed under a finite evaluation set and finite number of seeds, it is necessary to quantify the uncertainty of the estimate. The goal of Multi-Bootstrap is to estimate the distribution of the error in this estimate, θ̂ − θ, in order to compute confidence intervals and test hypotheses about θ, such as whether it is above some threshold of interest. Below, we describe a few common experimental designs in NLP that can be studied with these tools.
Design 1: Comparison to a Fixed Baseline. In many use cases, we want to compare BERT’s behavior to that of a single, fixed baseline. For instance, does BERT encode information about syntax as a feature-engineered model would (Tenney et al., 2019; Hewitt & Manning, 2019)? Does it encode social stereotypes, and how does it compare to human biases (Nadeem et al., 2020)? Does it encode world knowledge, similarly to explicit knowledge bases (Petroni et al., 2019)? Does another model such as RoBERTa (Liu et al., 2019) outperform BERT on common tasks such as those from the GLUE benchmark?
In all these cases, we compare MultiBERTs to some external baseline of which we only have a single estimate (e.g., random or human performance), or against an existing model that is not derived from the MultiBERTs checkpoints. We treat the baseline as fixed, and assess only the uncertainty that arises from MultiBERTs’ random seeds and the test examples.
Design 2: Paired Samples. Alternatively, we might seek to assess the effectiveness of a specific intervention on model behavior. In such studies, an intervention is proposed (e.g., representation learning via a specific intermediate task, or a specific architecture change) which can be applied to any pre-trained BERT checkpoint. The question is whether the procedure results in an improvement over the original BERT pre-training method: does the intervention reliably produce the desired effect, or is the observed effect due to the idiosyncracies of a particular model artifact? Examples of such studies include: Does intermediate tuning on NLI after pre-training make models more robust across language understanding tasks (Phang et al., 2018)? Does pruning attention heads degrade model performance on downstream tasks (Voita et al., 2019)? Does augmenting BERT with information about semantic roles improve performance on benchmark tasks (Zhang et al., 2020)?
We refer to studies like the above as paired since each instance of the baseline model fs (which does not receive the intervention) can be paired with an instance of the proposed model f ′s (which receives the stated intervention) such that fs and f ′ s are based on the same pretrained checkpoint produced using the same seed. Denoting θf and θf ′ as the expected performance defined above for the baseline and intervention model respectively, our goal is to test hypotheses about the true difference in performance δ = θf ′ − θf using the estimated difference δ̂ = θ̂f ′ − θ̂f .
In a paired study, Multi-Bootstrap allows us to estimate both of the errors θ̂f − θf and θ̂f ′ − θf ′ , as well as the correlation between the two. Together, these allow us to approximate the distribution of the overall estimation error δ̂ − δ = (θ̂f − θ̂f ′) − (θf − θf ′), between the estimate δ̂ and the truth δ. With this, we can compute confidence intervals for δ, the true average effect of the intervention on performance over seeds, and test hypotheses about δ, as well.
Design 3: Unpaired Samples. Finally, we might seek to compare a number of seeds for both the intervention and baseline models, but may not expect them to be aligned in their dependence on the seed. For example, the second model may use a different architecture so that they cannot be built from the same checkpoints, or the models may be generated from entirely separate initialization schemes. We refer to such studies as unpaired. Like in a paired study, the Multi-Bootstrap allows us to estimate the errors θ̂f − θf and θ̂f ′ − θf ′ ; however, in an unpaired study, we cannot estimate the correlation between the errors. Thus, we assume that the correlation is zero. This will give a conservative estimate of the error (θ̂f − θ̂f ′) − (θf − θf ′), as long as θ̂f − θf and θ̂f ′ − θf ′ are
not negatively correlated. Since there is little reason to believe that the random seeds used for two different models would induce a negative correlation between the models’ performance, we take this assumption to be relatively safe.
Hypothesis Testing. Given the measured uncertainty, we recommend testing whether or not the difference is meaningfully different from some arbitrary predefined threshold (i.e., 0 in the typical case). Specifically, we are often interested in rejecting the null hypothesis that the intervention does not improve over the baseline model, i.e.,
H0 : δ ≤ 0 (1) in a statistically rigorous way. This can be done with the Multi-Bootstrap procedure described below.

heading: 3.2 MULTI-BOOTSTRAP PROCEDURE
text: The Multi-Bootstrap is a non-parametric bootstrapping procedure that allows us to estimate the distribution of the error θ̂ − θ over the seeds and test instances. The algorithm supports both paired and unpaired study designs, differentiating the two settings only in the way the sampling is performed.
To keep the presentation simple, we will assume that the performance L(s) is an average of a perexample metric ℓ(x, y, fs) over the distribution D of (X,Y ), such as accuracy or the log likelihood, and L̂(s) is similarly an empirical average with the observed nx test examples,
L(s) = ED[ℓ(X,Y, fs)], and L̂(s) = 1
nx
nx∑
i=1
ℓ(Xi, Yi, fs).s
We note that the mapping D 7→ L(s) is linear in D, which is required for our result in Theorem 1. However, we conjecture that this is an artifact of the proof; like most bootstrap methods, the method here likely generalizes to any performance metric which behaves asymptotically like a linear mapping of D, including AUC, BLEU score (Papineni et al., 2002), and expected calibration error.
Building on the rich literature on bootstrap methods (e.g., Efron & Tibshirani, 1994), the MultiBootstrap is a new procedure which accounts for the way that the combined randomness from the seeds and test set creates error in the estimate θ̂. The statistical underpinnings of this approach have theoretical and methodological connections to inference procedures for two-sample tests (Van der Vaart, 2000), where the samples from each population are independent. However, in those settings, the test statistics naturally differ as a result of the scientific question at hand.
In our procedure, we generate a bootstrap sample from the full sample with replacement separately over both the randomness from the pre-training seed s and from the test set (X,Y ). That is, we generate a sample of pre-training seeds (S∗1 , S ∗ 2 , . . . , S ∗ ns) with each S ∗ j drawn randomly with replacement from the pre-training seeds, and we generate a test set sample ((X∗1 , Y ∗ 1 ), (X ∗ 2 , Y ∗ 2 ), . . . , (X ∗ nx , Y ∗ nx)) with each (X,Y ) pair drawn randomly with replacement from the full test set. Then, we compute the bootstrap estimate θ̂∗ as
θ̂∗ = 1
ns
ns∑
j=1
L̂(S∗j ), where L̂ ∗(s) =
1
nx
nx∑
i=1
ℓ(X∗i , Y ∗ i , fs).
To illustrate the procedure, we present a minimal Python implementation in Appendix A. For sufficiently large nx and ns, the distribution of the estimation error θ̂ − θ is approximated well by the distribution of θ̂∗ − θ̂ over re-draws of the bootstrap samples, as stated precisely in Theorem 1. Theorem 1. Assume that E[ℓ2(X,Y, fS)] < ∞. Furthermore, assume that for each s, E[ℓ2(X,Y, fs)] < ∞, and for almost every (x, y) pair, E[ℓ2(X,Y, fS) | X = x, Y = y] < ∞. Let n = nx +ns, and assume that 0 < ps = ns/n < 1 stays fixed (up to rounding error) as n → ∞. Then, there exists 0 < σ2 < ∞ such that √n(θ̂ − θ) d→ G with G ∼ N (0, σ2). Furthermore, conditionally on ((X1, Y1), (X2, Y2), . . . ), √ n(θ̂∗ − θ̂) d→ G.
The proof of Theorem 1 is in Appendix B, along with a comment on the rate of convergence for the approximation error. The challenge with applying existing theory to our method is that while the
seeds and data points are each marginally iid, the observed losses depend on both, and therefore are not iid. Therefore, we need to handle this non-iid structure in our method and proof.
For nested sources of randomness (e.g., if for each pre-training seed s, we have estimates from multiple finetuning seeds), we average over all of the inner samples (finetuning seeds) in every bootstrap sample, motivated by Field & Welsh (2007)’s recommendations for bootstrapping clustered data.
Paired Samples (design 2, continued). In a paired design, the Multi-Bootstrap procedure can additionally tell us the joint distribution of θ̂f ′ − θf ′ and θ̂f − θf . To do so, one must use the same bootstrap samples of the seeds (S∗1 , S ∗ 2 , . . . , S ∗ ns) and test examples ((X∗1 , Y ∗ 1 ), (X ∗ 2 , Y ∗ 2 ), . . . , (X ∗ nx , Y ∗ nx)) for both models. Then, the correlation between the errors θ̂f ′ − θf ′ and θ̂f − θf is well approximated by the correlation between the bootstrap errors θ̂∗f ′ − θ∗f ′ and θ̂∗f − θ∗f . In particular, recall that we defined the difference in performance between the intervention f ′ and the baseline f to be δ, and defined its estimator to be δ̂. With the Multi-Bootstrap, we can estimate the bootstrapped difference δ̂∗ = θ̂∗f ′ − θ̂∗f . With this, the distribution of the estimation error δ̂ − δ is well approximated by the distribution of δ̂∗ − δ̂ over bootstrap samples.
Unpaired Samples (design 3, continued). For studies that do not match the paired format, we adapt the Multi-Bootstrap procedure so that, instead of sampling a single pre-training seed that is shared between f and f ′, we sample pre-training seeds for each one independently. The remainder of the algorithm proceeds as in the paired case. Relative to the paired design discussed above, this additionally assumes that the errors due to differences in pre-training seed between θ̂f ′ − θf ′ and θ̂f − θf are independent.
Comparison to a Fixed Baseline (design 1, continued). Often, we do not have access to multiple estimates of L(s), for example, when the baseline f against which we are comparing is an estimate of human performance for which only mean accuracy was reported, or when f is the performance of a previously-published model for which there only exists a single artifact or for which we do not
have direct access to model predictions. When we have only a point estimate θ̂f = L̂(S1) of θf for the baseline f with a single seed S1, we recommend using Multi-Bootstrap to compute a confidence interval around θf ′ and reporting where the given estimate of baseline performance falls within that distribution. An example of this is given in Figure 1, in which the distribution of MultiBERTs performance is compared to that from the single checkpoint of the original BERT release. In general such results should be interpreted conservatively, as we cannot make any claims about the variance of the baseline model.
Hypothesis Testing. A valid p-value for the hypothesis test described in Equation 1 is the fraction
of bootstrap samples from the above procedure for which the estimate δ̂ is negative.

heading: 4 APPLICATION: GENDER BIAS IN COREFERENCE SYSTEMS
text: We present a case study to illustrate how Multi BERTs and the Multi-Bootstrap can help us draw more robust conclusions about model behavior. The use case is based on a particular property of BERT: gendered correlations. For a particular measure of gender bias, we take a single BERT checkpoint and measure a value of 0.35. We then apply an intervention, foo, designed to reduce this correlation, and measure 0.25. In an effort to do even better, we create a whole new checkpoint by applying the foo procedure from the very beginning of pre-training. On this checkpoint, we measure 0.3. How does one make sense of this result?
As a concrete example, we analyze gender bias in coreference systems (Rudinger et al., 2018) and showing how MultiBERTs and the Multi-Bootstrap can help us understand the effect of our intervention, counterfactual data augmentation (CDA). We follow a set-up similar to Webster et al.
(2020), which augments the BERT pretraining data with counterfactual sentences created by randomly swapping binary-gendered pronouns. The goal is to weaken the correlation between gendered pronouns and other words such as occupation terms (e.g., doctor, nurse). We compare our baseline MultiBERTs models to two strategies for CDA. In the first (CDA-incr), we continue pre-training each MultiBERTs model for an additional 50K steps on the counterfactual data of Webster et al. (2020). In the second, we train BERT models from scratch (CDA-full) on the same dataset.
The Winogender dataset consists of template sentences covering 60 occupation terms and instantiated with either male, female, or neutral pronouns. We follow Webster et al. (2020) and train a gold-mention coreference system using a two-layer feedforward network that takes span representations from a frozen BERT encoder as input and makes binary predictions for mention-referent pairs. The model is trained on OntoNotes (Hovy et al., 2006) and evaluated on the Winogender examples for both per-sentence accuracy and a bias score, defined as the Pearson correlation between the peroccupation bias score (Figure 4 of Rudinger et al. 2018) and the occupational gender statistics from the U.S. Bureau of Labor Statistics.6 For each pre-training run, we train five coreference models, using the same encoder but different random seeds to initialize the classifier weights and to shuffle the training data.

heading: 4.1 PAIRED ANALYSIS: CDA-INCR VS. BASE
text: We investigate the impact of the intervention on performance and bias. Overall accuracy is fairly consistent across pre-training seeds, at 62.6±1.2% for the base model, with only a small and not statistically significant change under CDA-incr (Table 1). However, as shown in Figure 3, there is considerable variation in bias correlation, with r values between 0.2 and 0.8 depending on pretraining seed.7 The range for CDA-incr overlaps somewhat, with values between 0.0 and 0.4; however, because the incremental CDA is an intervention on each base checkpoint, we can look at the individual seeds and see that in most cases there appears to be a significant improvement. A paired Multi-Bootstrap allows us to quantify this and further account for noise due to the finite evaluation
6We use the occupation data as distributed with the Winogender dataset, https://github.com/ rudinger/winogender-schemas.
7Some of this variation is due to the classifier training, but on this task there is a large intrinsic contribution from the pretraining seed. See Appendix D for a detailed analysis.
sample of 60 occupations. The results are shown in Table 1, which show that CDA-incr significantly reduces bias by δ̂ = −0.162 with p = 0.001.

heading: 4.2 UNPAIRED ANALYSIS: CDA-FULL VS. CDA-INCR
text: We can also test if we get any additional benefit from running the entire pre-training with counterfactually-augmented data. Similar to MultiBERTs, we trained 25 CDA-full checkpoints for 2M steps on the CDA dataset.8 Because these are entirely new checkpoints, independent from the base MultiBERTs runs, we use an unpaired version of the Multi-Bootstrap, which uses the same set of examples but independently samples seeds from the CDA-incr and CDA-full runs. As shown in Table 2, overall accuracy does not change appreciably (0.622 vs. 0.623, p = 0.416), while bias correlation seems to decrease but not significantly (0.256 vs 0.192, δ = -0.064 with p = 0.132).
As an ablation, we also experiment with sampling over either only seeds (taking the set of examples, i.e. occupations, as fixed), or over examples (taking the set of 25 seeds as fixed). As shown in Table 2, we find lower p-values (0.005 and 0.053) in both cases—showing that failing to account for finite samples along either dimension could lead to overconfident conclusions.
In Appendix E, we present two additional examples: a paired study where we increase pretraining time from 1M to 2M steps, as well as an unpaired comparison to the original bert-base-uncased checkpoint.

heading: 5 CONCLUSION
text: To make progress on language model pre-training, it is essential to distinguish between the properties of specific model artifacts and those of the training procedures that generated them. To this end, we have presented two resources: MultiBERTs, a set of 25 model checkpoints to support robust research on BERT, and the Multi-Bootstrap, a non-parametric statistical method to estimate the uncertainty of model comparisons across multiple training seeds. We demonstrated the utility of these resources by showing how to quantify the effect of an intervention to reduce a type of gender bias in coreference systems built on BERT. We hope that the release of multiple checkpoints and the use of principled hypothesis testing will become standard practices in research on pre-trained language models.
8Following Webster et al. (2020), we use 20 masks per sequence instead of the 80 from Devlin et al. (2019).

heading: A MINIMAL IMPLEMENTATION OF THE MULTI-BOOTSTRAP
text: Below, we present a simplified Python implementation of the Multi-Bootstrap algorithm presented in Section 3.2. It describes a single-sided version of the procedure, which could be used, e.g., to test that a model’s performance is greater than 0. The input is a matrix of predictions where row indices correspond to test examples and column indices to random seeds. The functions returns an array of
nboot samples [θ̂1, . . . , θ̂nboot ].
1 def multibootstrap(predictions, labels, metric_fun, nboot): 2 """ 3 Generates bootstrap samples of a model’s performance.
4 5 Input: 6 predictions: 2D Numpy array with the predictions for different seeds. 7 labels: 1D Numpy array with the labels. 8 metric_fun: Python function. Takes a pair of arrays as input, and
returns a metric or loss.
9 nboot: Number of bootstrap samples to generate.
10 11 Output: 12 Numpy array with nboot samples.
13 14 """ 15 # Checks the data format. 16 n_samples, n_seeds = predictions.shape 17 assert labels.shape == (n_samples,)
18 19 thetas = np.zeros(nboot) 20 for boot_ix in range(nboot): 21 # Samples n_samples test examples and n_seeds pre-training seeds. 22 x_samples = np.random.choice(n_samples, size=n_samples, replace=True) 23 s_samples = np.random.choice(n_seeds, size=n_seeds, replace=True)
24 25 # Computes the metric over the bootstrapping samples. 26 sampled_predictions = predictions[np.ix_(x_samples, s_samples)] 27 sampled_labels = labels[x_samples] 28 sampled_metrics = [ 29 metric_fun(sampled_predictions[:,j], sampled_labels) 30 for j in range(n_seeds) 31 ]
32 33 # Averages over the random seeds. 34 thetas[boot_ix] = np.mean(sampled_metrics)
35 36 return thetas
We provide the complete version of the algorithm on our repository http://[anonymized]. Our implementation is optimized and supports all the experiment designs described in Section 3, including paired and unpaired analysis as well as multiple fine-tuning runs for each pretraining seed.

heading: B PROOF OF THEOREM 1
text: Before giving the proof, we define some useful notation that will simplify the argument considerably. We let Dn be the empirical measure over the nx observations (Zi = (Xi, Yi)) n i=1, and Mn be the empirical measure over the ns observations (Sj) n j=1. For a function f : V → R and a distribution P over V , we will use the shorthand Pf to denote the expectation of f under P ,
Pf = EV∼P [f(V )].
For example, this allows us to write
θ = DMℓ = EZ∼DES∼M ℓ(Z, fS), and θ̂ = DnMnℓ = 1
nx
nx∑
i=1
1
ns
ns∑
j=1
ℓ(Zi, fSj ).
For the bootstrapped distributions, let D∗n denote the distribution over the bootstrap data samples (Z∗1 , Z ∗ 2 , . . . , Z ∗ nx) and M ∗ n denote the distribution over the bootstrapped seed samples, (S∗1 , S ∗ 2 , . . . , S ∗ ns), both conditional on the observed samples (Zi) nx i=1 and (Sj) ns j=1. Note that the empirical average over a bootstrapped sample
1
nx
nx∑
i=1
1
ns
ns∑
j=1
ℓ(Z∗i , fS∗j )
can be written as 1
nx
nx∑
i=1
1
ns
ns∑
j=1
AiBjℓ(Zi, fSj ),
where Ai is the number of times Zi appears in the bootstrapped sample (Z ∗ k) nx k=1, and Bj is the number of times Sj appears in the bootstrapped sample (S ∗ k) ns k=1. With this in mind, we will abuse notation, and also denote D∗n as the distribution over the Ai and M ∗ n as the distribution over the Bj . Finally, we will use E∗ and Var∗ to denote the expectation and variance of random variables defined with respect to D∗n or M ∗ n, conditional on Dn and Mn. We will use P to denote the distribution P = D×M . Throughout, all assertions made with respect to random variables made without a note about their probability of occurrence hold P -almost surely.
Proof. The challenge with applying existing theory to our method is that because the performance metric (ℓ(Zi, fSj ) nx i=1 over the nx observations for a given seed Sj all depend on the same Sj , they are not independent. Similarly for the performance on a given observation, over seeds. Therefore, we need to handle this non-iid structure in our proof for the multi-bootstrap.
There are conceptually three steps to our proof that allow us to do just that. The first is to show that
θ̂ has an asymptotically linear representation as √ n(θ̂ − θ) = √n(Dn −D)Mℓ+ √ n(Mn −M)Dℓ+ oP (1). (2)
The second is to show that conditional on Dn and Mn the multi-bootstrapped statistic θ̂ ∗ ∆= D∗nM ∗ nℓ has an asymptotically linear representation as √ n(θ̂∗ − θ̂) = √n(D◦n −Dn)Mℓ+ √ n(M◦n −Mn)Dℓ+ oP∗(1), (3) where D◦n and M ◦ n are multiplier bootstrap samples coupled to the bootstrap D ∗ n and M ∗ n which we define formally in the beginning of Step 2. The third step is to use standard results for the multiplier bootstrap of the mean of iid data to show that the distributions of the above linearized statistics converge to the same limit.
Because we have assumed that ℓ(Z, fS) < ∞, E[ℓ(Z, fS) | S] < ∞, and E[ℓ(Z, fS) | Z] < ∞, Fubini’s theorem allows us to switch the order of integration over Z and S as needed.
We will assume that DMℓ(X,Y, fS) = 0. This is without loss of generality, because adding and subtracting √ nDMℓ to the bootstrap expression gives
√ n(θ̂∗ − θ̂) = √n(D∗nM∗nℓ−DnMnℓ)
= √ n(D∗nM ∗ nℓ−DMℓ+DMℓ−DnMnℓ) = √ n(D∗nM ∗ n(ℓ−DMℓ)−DnMn(ℓ−DMℓ)),
so if we prove that the result holds with the mean zero assumption, it will imply that the result holds for ℓ with a nonzero mean.
This theorem guarantees consistency of the Multi-Bootstrap estimates. One question that comes up is whether it is possible to get meaningful / tight rates of convergence for the approximation. Unfortunately, getting OP (1/n) convergence as found in many bootstrap methods (Van der Vaart, 2000) is difficult without the use of Edgeworth expansions, by which the Multi-Bootstrap is not welladapted to analysis. That said, many of the remainder terms already have variance of order O(1/n), or could easily be adapted to the same, suggesting an OP (1/ √ n) convergence. The main difficulty, however, is showing rates of convergence for the strong law on separately exchangeable arrays (see the proof of Lemmas 2, 4-5). Showing a weaker notion of convergence, such as in probability, may perhaps allow one to show that the remainder is OP (1/ √ n), however the adaptation of the aforementioned Lemmas is nontrivial.
Step 1 Recalling that θ̂ ∆ = DnMnℓ and θ ∆ = DMℓ, we can expand √ n(θ̂ − θ) as follows, √ n(DnMnℓ−DMℓ) = √ n(DnMnℓ−DMnℓ+DMnℓ−DMℓ)
= √ n((Dn −D)Mnℓ+D(Mn −M)ℓ) = √ n((Dn −D)Mnℓ+ (Dn −D)Mℓ− (Dn −D)Mℓ+D(Mn −M)ℓ) = √ n((Dn −D)Mℓ+ (Dn −D)(Mn −M)ℓ+D(Mn −M)ℓ)
The following lemma shows that √ n(Dn −D)(Mn −M)ℓ is a lower order term.
Lemma 1. Under the assumptions of Theorem 1, √ n(Dn −D)(Mn −M)ℓ = oP (1).
Therefore,
√ n(DnMnℓ−DMℓ) =
1√ 1− ps √ nx(Dn −D)Mℓ+ 1√ ps √ ns(Mn −M)Dℓ+ oP (1).
Step 2 One of the challenges with working with the bootstrap sample D∗n and M ∗ n is that the induced per-sample weights {Ai}nxi=1 and {Bj}nsj=1 do not have independent components, because they each follow a multinomial distribution over nx items and ns items, respectively. However, they are close enough to independent that we can define a coupled set of random variables {A◦i }nxi=1 and {B◦j }nsj=1 that do have independent components, but behave similarly enough to {Ai} and {Bj} that using these weights has a negligible effect on distribution of the bootstrapped estimator, as described concretely below.
First, we discuss the coupled multiplier bootstrap sample D◦n and M ◦ n. The creation of this sequence, called “Poissonization” is a standard technique for proving results about the empirical bootstrap that require independence of the bootstrap weights (van der Vaart et al., 1996). We describe this for D◦n as the idea is identical for M◦n. Because our goal is to couple this distribution to D ∗ n, we define it on the same sample space, and extend the distribution P ∗, expectation E∗ and variance Var∗ to be over D◦n and M ◦ n, conditionally on Dn and Mn, as with D ∗ n and M ∗ n. To construct the distribution D◦n, from the empirical distribution Dn and a bootstrap sample D ∗ n, start with the distribution D∗n and modify it as follows: We draw a Poisson random variable Nnx with mean nx. If Nnx > nx, then we sample Nnx −nx iid observations from Dn, with replacement, and add them to the bootstrap sample initialized with D∗n to produce the distribution D ◦ n. If Nnx < nx, we sample nx − Nnx observations from D∗n, without replacement, and remove them from the bootstrap sample to produce the distribution D◦n. If Nnx = nx, then D ◦ n = D ∗ n. Recalling that Ai is the number of times the i-th sample is included in D ∗ n, similarly define A ◦ i as the number of times the i-th sample is included in D◦n. Note that by the properties of the Poisson distribution, A◦i ∼ Poisson(1), and {A◦i }nxi=1 are independent. Note that the natural normalization for D◦n would be Nnx . However, it will be useful to maintain the normalization by nx, so abusing notation, for a function f(z), we will say that D◦nf = 1 nx ∑nx i=1 A ◦ i f(Zi).
Define θ̂◦ as the following empirical estimator of θ under the distribution D◦n ×M◦n,
θ̂◦ = D◦nM ◦ nℓ =
1
nx
nx∑
i=1
1
ns
ns∑
j=1
A◦iB ◦ j ℓ(Zi, fSj ).
Lemma 2 shows that √ n(θ̂∗ − θ̂◦) = oP∗(1), and so √ n(θ̂∗ − θ) = √n(θ̂◦ − θ) + oP∗(1). Lemma 2. Under the assumptions of Theorem 1, and that DMℓ = 0, √ n(θ̂∗ − θ̂◦) = oP∗(1).
With this, the expansion of √ n(θ̂◦ − θ̂) begins mutatis mutandis the same as in Step 1, to get that
√ n(θ̂◦ − θ̂) = 1√ 1− ps √ nx(D ◦ n −Dn)Mnℓ+ √ n(D◦n −Dn)(M◦n −Mn)ℓ
+ 1√ ps √ ns(M ◦ n −Mn)Dnℓ.
As with Step 1, we provide Lemma 3 showing that the remainder term √ n(D◦n −Dn)(M◦n −Mn)ℓ will be lower order. Lemma 3. Under the assumptions of Theorem 1, √ n(D◦n −Dn)(M◦n −Mn)ℓ = oP∗(1).
Therefore, √ n(D◦nM ◦ nℓ−DnMnℓ) =
1√ 1− ps √ nx(D ◦ n −Dn)Mnℓ+ 1√ ps √ ns(M ◦ n −Mn)Dnℓ+ oP∗(1).
Then, to write √ n(θ̂∗−θ̂) in terms of √ns(M◦n−Mn)Dℓ as wanted in Eq. (3), instead of √ ns(M ◦ n− Mn)Dnℓ, we must additionally show that the functional has enough continuity that the error term√ ns(M ◦ n −Mn)(Dn −D)ℓ is lower order. The following lemma shows exactly this.
Lemma 4. Under the assumptions of Theorem 1, conditionally on the sequences Z1, Z2, . . . and S1, S2, . . . ,
(a) √ n(D◦n −Dn)(Mn −M)ℓ = oP∗(1), and (b) √ n(Dn −D)(M◦n −Mn)ℓ = oP∗(1).
Altogether, these imply that
√ n(D∗nM ∗ nℓ−DnMnℓ) =
1√ 1− ps √ nx(D ◦ n −Dn)Mℓ+ 1√ ps √ ns(M ◦ n −Mn)Dℓ+ oP∗(1).
Step 3 Noting that Mℓ(·, fS) = ED×M [ℓ(·, fS) | Z = ·] is a real-valued random variable with finite variance (similarly for Dℓ(Z, ·)), and recalling that the nx samples used for Dn and ns samples for Mn satisfy n = nx/(1 − ps) and n = ns/ps, for 0 < ps < 1, the conventional central limit theorem shows that for some positive semi-definite matrix Σ ∈ R2×2, and G ∼ N (0,Σ),
√ n ( (Dn −D)Mℓ (Mn −M)Dℓ ) =
( 1
1−ps √ nx(Dn −D)Mℓ
1 ps
√ ns(Mn −M)Dℓ
) d→ G.
Note that Dn and Mn are independent, so G is, in fact, a diagonal matrix.
Additionally, the conditional multiplier CLT (van der Vaart et al., 1996, Lemma 2.9.5, pg. 181) implies that conditionally on Z1, Z2, . . . and S1, S2, . . . ,
√ n ( (D∗n −Dn)Mℓ (M∗n −Mn)Dℓ ) d→ G.
Finally, applying the delta method (see Theorem 23.5 from Van der Vaart (2000)) along with the results from Steps 1 and 2 shows that the distributions of √ n(θ̂ − θ) and √n(θ̂∗ − θ̂) converge to N (0, σ2), where σ2 = Σ11/(1− ps) + Σ22/ps.

heading: B.1 PROOF OF LEMMA 1
text: Fix ǫ > 0. Note that E[(Dn −D)(Mn −M)ℓ] = 0, so by Chebyshev’s inequality,
P ( |√n(Dn −D)(Mn −M)ℓ| > ǫ ) ≤ Var(
√ n(Dn −D)(Mn −M)ℓ)
ǫ2 .
Therefore, it suffices to show that limn→∞ Var( √ n(Dn−D)(Mn−M)ℓ) = 0. To do so, we apply the law of total variance, conditioning on Dn, and bound the resulting expression by C/n. Var( √ n(Dn −D)(Mn −M)ℓ) = nE[Var((Dn −D)(Mn −M)ℓ | Dn)]
+ nVar(E[(Dn −D)(Mn −M)ℓ | Dn]) = nE[Var((Dn −D)(Mn −M)ℓ | Dn)] = nE[Var((Mn −M)(Dn −D)ℓ | Dn)]
= E   n n2s ns∑
j=1
Var((Dn −D)ℓ(·, fSj ) | Dn)
 
= E
[ n
ns Var((Dn −D)ℓ(·, fS1) | Dn)
]
= E   1
ps E
    1 nx nx∑
i=1
ℓ(Zi, fS1)− E[ℓ(Zi, fS1) | S1]
  2
| {Zi}nxi=1    
= E   1
ps   1 nx nx∑
i=1
ℓ(Zi, fS1)− E[ℓ(Zi, fS1) | S1]
  2  
= E   1 psn2x nx∑
i=1
nx∑
k=1
(ℓ(Zi, fS1)− E[ℓ(Zi, fS1) | S1])(ℓ(Zk, fS1)− E[ℓ(Zk, fS1) | S1])
 
= E   1 psn2x nx∑
i=1
(ℓ(Zi, fS1)− E[ℓ(Zi, fS1) | S1])2  
= 1
ps(1− ps)n E
[ (ℓ(Z1, fS1)− E[ℓ(Z1, fS1) | S1])2 ] ≤ C
n → 0.

heading: B.2 PROOF OF LEMMA 2
text: First, note the following representation for θ̂∗ − θ̂◦:
θ̂∗ − θ̂◦ = 1 nx
nx∑
i=1
1
ns
ns∑
j=1
AiBjℓ(Zi, fSj )− 1
nx
nx∑
i=1
1
ns
ns∑
j=1
A◦iB ◦ j ℓ(Zi, fSj )
= 1
ns
ns∑
j=1
(Bj −B◦j ) nx nx∑
i=1
Aiℓ(Zi, fSj )
︸ ︷︷ ︸ ∆ =I1
+ 1
nx
nx∑
i=1
(Ai −A◦i ) ns ns∑
j=1
B◦j ℓ(Zi, fSj )
︸ ︷︷ ︸ ∆ =I2
.
Let ǫ > 0. Noting that E∗[I1] = E∗[I2] = 0, applying Chebyshev’s inequality gives
P ∗ (√ n|θ̂∗ − θ̂◦| > ǫ ) ≤ nVar ∗(θ̂∗ − θ̂◦) ǫ2 ≤ 2nVar ∗(I1) + Var ∗(I2) ǫ2
It suffices to show that nVar∗(I1) → 0 and nVar∗(I2) → 0. The arguments for each term are mutatis mutandis the same, and so we proceed by showing the proof for I2.
By the law of total variance,
Var∗(I2) = Var ∗(E∗[I2 | {Bj}nsj=1]) + E∗[Var∗(I2 | {Bj}nsj=1)].
Because E∗[Ai] = E∗[A◦i ] and {Bj}nsj=1 ⊥ Ai, A◦i , it follows that E∗[I2 | {Bj}nsj=1] = 0. Taking the remaining term and re-organizing the sums in I2,
Var∗(I2) = E ∗  Var ∗   1
nx
nx∑
i=1
(Ai −A◦i )   1 ns ns∑
j=1
Bjℓ(Zi, fSj )   | {Bj}nsj=1     . (4)
Next, we apply the law of total variance again, conditioning on Nnx = ∑ i A ◦ i . First,
E ∗[I2 | Nnx , {Bj}nsj=1] = Nnx − nx nx 1 nx
nx∑
i=1
1
ns
ns∑
j=1
Bjℓ(Zi, fSj ),
and so
Var∗ ( E ∗[I2 | Nnx , {Bj}nsj=1] | {Bj}nsj=1 ) = 1
nx   1 nx nx∑
i=1
1
ns
ns∑
j=1
Bjℓ(Zi, fSj )
  2
Then, conditionally on Nnx (and {Bj}), I2 is the (centered) empirical average of |Nn − n| samples from a finite population of size n, rescaled by |Nn − n|/n. Therefore, applying Theorem 2.2 of Cochran (2007) gives the conditional variance as
|Nnx − nx| n2x
  1
nx − 1
nx∑
i=1
  1 ns ns∑
j=1
Bjℓ(Zi, fSj )
  2
− nx nx − 1   1 nx nx∑
i=1
1
ns
ns∑
j=1
Bjℓ(Zi, fSj )
  2  
︸ ︷︷ ︸ ∆ =V 2
.
To take the expectation over Nnx , notice that because E ∗[Nnx ] = nx, this is the mean absolute deviation (MAD) of Nnx . Using the expression for the MAD of a Poisson variable from Ramasubban (1958) gives
E ∗|Nnx − nx| = 2nx nnxx exp(−nx) nx! ,
and using Stirling’s approximation, this is bounded by C √ nx, for some 0 < C < ∞.
Combining this with the above term for the variance of the conditional expectation, we have
Var∗   1
nx
nx∑
i=1
(Ai −A◦i )   1 ns ns∑
j=1
Bjℓ(Zi, fSj )   | {Bj}nsj=1   ≤ 1
nx   1 nx nx∑
i=1
1
ns
ns∑
j=1
Bjℓ(Zi, fSj )
  2
+ 1
n1.5x V 2. (5)
Noting that E∗[B2j ] = E ∗[BjBk] = 1, we get the following bound:
Var∗(I2) ≤ 1
nx   1 nx nx∑
i=1
1
ns
ns∑
j=1
ℓ(Zi, fSj )
  2 + 1
n1.5x V̄ 2,
where
V̄ 2 = 1
nx − 1
nx∑
i=1
  1 ns ns∑
j=1
ℓ(Zi, fSj )
  2
− nx nx − 1   1 nx nx∑
i=1
1
ns
ns∑
j=1
ℓ(Zi, fSj )
  2
.
Because of the assumption that DMℓ = 0, the SLLN adapted to separately exchangeable arrays (Rieders, 1991, Theorem 1.4) implies that
lim n→∞
1
nx
nx∑
i=1
1
ns
ns∑
j=1
ℓ(Zi, fSj ) = 0,
almost surely. Therefore, the first term of (5) is o(1/n).
Note that V̄ 2 is the empirical variance of the conditional expectation of ℓ(Zi, fSj ) given {Zi}ni=1. Therefore, the law of total variance shows that
V̄ 2 ≤ 1 nx 1 ns
nx∑
i=1
ns∑
j=1
ℓ2(Zi, fSj )−   1 nx 1 ns nx∑
i=1
ns∑
j=1
ℓ(Zi, fSj )
  2
.
By the SLLN adapted to separately exchangeable arrays (Rieders, 1991, Theorem 1.4), both of the terms converge almost surely to DMℓ2 < ∞ and (DMℓ)2, respectively. and therefore,
lim n→∞ nVar∗(Is) ≤ lim n→∞
n
nx   1 nx nx∑
i=1
1
ns
ns∑
j=1
ℓ(Zi, fSj )
  2 + n
n1.5x V̄ 2 = 0.

heading: B.3 PROOF OF LEMMA 3
text: As with Lemma 1, the main idea of the proof is to apply Chebyshev’s inequality, and show that the variance tends to zero. Indeed, choosing an arbitrary ǫ > 0,
P ∗ ( |√n(D◦n −Dn)(M◦n −Mn)ℓ| ≥ ǫ ) ≤ Var
∗ (√n(D◦n −Dn)(M◦n −Mn)ℓ )
ǫ2 .
Therefore, it suffices to show that the variance in the above display goes to zero. To do this, we start by re-writing the expression in terms of A◦i and B ◦ j , and then apply the law of total variance.
Var∗ (√ n(D◦n −Dn)(M◦n −Mn)ℓ ) = nVar∗   1 nxns nx∑
i=1
ns∑
j=1
(A◦i − 1)(B◦j − 1)ℓ(Zi, fSj )
 
= nVar∗  E∗   1 nxns nx∑
i=1
ns∑
j=1
(A◦i − 1)(B◦j − 1)ℓ(Zi, fSj ) | {A◦i }nxi=1
   
+ nE∗  Var∗   1 nxns nx∑
i=1
ns∑
j=1
(A◦i − 1)(B◦j − 1)ℓ(Zi, fSj ) | {A◦i }nxi=1
    .
Because {B◦j }nsj=1 are independent of {A◦i }nxi=1, and have mean 1, the conditional expectation in the first term is 0 almost surely. Expanding out the second term, using that Var∗(B◦j ) = 1, and that the {B◦j }nsj=1 are uncorrelated,
nE∗  Var∗   1 nxns nx∑
i=1
ns∑
j=1
(A◦i − 1)(B◦j − 1)ℓ(Zi, fSj ) | {Ai}nxi=1
   
= nE∗   1
n2s
ns∑
j=1
Var∗  (B◦j − 1) 1
nx
nx∑
i=1
(A◦i − 1)ℓ(Zi, fSj ) | {A◦i }nxi=1
   
= nE∗   1
n2s
ns∑
j=1
  1 nx nx∑
i=1
(A◦i − 1)ℓ(Zi, fSj )
  2  
= nE∗   1 n2s ns∑
j=1
1
n2x
nx∑
i=1
nx∑
k=1
(A◦i − 1)(A◦k − 1)ℓ(Zi, fSj )ℓ(Zk, fSj )
  .
Now, noting that Var∗(A◦i ) = 1, and that the {A◦i }nxi=1 are uncorrelated, this simplifies to
nE∗   1 n2s ns∑
j=1
1
n2x
nx∑
i=1
(A◦i − 1)2ℓ2(Zi, fSj )
  = n
nsnx
1
ns
ns∑
j=1
1
nx
nx∑
i=1
ℓ2(Zi, fSj ).
Because ED×M [ℓ2(Z, fS)] < ∞, the SLLN adapted to separately exchangeable arrays (Rieders, 1991, Theorem 1.4) implies that this converges almost surely to 0.

heading: B.4 PROOF OF LEMMA 4
text: We prove (a) of the Lemma, as (b) follows from applying Fubini’s theorem and following mutatis mutandis the same argument. Without loss of generality, we will assume that ℓ(Zi, fSj ) ≥ 0. Because Var(ℓ(Zi, fSj )) < ∞, we can always decompose ℓ(·, ·) into a positive and negative part, and show that the result holds for each individually.
Once again, we prove (a) by turning to Chebyshev’s inequality. Fix ǫ > 0, and observe that
P ∗ ( |√n(D◦n −Dn)(Mn −M)ℓ| > ǫ ) ≤ Var
∗ (√n(D◦n −Dn)(Mn −M) )
ǫ2 ,
so it is sufficient to show that Var∗ (√ n(D◦n −Dn)(Mn −M) ) → 0. Writing the above in terms of A◦i , we have
Var∗ (√ n(D◦n −Dn)(Mn −M) )
= Var∗   √ n
nx
nx∑
i=1
(A◦i − 1)   1 ns ns∑
j=1
ℓ(Zi, fSj )− E[ℓ(Zi, fSj ) | Zi]
   
= n
n2x
nx∑
i=1
Var∗ (A◦i − 1)   1 ns ns∑
j=1
ℓ(Zi, fSj )− E[ℓ(Zi, fSj ) | Zi]
  2
= n
n2x
nx∑
i=1
  1 ns ns∑
j=1
ℓ(Zi, fSj )− E[ℓ(Zi, fSj ) | Zi]
  2
.
Now, we want to show that the last display converges almost surely to 0. Notice that each term within the outer sum will obviously converge due to the SLLN. Showing that the outer sum also converges almost surely is technically difficult, but conceptually follows the same argument used to prove the SLLN (specifically, we follow the one done elegantly by Etemadi (1981); Luzia (2018) provides a more detailed account of this proof technique that is helpful for developing a deeper understanding).
We show the following version of almost sure convergence: that for any ǫ > 0,
P   n
n2x
nx∑
i=1
  1 ns ns∑
j=1
ℓ(Zi, fSj )− E[ℓ(Zi, fSj ) | Sj ]
  2
> ǫ i.o.   = 0,
where i.o. stands for infinitely often.
Define the shorthand Lij = ℓ(Zi, fSj ) and let L̄ij = Lij1{Lij < ij} be a truncated version of Lij . The proof of Theorem 2 of Etemadi (1981) implies that P (L̄ij 6= Lij i.o.) = 0, because the assumption Var(Lij) < ∞ implies the assumption used in Etemadi (1981), and independence of {Lij}i,j is not needed for this result. Therefore,
1
nx
nx∑
i=1
  1 ns ns∑
j=1
Lij − L̄ij
  2
a.s.→ 0, and 1 nx
nx∑
i=1
  1 ns ns∑
j=1
E[Lij | Zi]− E[L̄ij | Zi]
  2
a.s.→ 0.
Together, these imply that if we can prove that the truncated sum converges, ie.,
1
nx
n∑
i=1
  1 ns ns∑
j=1
L̄ij − E[L̄ij | Zi]
  2
a.s.→ 0, (6)
this is sufficient to show that the un-truncated version converges almost surely.
To prove (6), we show two things: first, that there is a subsequence kn such that (6) holds when restricted to the subsequence, and then we show that the sequence is a Cauchy sequence, which together imply the result.
Let α > 1 and let kn = α n. For convenience, denote knx as the number of data samples and kns as the number of seed samples when knx + kns = kn total samples are drawn. We will ignore integer rounding issues, and assume knx = (1− ps)αn, and kns = psαn. The following lemma shows that the subsequence defined by kn converges almost surely. Lemma 5. Let α > 1, and kn = α n. Under the assumptions of Theorem 1 and that Lij ≥ 0
P   1
knxk2ns
knx∑
i=1
  kns∑
j=1
L̄ij − E[L̄ij | Zi]
  2
> ǫ i.o.   = 0.
We now must show that the sequence in (6) is a Cauchy sequence. Note that the SLLN implies that
1
nx
nx∑
i=1
E[L̄ij | Zi]2 a.s.→ E[E[L̄ij | Zi]2],
and the LLN for exchangeable arrays (Rieders, 1991, Theorem 1.4) implies that
1
nx
nx∑
i=1
1
ns
ns∑
j=1
L̄ijE[L̄ij | Zi] a.s.→ E[E[L̄ij | Zi]2].
Therefore,
1
knxk2ns
knx∑
i=1
  kns∑
j=1
L̄ij
  2
a.s.→ E[E[L̄ij | Zi]2]. (7)
Notice that because L̄ij ≥ 0, the sum ∑nx
i=1 (∑ns j=1 L̄ij )2 is monotone increasing in ns and nx.
With this in mind, for any m > 0, let n be such that kn ≤ m < kn+1. Then, by the montonicity, (
kn kn+1 1 kn
)3 knx∑
i=1
  kns∑
j=1
L̄ij
  2 ≤ ∑(1−ps)m i=1 (∑psm j=1 L̄ij )2
p2s(1− ps)m3 ≤ ( kn+1 kn 1 kn+1 )3 k(n+1)x∑
i=1
  k(n+1)s∑
j=1
L̄ij
  2
.
From (7), the left hand side converges to 1α3E[E[L̄ij | Zi]2], and the right hand side converges to α3E[E[L̄ij | Zi]2]. Because α is arbitrary, this proves that the sequence

∑(1−ps)m i=1 (∑psm j=1 L̄ij )2
p2s(1− ps)m3
 
m=1,...
is almost surely Cauchy. Together with Lemma 5, this implies (6).

heading: B.5 PROOF OF LEMMA 5
text: We will show that
∞∑
n=1
P   1
knxk2ns
knx∑
i=1
  kns∑
j=1
L̄ij − E[L̄ij | Zi]
  2
> ǫ   < ∞.
This, along with the first Borel-Cantelli lemma (Émile Borel, 1909; Cantelli, 1917) implies the result.
Applying Markov’s inequality and using the fact that L̄ij and L̄ih are independent conditional on Zi gives
∞∑
n=1
P   1
knxk2ns
knx∑
i=1
  kns∑
j=1
L̄ij − E[L̄ij | Zi]
  2
> ǫ  
≤ 1 ǫ
∞∑
n=1
E   1
knxk2ns
knx∑
i=1
  kns∑
j=1
L̄ij − E[L̄ij | Zi]
  2  
= 1
ǫ
∞∑
n=1
1
knxk2ns
knx∑
i=1
kns∑
j=1
E [( L̄ij − E[L̄ij | Zi] )2]
≤ 1 ǫ
∞∑
n=1
1
knxk2ns
knx∑
i=1
kns∑
j=1
E[L̄2ij ],
where the last line follows from the law of total variance. To simplify the remaining algebra, we will use a . b to denote that there is some constant 0 < c < ∞ such that a < cb. Continuing, we have
1
ǫ
∞∑
n=1
1
knxk2ns
knx∑
i=1
kns∑
j=1
E[L̄2ij ] . 1
ǫ
∞∑
n=1
knx∑
i=1
kns∑
j=1
1
k3n E[L̄2ij ]
= 1
ǫ
∞∑
i=1
∞∑
j=1
E[L̄2ij ]
∞∑
n=n(i,j)
1
α3n
. 1
ǫ
∞∑
i=1
∞∑
j=1
1
max{i/(1− ps), j/ps}3 E[L̄2ij ]
. 1
ǫ
∞∑
i=1
∞∑
j=1
1
max{i, j}3E[L̄ 2 ij ]
= 1
ǫ
∞∑
i=1
∞∑
j=1
1
max{i, j}3E[L̄ 2 ij ]
where n(i, j) is shorthand for n(i, j) = logα max{i/(1− ps), j/ps} is the first n such that knx ≥ i and kns ≥ j. Now, define Q as the distribution of L11 induced by Z1 and S1. Additionally, split the inner sum into two pieces, one for when j < i and so max{i, j} = i and one for when j ≥ i and so max{i, j} = j.
1
ǫ
∞∑
i=1
∞∑
j=1
1
max{i, j}3E[L̄ 2 ij ] =
1
ǫ
∞∑
i=1
  i∑
j=1
1
i3
∫ ij
0
x2 dQ(x) +
∞∑
j=i
∫ ij
0
x2 dQ(x)
 
= 1
ǫ
∞∑
i=1
  i−1∑
j=1
1
i3
ij∑
k=1
∫ k
k−1 x2 dQ(x) +
∞∑
j=i
ij∑
k=1
∫ k
k−1 x2 dQ(x)
 
switching the order of the indices over j and k, using that 1 ≤ k ≤ ij and the constraints on j relative to i,
1
ǫ
∞∑
i=1
  i−1∑
j=1
1
i3
ij∑
k=1
∫ k
k−1 x2 dQ(x) +
∞∑
j=i
ij∑
k=1
∫ k
k−1 x2 dQ(x)
 
. 1
ǫ
∞∑
i=1
  i2−1∑
k=1
(i− k/i) i3
∫ k
k−1 x2 dQ(x) +
∞∑
k=1
∞∑
j=max{i,k/i}
1
j3
∫ k
k−1 x2 dQ(x)
 
. 1
ǫ
∞∑
i=1
  i2−1∑
k=1
(i− k/i) i3
∫ k
k−1 x2 dQ(x) +
∞∑
k=1
1 max{i, k/i}2 ∫ k k−1 x2 dQ(x)
  .
Switching the order of summation over i and k, and separating out the terms where k/i < i and k/i ≥ i,
1
ǫ
∞∑
i=1
  i2−1∑
k=1
(i− k/i) i3
∫ k
k−1 x2 dQ(x) +
∞∑
k=1
1 max{i, k/i}2 ∫ k k−1 x2 dQ(x)
 
= 1
ǫ
∞∑
k=1
(∫ k
k−1 x2 dQ(x)
)  √ k+1∑
i=1
(i− k/i) i3 +
∞∑
i= √ k
1 i2 +
√ k∑
i=1
i2
k2
 
. 1
ǫ
∞∑
k=1
1√ k
(∫ k
k−1 x2 dQ(x)
)
. 1
ǫ
∞∑
k=1
(∫ k
k−1
x2√ x dQ(x)
)
. 1
ǫ
∫ ∞
0
x1.5 dQ(x) < ∞.
C INSTANCE-LEVEL AGREEMENT OF MULTIBERTS ON GLUE
We present additional performance experiments to complement Section 2.
Table 3 shows per-example agreement rates on GLUE predictions between pairs of models pretrained with a single seed (“same”) and pairs pre-trained with different seeds (“diff”); in all cases, models are fine-tuned with different seeds. With the exception of RTE, we see high agreement (over 90%) on test examples drawn from the same distribution as the training data, and note that agreement is 1–2% lower on average for the predictions of models pre-trained on different seeds compared to models pre-trained on the same seed. However, this discrepancy becomes significantly more pronounced if we look at out-of-domain “challenge sets” which feature a different data distribution from the training set. For example, if we evaluate our MNLI models on the anti-sterotypical examples from HANS (McCoy et al., 2019), we see agreement drop from 88% to 82% when comparing across pre-training seeds. Figure 4 shows how this can affect overall accuracy, which can vary over a range of nearly 20% depending on the pre-training seed. Such results underscore the need to evaluate multiple pre-training runs, especially when evaluating a model’s ability to generalize outside of its training distribution.

heading: D CROSS-SEED VARIATION
text: Figure 5 shows variation in Winogender bias correlation (S4) between each MultiBERTs pretraining seed. Each box shows the distribution over five runs, and some of the variation between seeds may simple be due to variation in training the coreference model. If we average the scores for each seed then look at the distribution of this per-seed average score, we get 0.45±0.11. What if pretraining didn’t matter? If we ignore the seed and randomly sample sets of five runs from this set with replacement, we get scores of 0.45±0.05 - telling us that most of the variance can only be explained by differences between the pretraining checkpoints.
We can confirm this by taking a subset of our pretraining seeds and training additional 25 randomlyinitialized coreference models. Figure 6 shows the result: seeds 0, 2, 3, and 4 appear closer together than in Figure 5, but seed 1 clearly has different properties with respect to our Winogender metric. We can confirm this with an unpaired multibootstrap analysis, taking seed 0 as base and seed 1 as experiment: we observe a significant effect of δ = 0.203 (p = 0.009), as shown in Table 4.

heading: E CASE STUDY: MULTIBERTS VS. ORIGINAL BERT
text: We now discuss challenges in reproducing the performance of the original BERT checkpoint, using the Multi-Bootstrap procedure.
The original bert-base-uncased checkpoint appears to be an outlier when viewed against the distribution of scores obtained using the MultiBERTs reproductions. Specifically, in reproducing the training recipe of Devlin et al. (2019), we found it difficult to simultaneously match performance on all tasks using a single set of hyperparameters. Devlin et al. (2019) reports training for 1M steps. However, as shown in Figure 1 and 2, models pre-trained for 1M steps matched the original checkpoint on SQuAD but lagged behind on GLUE tasks; if pre-training continues to 2M steps, GLUE performance matches the original checkpoint but SQuAD performance is significantly higher.
The above observations suggest two separate but related hypotheses (below) about the BERT pretraining procedure.
1. On most tasks, running BERT pre-training for 2M steps produces better models than 1M steps.
2. The MultiBERTs training procedure outperforms the original BERT procedure on SQuAD.
Let us use the Multi-Bootstrap to test these hypotheses.

heading: E.1 HOW MANY STEPS TO PRETRAIN?
text: Let f be the predictor induced by the BERT pre-training procedure using the default 1M steps, and let f ′ be the predictor resulting from the proposed intervention of training to 2M steps. From a glance at the histograms in Figure 8, we can see that MNLI appears to be a case where 2M is generally better, while MRPC and RTE appear less conclusive. With the MultiBERTs, we can test the significance of the results. The results are shown in Table 5. We find that MNLI conclusively performs better (δ = 0.007 with p < 0.001) with 2M steps; for RTE and MRPC we cannot reject the null hypothesis of no difference (p = 0.14 and p = 0.56 respectively).
As an example of the utility of this procedure, Figure 7 shows the distribution of individual samples of L for the intervention f ′ and baseline f from this bootstrap procedure (which we denote as L′ and L, respectively). The distributions overlap significantly, but the samples are highly correlated due to the paired sampling, and we find that individual samples of the difference (L′−L) are nearly always positive.

heading: E.2 DOES THE MULTIBERTS PROCEDURE OUTPERFORM ORIGINAL BERT ON SQUAD?
text: To test our second hypothesis, i.e., that the MultiBERTs procedure outperforms original BERT on SQuAD, we must use the unpaired Multi-Bootstrap procedure. In particular, we are limited to the case in which we only have a point estimate of L′(S), because we only have a single estimate of the performance of our baseline model f ′ (the original BERT checkpoint). However, the MultiBootstrap procedure still allows us to estimate variance across our MultiBERTs seeds and across the examples in the evaluation set. On SQuAD 2.0, we find that the MultiBERTs models trained for 2M
steps outperform original BERT with a 95% confidence range of 1.9% to 2.9% and p < 0.001 for the null hypothesis, corroborating our intuition from Figure 2.

Review summaries:
Many tasks in contemporary NLP begin by building off of a large language model.  This can cast the downstream task as a sort of fine-tuning experiment, whereby the results are heavily influenced by conditioning on the starting point of a single pre-trained version of an LLM.  In this work, the authors take BERT as an example, and ask how much a specific artifact as a draw from the distribution over (model weights, initialization scheme, training data, loss function) affects downstream tasks built upon it.    The authors provide a wide variety of BERT models that are varied in their training and initalization.  They define a bootstrap procedure for the scenario where multiple instantiations of base models are available, and tie their findings together in a case study of gender bias in coreference resolution.