Accounting for Variance in Machine Learning Benchmarks

Introduction: trustworthy benchmarks account for fluctuations

Machine learning increasingly relies upon empirical evidence to validate publications or efficacy. The value of a new method or algorithm $A$ is often established by empirical benchmarks comparing it to prior work. Although such benchmarks are built on quantitative measures of performance, uncontrolled factors can impact these measures and dominate the meaningful difference between the methods. In particular, recent studies have shown that loose choices of hyper-parameters lead to non-reproducible benchmarks and unfair comparisons Raff (2019; 2021); Lucic et al. (2018); Henderson et al. (2018); Kadlec et al. (2017); Melis et al. (2018); Bouthillier et al. (2019); Reimers & Gurevych (2017); Gorman & Bedrick (2019). Properly accounting for these factors may go as far as changing the conclusions for the comparison, as shown for recommender systems Dacrema et al. (2019), neural architecture pruning Blalock et al. (2020), and metric learning Musgrave et al. (2020).

The steady increase in complexity –e.g. neural-network depth– and number of hyper-parameters of learning pipelines increases computational costs of models, making brute-force approaches prohibitive. Indeed, robust conclusions on comparative performance of models $A$ and $B$ would require multiple training of the full learning pipelines, including hyper-parameter optimization and random seeding. Unfortunately, since the computational budget of most researchers can afford only a small number of model fits Bouthillier & Varoquaux (2020), many sources of variances are not probed via repeated experiments. Rather, sampling several model initializations is often considered to give enough evidence. As we will show, there are other, larger, sources of uncontrolled variation and the risk is that conclusions are driven by differences due to arbitrary factors, such as data order, rather than model improvements.

The seminal work of Dietterich (1998) studied statistical tests for comparison of supervised classification learning algorithms focusing on variance due to data sampling. Following works Nadeau & Bengio (2000); Bouckaert & Frank (2004) perpetuated this focus, including a series of work in NLP Riezler & Maxwell (2005); Taylor Berg-Kirkpatrick & Klein (2012); Anders Sogaard & Alonso (2014) which ignored variance extrinsic to data sampling. Most of these works recommended the use of paired tests to mitigate the issue of extrinsic sources of variation, but Hothorn et al. (2005) then proposed a theoretical framework encompassing all sources of variation. This framework addressed the issue of extrinsic sources of variation by marginalizing all of them, including the hyper-parameter optimization process. These prior works need to be confronted to the current practice in machine learning, in particular deep learning, where 1) the machine-learning pipelines has a large number of hyper-parameters, including to define the architecture, set by uncontrolled procedures, sometimes manually, 2) the cost of fitting a model is so high that train/validation/test splits are used instead of cross-validation, or nested cross-validation that encompasses hyper-parameter optimization (Bouthillier & Varoquaux, 2020).

In Section 2, we study the different source of variation of a benchmark, to outline which factors contribute markedly to uncontrolled fluctuations in the measured performance. Section 3 discusses estimation the performance of a pipeline and its uncontrolled variations with a limited budget. In particular we discuss this estimation when hyper-parameter optimization is run only once. Recent studies emphasized that model comparisons with uncontrolled hyper-parameter optimization is a burning issue Lucic et al. (2018); Henderson et al. (2018); Kadlec et al. (2017); Melis et al. (2018); Bouthillier et al. (2019); here we frame it in a statistical context, with explicit bias and variance to measure the loss of reliability that it incurs. In Section 4, we discuss criterion using these estimates to conclude on whether to accept algorithm $A$ as a meaningful improvement over algorithm $B$, and the error rates that they incur in the face of noise.

Based on our results, we issue in Section 5 the following recommendations:

As many sources of variation as possible should be randomized whenever possible. These include weight initialization, data sampling, random data augmentation and the whole hyperparameter optimization. This helps decreasing the standard error of the average performance estimation, enhancing precision of benchmarks.

Deciding of whether the benchmarks give evidence that one algorithm outperforms another should not build solely on comparing average performance but account for variance. We propose a simple decision criterion based on requiring a high-enough probability that in one run an algorithm outperforms another.

Resampling techniques such as out-of-bootstrap should be favored instead of fixed held-out test sets to improve capacity of detecting small improvements.

Before concluding, we outline a few additional considerations for benchmarking in Section 6.

The variance in ML benchmarks

Machine-learning benchmarks run a complete learning pipeline on a finite dataset to estimate its performance. This performance value should be considered the realization of a random variable. Indeed the dataset is itself a random sample from the full data distribution. In addition, a typical learning pipeline has additional sources of uncontrolled fluctuations, as we will highlight below. A proper evaluation and comparison between pipelines should thus account for the distributions of such metrics.

Here we extend the formalism of Hothorn et al. (2005) to model the different sources of variation in a machine-learning pipeline and that impact performance measures. In particular, we go beyond prior works by accounting for the choice of hyperparameters in a probabilistic model of the whole experimental benchmark. Indeed, choosing good hyperparameters –including details of a neural architecture– is crucial to the performance of a pipeline. Yet these hyperparameters come with uncontrolled noise, whether they are set manually or with an automated procedure.

The training procedure builds a predictor given a training set $S^{t}$. But since it requires specifying hyperparameters $\lambda$, a complete learning pipeline has to tune all of these. A complete pipeline will involve a hyper-parameter optimization procedure, which will strive to find a value of $\lambda$ that minimizes objective

After hyperparameters have been tuned, it is often customary to retrain the predictor using the full data $S^{tv}$. The complete learning pipeline $\mathcal{P}$ will finally return a single predictor:

The full learning procedure $\mathcal{P}$ described above yields a model $h^{*}$ $\textstyle\kern-0.2pt\bigwedge\kern-0.2pt$ . We now must define a metric that we can use to evaluate the performance of this model with statistical tests. For simplicity, we will use the same evaluation metric $e$ on which we based hyperparameter optimization. The expected risk obtained by applying the full learning pipeline $\mathcal{P}$ to datasets $S^{tv}\sim\mathcal{D}^{n}$ of size $n$ is:

where the expectation is also over the random sources $\xi$ that affect the learning procedure (initialization, ordering, data-augmentation) and hyperparameter optimization.

As we only have access to a single finite dataset $S$, the performance of the learning pipeline can be evaluated as the following expectation over splits:

2 Empirical evaluation of variance in benchmarks

We conducted thorough experiments to probe the different sources of variance in machine learning benchmarks.

We selected i) the CIFAR10 (Krizhevsky et al., 2009) image classification with VGG11 (Simonyan & Zisserman, 2014), ii) PascalVOC (Everingham et al., ) image segmentation using an FCN Long et al. (2014) with a ResNet18 He et al. (2015a) backbone pretrained on imagenet Deng et al. (2009), iii-iv) Glue Wang et al. (2019) SST-2 Socher et al. (2013) and RTE Bentivogli et al. (2009) tasks with BERT Devlin et al. (2018) and v) peptide to major histocompatibility class I (MHC I) binding predictions with a shallow MLP. All details on default hyperparameters used and the computational environments –which used $\sim 8$ GPU years– can be found in Appendix D.

For the sources of variance from the learning procedure ($\xi_{O}$), we identified: i) the data sampling, ii) data augmentation procedures, iii) model initialization, iv) dropout, and v) data visit order in stochastic gradient descent. We model the data-sampling variance as resulting from training the model on a finite dataset $S$ of size $n$, sampled from an unknown true distribution. $S\sim\mathcal{D}^{n}$ is thus a random variable, the standard source of variance considered in statistical learning. Since we have a single finite dataset in practice, we evaluate this variance by repeatedly generating a train set from bootstrap replicates of the data and measuring the out-of-bootstrap error Hothorn et al. (2005)The more common alternative in machine learning is to use cross-validation, but the latter is less amenable to various sample sizes. Bootstrapping is discussed in more detail in Appendix B..

We first fixed hyperparameters to pre-selected reasonable choicesThis choice is detailed in Appendix D.. Then, iteratively for each sources of variance, we randomized the seeds 200 times, while keeping all other sources fixed to initial values. Moreover, we measured the numerical noise with 200 training runs with all fixed seeds.

Figure 1 presents the individual variances due to sources from within the learning algorithms. Bootstrapping data stands out as the most important source of variance. In contrast, model initialization generally is less than 50% of the variance of bootstrap, on par with the visit order of stochastic gradient descent. Note that these different contributions to the variance are not independent, the total variance cannot be obtained by simply adding them up.

For classification, a simple binomial can be used to model the sampling noise in the measure of the prediction accuracy of a trained pipeline on the test set. Indeed, if the pipeline has a chance $\tau$ of giving the wrong answer on a sample, makes i.i.d. errors, and is measured on $n$ samples, the observed measure follows a binomial distribution of location parameter $\tau$ with $n$ degrees of freedom. If errors are correlated, not i.i.d., the degrees of freedom are smaller and the distribution is wider. Figure 2 compares standard deviations of the performance measure given by this simple binomial model to those observed when bootstrapping the data on the three classification case studies. The match between the model and the empirical results suggest that the variance due to data sampling is well explained by the limited statistical power in the test set to estimate the true performance.

To study the $\xi_{H}$ sources of variation, we chose three of the most popular hyperparameter optimization methods: i) random search, ii) grid search, and iii) Bayesian optimization. While grid-search in itself has no random parameters, the specific choice of the parameter range is arbitrary and can be an uncontrolled source of variance (e.g., does the grid size step by powers of 2, 10, or increments of 0.25 or 0.5). We study this variance with a noisy grid search, perturbing slightly the parameter ranges (details in Appendix E).

For each of these tuning methods, we held all $\xi_{O}$ fixed to random values and executed 20 independent hyperparameter optimization procedures up to a budget of 200 trials. This way, all the observed variance across the hyperparameter optimization procedures is strictly due to $\xi_{H}$. We were careful to design the search space so that it covers the optimal hyperparameter values (as stated in original studies) while being large enough to cover suboptimal values as well.

Results in figure 1 show that hyperparameter choice induces a sizable amount of variance, not negligible in comparison to the other factors. The full optimization curves of the 320 HPO procedures are presented in Appendix F. The three hyperparameter optimization methods induce on average as much variance as the commonly studied weights initialization. These results motivate further investigation the cost of ignoring the variance due to hyperparameter optimization.

For a given case study, the total variance due to arbitrary choices and sampling noise revealed by our study can be put in perspective with the published improvements in the state-of-the-art. Figure 3 shows that this variance is on the order of magnitude of the individual increments. In other words, the variance is not small compared to the differences between pipelines. It must be accounted for when benchmarking pipelines.

But before delving into this, we will explain why many splits help estimating the expected empirical risk ($\hat{R}_{\mathcal{P}}$).

The majority of machine-learning benchmarks are built with fixed training and test sets. The rationale behind this design, is that learning algorithms should be compared on the same grounds, thus on the same sets of examples for training and testing. While the rationale is valid, it disregards the fact that the fundamental ground of comparison is the true distribution from which the sets were sampled. This finite set is used to compute the expected empirical risk ($\hat{R}_{\mathcal{P}}$ Eq 2.1), failing to compute the expected risk ($R_{\mathcal{P}}$ Eq 4) on the whole distribution. This empirical risk is therefore a noisy measure, it has some uncertainty because the risk on a particular test set gives limited information on what would be the risk on new data. This uncertainty due to data sampling is not small compared to typical improvements or other sources of variation, as revealed by our study in the previous section. In particular, figure 2 suggests that the size of the test set can be a limiting factor.

When comparing two learning algorithms $A$ and $B$, we estimate their expected empirical risks $\hat{R}_{\mathcal{P}}$ with $\hat{\mu}_{(k)}$, a noisy measure. The uncertainty of this measure is represented by the standard error $\frac{\sigma}{\sqrt{k}}$ under the normal assumptionOur extensive numerical experiments show that a normal distribution is well suited for the fluctuations of the risk–figure G.3 of $\hat{R}_{e}$. This uncertainty is an important aspect of the comparison, for instance it appears in statistical tests used to draw a conclusion in the face of a noisy evidence. For instance, a z-test states that a difference of expected empirical risk between $A$ and $B$ of at least $z_{0.05}\sqrt{\frac{\sigma_{A}^{2}+\sigma_{B}^{2}}{k}}$ must be observed to control false detections at a rate of 95%. In other words, a difference smaller than this value could be due to noise alone, e.g. different sets of random splits may lead to different conclusions.

2 Bias and variance of estimators depends on whether they account for all sources of variation

Note that $T$ does not appear in these equations. Yet it controls $\operatorname{HOpt}$’s runtime cost ($T$ trials to determine $\hat{\lambda^{*}}$), and thus the variance $\sigma^{2}$ is a function of $T$.

2.2 Biased estimator: fixing HOpt

When we fix sources of variation $\xi$ to arbitrary values (e.g. random seed), we are conditioning the distribution of $\hat{R}_{e}$ on some arbitrary $\xi$. Intuitively, holding fix some sources of variations should reduce the variance of the whole process. What our intuition fails to grasp however, is that this conditioning to arbitrary $\xi$ induces a correlation between the trainings which in turns increases the variance of the estimator. Indeed, a sum of correlated variables increases with the strength of the correlations.

We present results from a subset of the tasks in Figure 5 (all tasks are presented in Figure H.4). Randomizing weights initialization only (FixedHOptEst($k$,init)) provides only a small improvement with $k>1$. In the task where it best performs (Glue-RTE), it converges to the equivalent of $\hat{\mu}_{(k=2)}$. This is an important result since it corresponds to the predominant approach used in the literature today. Bootstrapping with FixedHOptEst($k$,Data) improves the standard error for all tasks, converging to equivalent of $\hat{\mu}_{(k=2)}$ to $\hat{\mu}_{(k=10)}$. Still, the biased estimator including all sources of variations excluding hyperparameter optimization FixedHOptEst($k$,All) is by far the best estimator after the ideal estimator, converging to equivalent of $\hat{\mu}_{(k=2)}$ to $\hat{\mu}_{(k=100)}$.

This shows that accounting for all sources of variation reduces the likelihood of error in a computationally achievable manner. IdealEst($k=100$) takes 1 070 hours to compute, compared to only 21 hours for each FixedHOptEst($k=100$). Our study paid the high computational cost of multiple rounds of FixedHOptEst($k$,All), and the cost of IdealEst($k$) for a total of 6.4 GPU years to show that FixedHOptEst($k$,All) is better than the status-quo and a satisfying option for statistical model comparisons without these prohibitive costs.

Accounting for variance to draw reliable conclusions

Given an estimate of the performance of two learning pipelines and their variance, are these two pipelines different in a meaningful way? We first formalize common practices to draw such conclusions, then characterize their error rates.

A typical criterion to conclude that one algorithm is superior to another is that one reaches a performance superior to another by some (often implicit) threshold $\delta$. The choice of the threshold $\delta$ can be arbitrary, but a reasonable one is to consider previous accepted improvements, e.g. improvements in Figure 3.

This difference in performance is sometimes computed across a single run of the two pipelines, but a better practice used in the deep-learning community is to average multiple seeds Bouthillier & Varoquaux (2020). Typically hyperparameter optimization is performed for each learning algorithm and then several weights initializations or other sources of fluctuation are sampled, giving $k$ estimates of the risk $\hat{R}_{e}$ – note that these are biased as detailed in subsubsection 3.2.2. If an algorithm $A$ performs better than an algorithm $B$ by at least $\delta$ on average, it is considered as a better algorithm than $B$ for the task at hand. This approach does not account for false detections and thus can not easily distinguish between true impact and random chance.

Let $\hat{R}_{e}^{A}=\frac{1}{k}\sum_{i=1}^{k}\hat{R}_{ei}^{A}$, where $\hat{R}_{ei}^{A}$ is the empirical risk of algorithm $A$ on the $i$-th split, be the mean performance of algorithm $A$, and similarly for $B$. The decision whether $A$ outperforms $B$ is then determined by $({\hat{R}_{e}^{A}-\hat{R}_{e}^{B}>\delta})$.

The variance is not accounted for in the average comparison. We will now present a statistical test accounting for it. Both comparison methods will next be evaluated empirically using simulations based on our case studies.

2 Characterizing errors of these conclusion criteria

We now run an empirical study of the two conclusion criteria presented above, the popular comparison of average differences and our recommended probability of outperforming. We will re-use mean and variance estimates from subsection 3.3 with the ideal and biased estimators to simulate performances of trained algorithms so that we can measure the reliability of these conclusion criteria when using ideal or biased estimators.

For decisions based on comparing averages, we set $\delta=1.9952\sigma$ where $\sigma$ is the standard deviation measured in our case studies with the ideal estimator. The value 1.9952 is set by linear regression so that $\delta$ matches the average improvements obtained from paperswithcode.com. This provides a threshold $\delta$ representative of the published improvements. For the probability of outperforming, we use a threshold of $\gamma=0.75$ which we have observed to be robust across all case studies (See Appendix I).

Figure 6 reports results for different decision criteria, using the ideal estimator and the biased estimator, as the difference in performance of the algorithms $A$ and $B$ increases (x-axis). The x-axis is broken into three regions: 1) Leftmost is when $H_{0}$ is true (not-significant). 2) The grey middle when the result is significant, but not meaningful in our framework ($\cancel{H_{0}}\cancel{H_{1}}$). 3) The rightmost is when $H_{1}$ is true (significant and meaningful). The single point comparison leads to the worst decision by far. It suffers from both high false positives ($\approx 10\%$) and high false negatives ($\approx 75\%$). The average with $k=50$, on the other hand, is very conservative with low false positives ($<5\%$) but very high false negatives ($\approx 90\%$). Using the probability of outperforming leads to better balanced decisions, with a reasonable rate of false positives ($\approx 5\%$) on the left and a reasonable rate of false negatives on the right ($\approx 30\%$).

The main problem with the average comparison is the threshold. A t-test only differs from an average in that the threshold is computed based on the variance of the model performances and the sample size. It is this adjustment of the threshold based on the variance that allows better control on false negatives.

Our recommendations: good benchmarks with a budget

We now distill from the theoretical and empirical results of the previous sections a set of practical recommendations to benchmark machine-learning pipelines. Our recommendations are pragmatic in the sense that they are simple to implement and cater for limited computational budgets.

Fitting and evaluating a modern machine-learning pipeline comes with many arbitrary aspects, such as the choice of initializations or the data order. Benchmarking a pipeline given a specific instance of these choices will not give an evaluation that generalize to new data, even drawn from the same distribution. On the opposite, a benchmark that varies these arbitrary choices will not only evaluate the associated variance (section 2), but also reduce the error on the expected performance as they enable measures of performance on the test set that are less correlated (3). This counter-intuitive phenomenon is related to the variance reduction of bagging Breiman (1996a); Bühlmann et al. (2002), and helps characterizing better the expected behavior of a machine-learning pipeline, as opposed to a specific fit.

The subset of the data used as test set to validate an algorithm is arbitrary. As it is of a limited size, it comes with a limited estimation quality with regards to the performance of the algorithm on wider samples of the same data distribution (figure 2). Improvements smaller than this variance observed on a given test set will not generalize. Importantly, this variance is not negligible compared to typical published improvements or other sources of variance (figures 1 and 3). For pipeline comparisons with more statistical power, it is useful to draw multiple tests, for instance generating random splits with a out-of-bootstrap scheme (detailed in appendix B).

Additional considerations

There are many aspects of benchmarks which our study has not addressed. For completeness, we discuss them here.

Our framework provides value when the user can control the model training process and source of variation. In cases where models are given but not under our control (e.g., purchased via API or a competition), the only source of variation left is the data used to test the model. Our framework and analysis does not apply to such scenarios.

We focused on comparing two learning algorithms. Benchmarks – and competitions in particular – commonly involve large number of learning algorithms that are being compared. Part of our results carry over unchanged in such settings, in particular those related to variance and performance estimation. With regards to reaching a well-controlled decision, a new challenge comes from multiple comparisons when there are many algorithms. A possible alley would be to adjust the decision threshold $\gamma$, raising it with a correction for multiple comparisons (e.g. Bonferroni) Dudoit et al. (2003). However, as the number gets larger, the correction becomes stringent. In competitions where the number of contestants can reach hundreds, the choice of a winner comes necessarily with some arbitrariness: a different choice of test sets might have led to a slightly modified ranking.

Comparison over multiple datasets is often used to accumulate evidence that one algorithm outperforms another one. The challenge is to account for different errors, in particular different levels of variance, on each dataset.

Demšar (2006) recommended Wilcoxon signed ranks test or Friedman tests to compare classifiers across multiple datasets. These recommendations are however hardly applicable on small sets of datasets – machine learning works typically include as few as 3 to 5 datasets Bouthillier & Varoquaux (2020). The number of datasets corresponds to the sample size of these tests, and such a small sample size leads to tests of very limited statistical power.

Dror et al. (2017) propose to accept methods that give improvements on all datasets, controlling for multiple comparisons. As opposed to Demšar (2006)’s recommendation, this approach performs well with a small number of datasets. On the other hand, a large number of datasets will increase significantly the severity of the family-wise error-rate correction, making Demšar’s recommendations more favorable.

We focused on model performance, but model evaluation in practice can include other metrics such as the training time to reach a performance level or the memory foot-print Reddi et al. (2020). Performance metrics are generally averages over samples which typically makes them amenable to a reasonable normality assumption.

Conclusion

We showed that fluctuations in the performance measured by machine-learning benchmarks arise from many different sources. In deep learning, most evaluations focus on the effect of random weight initialization, which actually contribute a small part of the variance, on par with residual fluctuations of hyperparameter choices after their optimization but much smaller than the variance due to perturbing the split of the data in train and test sets. Our study clearly shows that these factors must be accounted to give reliable benchmarks. For this purpose, we study estimators of benchmark variance as well as decision criterion to conclude on an improvement. Our findings outline recommendations to improve reliability of machine learning benchmarks: 1) randomize as many sources of variations as possible in the performance estimation; 2) prefer multiple random splits to fixed test sets; 3) account for the resulting variance when concluding on the benefit of an algorithm over another.

References

Appendix A Notes on reproducibility

Ensuring full reproducibility is often a tedious work. We provide here notes and remarks on the issues we encountered while working towards fully reproducible experiments.

To ensure proper study of the sources of variation it was necessary to control them close to perfection. For all tasks, we ran a pipeline of tests to ensure perfect reproducibility at execution and also at resumption. During the tests, each source of variation was varied with 5 different seeds, each executed 5 times. This ensured that the pipeline was reproducible for different seeds. Additionally, for each source of variation and for each seed, another training was executed but automatically interrupted after each epoch. The worker would then start the training of the next seed and iterate through the trainings for all seeds before resuming the first one. All these tests uncovered many bugs and typical reproducibility issues in machine learning. We report here some notes.

Although we did not measure the variance induced by different GPU architectures, we did observe that different GPU models would lead to different results. The CPU model had less impact on the Deep Learning tasks but the MLP-MHC task was sensitive to it. We therefore limited all tasks to specific computer architectures. We also observed issues when CUDA drivers were updated during preliminary experiments. We ensured all experiments were run using CUDA 10.2.

PyTorch versions lead to different results as well. We ran every Deep Learning experiments with PyTorch 1.2.0.

We implemented our data pipeline so that we could seed the iterators, the data augmentation objects and the splitting of the datasets. We had less control at the level of the models however. For PyTorch 1.2.0, the random number generator (RNG) must be seeded globally which makes it difficult to seed different parts separately. We seeded PyTorch’s global RNG for weight initialization at the beginning of the training process and then seeded PyTorch’s RNG for the dropout. Afterwards we checkpoint the RNG state so that we can restore the RNG states at resumption. We found that models with convolutionnal layers would not yield reproducible results unless we enabled cudnn.deterministic and disabled cudnn.benchmark.

We used the library RoBO Klein et al. (2017) for our Bayesian Optimizer. There was no support for seeding, we therefore resorted to seeding the global seed of python and numpy random number generators. We needed again to keep track of the RNG states and checkpoint them so that we can resume the Bayesian Optimizer without harming the reproducibility.

For one of our case study, image segmentation, we have been unable to make the learning pipeline perfectly reproducible. This is problematic because it prevents us from studying each source of variation in isolation. We thus trained our model with every seeds fixed across all 200 trainings and measured the variance we could not control. This is represented as the numerical noise in Figures 1 and G.3.

Appendix B Our bootstrap procedure

Cross-validation with different $k$ impacts the number of samples, it is not the case with not bootstrap. That means flexible sample sizes for statistical tests is hardly possible with cross-validation within affecting the training dataset sizes. Hothorn et al. (2005) focuses on the dataset sampling as the most important source of variation and marginalize out all other sources by taking the average performance over multiple runs for a given dataset. This increases even more the computational cost of the statistical tests.

We probe the effect of data sampling with bootstrap, specifically by bootstrapping to generate training sets and measuring the out-of-bootstrap error, as introduced by Breiman (1996b) in the context of bagging and generalized by Hothorn et al. (2005). For completeness, we formalize this use of the bootstrap to create training and test sets and how it can estimate the variance of performance measure due to data sampling on a finite dataset.

We assume we are seeking to generate sets of i.i.d. samples from true distribution $\mathcal{D}$. Ideally we would have access to $\mathcal{D}$ and could sample our finite datasets independently from it.

Instead we have one dataset $S\sim\mathcal{D}^{n}$ of finite size $n$ and need to sample independent datasets from it. A popular method in machine learning to estimate performance on a small dataset is cross-validation Bouckaert & Frank (2004); Dietterich (1998). This method however underestimates variance because of correlations induced by the process. We instead favor bootstrapping Efron (1979) as used by Hothorn et al. (2005) to simulate independent data sampling from the true distribution.

Where $S^{t}_{b}\sim S$ represents sampling the $b$-th training set with replacement from the set $S$. We then turn to out-of-bootstrapping to generate the held-out set. We use all remaining samples in $S\setminus S^{t}_{b}$ to sample $S^{o}_{b}$.

This procedure is represented as $(S^{tv},S^{o})\sim sp_{n,n^{\prime}}(S)$ in the empirical average risk $\hat{R}_{\mathcal{P}}(S,n,n^{\prime})$, end of Section 2.1.

Appendix C Statistical testing

We are interested in asserting whether a learning algorithm $A$ better performs than another learning algorithm $B$. Measuring the performance of these learning algorithms is not a deterministic process however and we may be deceived if noise is not accounted for. Because of the noise, we cannot know for sure whether a conclusion we draw is true, but using a statistical test, we can at least ensure a bounded rate of false positives (drawing $A>B$ while truth is $A\leq B$) and false negatives (drawing $A\leq B$ while truth is $A>B$). The capacity of a statistical test to identify true differences, that is, of correctly inferring $A>B$ when this is true, is called the statistical power of a test. The procedure we describe here seeks to avoid deception from false positives while providing a strong statistical power.

As shown in Section 3, randomizing as many sources of variance as possible in the learning pipelines help reduce the correlation and thus improve the reliability of the performance estimation. The simplest way to randomize as many as possible is to simply avoid seeding the random number generators. We list here sources of variations we faced in our case studies, but there exists many other sources of variations in diverse learning algorithms and tasks.

The data being used should ideally always be different samples from the true distribution of interest. In practice we only have access to a finite dataset and therefore the best we can do is random splits with cross-validation or out-of-bootstrap as described in Appendix B.

The ordering of the data can have a surprisingly important impact as can be observed in Figure 1.

Stochastic data augmentation should not be seeded, so that it follows a different sequence at each run.

Model initialization, e.g. weights initialization in neural networks, should be randomized across all trainings.

Learning algorithms sometimes include stochastic computations such as dropout in neural networks Srivastava et al. (2014), or samplings methods Kingma & Welling (2014); Maddison et al. (2017).

The optimization of the hyperparameters generally include stochasticity which should ideally be randomized. Running multiple hyperparameter optimizations may often be practically unaffordable. Tests may still be carried out while fixing the hyperparameters after a single hyperparameter optimization, but keep in mind the incurred degradation of the reliability of the conclusion as shown in Section 4.

C.2 Pairing

Pairing is optional but is highly recommended to increase statistical power. Avoiding seeding is the simplest solution for the randomization, but it is not the best solution. If possible, meticulously seeding all sources of variation with different random seeds at each run makes it possible to pair trainings of the algorithms so that we can conduct paired comparisons.

Pairing is a simple but powerful way of increasing the power of statistical tests, that is, enabling the reliable detection of difference with smaller sample sizes. Let $\sigma_{A}$ and $\sigma_{B}$ be the standard deviation of the performance metric of learning algorithms $A$ and $B$ respectively. If measures of $\hat{R}_{e}^{A}$ and $\hat{R}_{e}^{A}$ are not paired, the standard deviation of $\hat{R}_{e}^{A}-\hat{R}_{e}^{B}$ is then $\sigma_{A}+\sigma_{B}$. If we pair them, then we marginalize out sources of variance which results in a smaller variance $\sigma_{A-B}\leq\sigma_{A}+\sigma_{B}$. This reduction of variance makes it possible to reliably detect smaller differences without increasing the sample size.

To pair the learning algorithms, sources of variation should be randomized similarly for all of them. For instance, the random split of the dataset obtained from out-of-bootstrap should be used for both $A$ and $B$ when making a comparison. Suppose we plan to execute 10 runs of $A$ and $B$, then we should generate 10 different splits $\{(S_{1}^{tv},S_{1}^{o}),(S_{2}^{tv},S_{2}^{o}),\cdots,(S_{10}^{tv},S_{10}^{o})\}$ and train $A$ and $B$ on each. The performances $(\hat{R}_{ei}^{A},\hat{R}_{ei}^{B})$ would then be compared only on the corresponding splits $(S_{i}^{tv},S_{i}^{o})$. The same would apply to all other sources of variations. In practical terms, pairing $A$ and $B$ requires sampling seeds for each pairs, re-using the same seed for $A$ and $B$ in each pairs.

For some sources of variation it may not make sense to pair. This is the case for instance with weights initialization if $A$ and $B$ involve different neural network architectures. We can still pair. This would not help much, but would not hurt as well. In doubt, it is better to pair.

C.3 Sample size

We must first set the threshold $\gamma$ for our test. Based on our experiments in Section 4, we recommend a value of 0.75. We then set the desired rates of false positives and false negatives with $\alpha$ and $\beta$ respectively. Usual value for $\alpha$ is 0.05 while $\beta$ ranges from 0.05 to 0.2. We recommend $\beta=0.05$ for a strong statistical power.

The estimation of $P(A>B)$ is equivalent to a Mann–Whitney test Perme & Manevski (2019), thus we can use Noether’s sample size determination method for this type of test Noether (1987).

Where $\Phi^{-1}$ is the inverse cumulative function of the normal distribution.

Figure C.1 shows how the minimal sample size evolves with $\gamma$. Detecting $P(A>B)$ below $\gamma=0.6$ is unpractical, requiring more that 700 trainings below 0.55 for instance. For a threshold that is representative of the published improvements as presented in Figure 3, $\gamma=0.75$, the minimal sample size required to ensure a rate of 5% false negatives (as defined by $\beta=0.05$) is reasonably small; 29 trainings.

Let $\mbox{CI}_{\min}$ and $\mbox{CI}_{\max}$ be the lower and upper bounds of the confidence interval. We draw a conclusion based on the three following scenarios.

: Not statistically significant. No conclusion should be drawn as the result could be explained by noise alone.

: Statistically significant and meaningful. We can conclude that learning algorithm $A$ is better performing than $B$ in the conditions defined by the experiments.

Appendix D Case studies

CIFAR10 Krizhevsky et al. (2009) is a dataset of 60,000 32x32 color images selected from 80 million tiny images dataset Torralba et al. (2008), divided in 10 balanced classes. The original split contains 50,000 images for training and 10,000 images for testing. We applied random cropping and random horizontal flipping data augmentations.

The aggregation of all original training and testing samples are used for the bootstrap. To preserve the balance of the classes, we applied stratified bootstrap. For each class separately, we sampled with replacement 4,000 training samples, 1,000 for validation and 1,000 for testing. As for all tasks, we use out-of-bootstrap to ensure samples cannot be contained in more than one set.

We used VGG11 Simonyan & Zisserman (2014) with batch-normalization and no dropout. The weights are initialized with Glorot method based on a uniform distribution Glorot & Bengio (2010).

We focused on learning rate, weight decay, momentum and learning rate schedule. Batch-size was omitted to simplify the multi-model training on GPUs, so that memory usage was consistent and predictable across all hyperparameter settings. To ease the definition of the search space for the learning rate schedule, we used exponential decay instead of multi-step decay despite the wide use of the latter with similar tasks and models Simonyan & Zisserman (2014); Xie et al. (2019); Mahajan et al. (2018); Liu et al. (2018); He et al. (2016; 2015b). The former only require tuning of $\gamma$ while the later requires additionally selecting number of steps. Search space for all experiments and default values used for the variance experiments are presented in Table 2.

D.2 Glue-SST2 sentiment prediction with BERT

SST2 (Stanford Sentiment Treebank) Socher et al. (2013) is a binary classification task included in GLUE Wang et al. (2019). In this task, the input is a sentence from a collection of movie reviews, and the target is the associated sentiment (either positive or negative). The publicly available data contains around 68k entries.

We maintained the same size ratio between train/validation (i.e., 0.013) when performing the bootstrapping analysis. We performed standard out-of-bootstrap without conserving class balance since the original dataset is not balanced and ratios between classes vary from training and validation set in the original splits. The variable ratios of classes across bootstrap samples generate additional variance in our results, but is representative of the effect of generating a dataset that is not perfectly balanced.

We used the BERT Devlin et al. (2018) implementation provided by the Hugging Face Wolf et al. (2019) repository. BERT is a Transformer Vaswani et al. (2017) encoder pre-trained on the self-supervised Masked Language Model task Devlin et al. (2018). We chose BERT given its importance and influence in the NLP literature. It is worthy to note that the pre-training phase of BERT is also affected by sources of variations. Nevertheless, we didn’t investigate this phase given the amount of time (and resources) required to perform it. Instead, we always start from the (same) pre-trained model image provided by the Hugging Face Wolf et al. (2019) repository. Indeed, the weight initialization was only applied to the final classifier. The initialization method used is standard Gaussian with $0.0$ mean and standard deviation that depends on the related hyperparameter.

We ran a small-scale hyperparameter space exploration in order to select the hyperparameter search space to use in our experiments. As such, we decided to include the learning rate, weight decay and the standard deviation for the model parameter initialization (see Table 3). We fixed the dropout probability to the value of 0.1 as in the original BERT architecture. For the same reason, we fixed $\beta_{1}=0.9$ and $\beta_{2}=0.999$. Default values used for the variance experiments are also reported in Table 3. The model has been fine-tuned on SST2 for 3 epochs, with a batch size of 32. Training has been performed with mixed precision. Note that for weight decay we used the default value from the Hugging Face repository (i.e., $0.0$) even if this is outside of the hyperparameter search space. We confirmed that this makes no difference by looking at the results of the small-scale hyperparameter space exploration.

D.3 Glue-RTE entailment prediction with BERT

RTE (Recognizing Textual Entailment) Bentivogli et al. (2009) is a also a binary classification task included in GLUE Wang et al. (2019). The task is a collection of text fragment pairs, and the target is to predict if the first text fragment entails the second one. RTE dataset only contains around 2.5k entries.

In our bootstrapping analysis we maintained the train/validation ratio of 0.1. As for Glue-SST2, we used standard out-of-bootstrap and did not preserve original class ratios.

We used the BERT Devlin et al. (2018) model for RTE as well, trained in the same way specified in the SST-2 section. In particular, we used the same hyperparameters (see Table 3), same batch size, and we trained in the same mixed-precision environment. The model has been fine-tuned on RTE for 3 epochs.

D.4 PascalVOC image segmentation with ResNet Backbone

The PascalVOC segmentation task Everingham et al. entails generating pixel-level segmentations to classify each pixel in an image as one of 20 classes or background. This publicly available dataset contains 2913 images and associated ground truth segmentation labels. The original splits contains 2184 images for training and 729 for validation. Images were normalized and zero-padded to a final size of 512x512.

We used a train/validation ratio of 0.25 for our bootstrap analysis, generating training sets of 2184 images, validation and test sets of 729 images each. Since multiple classes can appear in a single image, the original dataset was not balanced, we thus used standard out-of-bootstrap for our experiments.

We used an FCN-16s Long et al. (2014) with a ResNet18 backbone He et al. (2015a) pretrained on ImageNet Deng et al. (2009). After exploring several possible backbones, ResNet18 was selected since it could be trained relatively quickly. We use weighted cross entropy, with only predictions within the original image boundary contributing to the loss. The model is optimized using SGD with momentum.

The metric used is the mean Intersection over Union (mIoU) of the twenty classes and the background class. The complement of the mIoU, the mean Jaccard Distance, is the metric minimized in all HPO experiments.

Certain hyperparmeters, such as the number of kernals, or the total number of layers, are part of the definition of the ResNet18 architecture. As a result, we explored key optimization hyperparameters including: learning rate, momentum, and weight decay. The hyperparameter ranges selected, as well as the default hyperparameters used in the variance experiments, can be found in table 5 and in table LABEL:table:hps-pascal-voc-default-appendix, respectively. A batch size of 16 was used for all experiments.

D.5 Major histocompatibility class I-associated peptide binding prediction with shallow MLP

The MLP-MHC is a regression task with the goal of predicting the relative binding affinity for a given peptide and major histocompatibility complex class I (MHC) allele pair. The major histocompatibility complex (MHC) class I proteins are present on the surface of most nucleated cells in jawed vertebrates Pearson et al. (2016). These proteins bind short peptides that arise from the degradation of intra-cellular proteins Pearson et al. (2016). The complex of peptide-MHC molecule is used by immune cells to recognize healthy cells and eliminate cancerous or infected cells, a mechanism studied in the development of immunotherapy and vaccines O’Donnell et al. (2018). The peptide binding prediction task is therefore at the base of the search for good vaccine and immunotherapy targets O’Donnell et al. (2018); Jurtz et al. (2017).

The input data is the concatenated pairs of sequences: the MHC allele and the peptide sequence. For the MHC alleles, we restricted the sequences to the binding pocket of the peptide, as seen in Jurtz et al. (2017). The prediction target is a normalized binding affinity score, as described in Jurtz et al. (2017); O’Donnell et al. (2018).

While both MHCflurry and NetMHCpan4 models use a BLOSUM62 encoding Henikoff & Henikoff (1992) for the amino acids, in we chose to instead encode the amino acids as one-hot as described in Nielsen et al. (2007).

The NetMHCpan4 model is trained on a manually filtered dataset from the immune epitope database Vita et al. (2019); Jurtz et al. (2017) that has been split into five folds used for cross-validation, available on the author’s website Jurtz et al. (2017).

In contrast, the MHCflurry model is trained on a custom multi-source dataset (available from Mendeley data and the O’Donnell et al. (2018) publication cite) and validated/tested on two external datasets from Pearson et al. (2016) and an HPV peptide dataset available at the same website as above.

We have three different sets for training, validating and testing. We thus performed bootstrapping separately on each set for every training and evaluation.

The model is a shallow MLP with one hidden layer from sklearn. We used the default setting for the non-linearity relu and weight initialization strategy (Glorot & Bengio (2010)). The following table (Table 9) offers some comparison points between our model and the NetMHCpan4 Jurtz et al. (2017) and MHCflurry O’Donnell et al. (2018) models.

While the MHCflurry model (O’Donnell et al., 2018) train only no the peptide sequences and uses ensembling to perform its predictions, training multiple models for each MHC allele, the NetMHCpan4 model (Jurtz et al., 2017) uses the allele sequence as input and trains one single model.

We chose to retain the strategy proposed by the NetMHCpan4 model, where a single model is trained for all alleles Jurtz et al. (2017). As a reference, MHCflurry uses ensembling to perform predictions; indeed, the authors report that for each MHC allele, an ensemble of 8-16 are selected from the 320 that were trained O’Donnell et al. (2018).

For the hyperparameter search, we selected hidden layer sizes between 20 and 400 (Table 6), to engulph a range slightly larger than the ones described by both Jurtz et al. (2017); O’Donnell et al. (2018). The second hyperparameter that was explored was the L2 regularisation parameter, for which a log-uniform range between 0 and 1 was explored.

We would like to state the goal of the present study was not to establish new state of the art (SOTA) on the MHC-peptide binding prediction task. However, we still report that when comparing the performance of our model to those of NetMHCpan4 and MHCflurry we found the performance of our model comparable. Briefly, for the results in Table 9, we used the existing pre-trained NetMHCpan4 and MHCflurry tools to predict the binding affinity of both datasets: the previously described HPV external test data (HPV) from O’Donnell et al. (2018) and the cross-validation test datasets from Jurtz et al. (2017) (NetMHC-CVsplits).

We would like to point out that since the MHCflurry model was published later than the NetMHCpan4 one, there is a high chance that the dataset from the cross-validation splits (NetMHC-CVsplits) may be contained in the dataset used to train the existing MHCflurry tool. The proper way to compare performances would be to re-train the MHCflurry model on each fold and test susequently its performance; however, since our goal is not to reach new SOTA on this task, we leave this experiment to be performed at a later time.

This would result in a likely overestimation of the performance of MHCflurry on this dataset, which we noted with the $*$ sign in Table 9.

A more in-depth study is necessary to compare in a more through way this performance with respect to the differences in model design, dataset encoding and other factors.

Appendix E Hyperparameter optimization algorithms

Let $a_{i}$, $b_{i}$ and $n$ be the hyperparameters of the grid search, where $a_{i}$ and $b_{i}$ are vectors of minima and maxima for each dimension of the search space, and $n$ is the number of values per dimension. We define $\Delta_{i}$ as the interval between each value on dimension $i$. A point on the grid is defined by $p_{ij}=a_{i}+\Delta_{i}(j-1)$. Grid search is simply the evaluation of $r(\lambda)$ from Equation 2 on all possible combinations of values $p_{ij}$.

E.2 Noisy Grid Search

Grid search is a fully deterministic algorithm. Yet, it is highly sensitive to the design of the grid. To provide a variance estimate of similar choices of the grid and to be able to distinguish lucky grid, we consider a noisy version of grid search.

This provides us a variance estimate of grid search that we can compare against non-deterministic hyperparameter optimization algorithms.

E.3 Random Search

The search space of random search will be increased by $\pm\frac{\Delta_{i}}{2}$ as defined for the noisy grid search to ensure that they both cover the same search space. For all hyperparameters, the values are sampled from a uniform $p_{i}\sim U(a_{i}-\frac{\Delta_{i}}{2},b_{i}+\frac{\Delta_{i}}{2})$. For learning rate and weight decay, values are sampled uniformly in the logarithmic space.

Appendix F Hyperparameter optimization results

Figure F.2 presents the optimization curves of the hyperparameter optimization executions in Section 2.2.

Appendix G Normality of performance distributions in the case studies

Figure G.3 presents the Shapiro-Wilk test of normality on all our results on sources of variations.

Appendix H Randomizing more sources of variance increase the quality of the estimator

Figure 5 only presented the Glue-RTE and CIFAR10 tasks. We provide here a complete picture of the standard deviation of the different estimators in Figure H.4. We further present a decomposition of the mean-squared-error in Figure H.5 to help understand why accounting for more sources of variations improves the mean-squared-error of the biased estimators.

Appendix I Analysis of robustness of comparison methods

In addition to simulations described in Section 4.2, we executed experiments in which we varied the sample size and the threshold $\gamma$. To select the threshold of the average, we converted $\gamma$ into the equivalent performance difference $(\delta=\Phi^{-1}(\gamma)\sigma)$. Results are presented in Figure I.6