mathx”17 e_ s o ¿1—m E#1_#1#2*#4#3#4[#3]#4 mm#1#2 — #2 e_ s o ¿1—m
P#1_#1#2*#4#3#4[#3]#4 mm#1#2 — #2

# Measuring model variability using robust non-parametric testing††thanks:
This work was partially supported by the Statistical Inference Generates
kNowledge for Artificial Learners (SIGNAL) program at Pacific Northwest
National Laboratory (PNNL), as well as by the Mathematics for Artificial
Reasoning in Science (MARS) initiative via the Laboratory Directed Research
and Development (LDRD) Program at PNNL. PNNL is a multi-program national
laboratory operated for the U.S. Department of Energy (DOE) by Battelle
Memorial Institute under Contract No. DE-AC05-76RL0-1830.

Sinjini Banerjee Corresponding author<EMAIL_ADDRESS>Department of
Electrical and Computer Engineering, Rutgers University  Tim Marrinan  Pacific
Northwest National Lab  Reilly Cannon  Pacific Northwest National Lab  Tony
Chiang  Pacific Northwest National Lab  Anand D. Sarwate  Department of
Electrical and Computer Engineering, Rutgers University

###### Abstract

Training a deep neural network often involves stochastic optimization, meaning
each run will produce a different model. The seed used to initialize random
elements of the optimization procedure heavily influences the quality of a
trained model, but this may be obscure from many commonly reported summary
statistics, like accuracy. For example, two classification models with
identical validation accuracy may disagree on many predicted labels,
suggesting that the models have learned two fundamentally different functions.
Thus, optimizing over the space of random seeds can lead to higher-quality
models, much like optimizing over other hyper-parameters, such as learning
rate or regularization parameters, can lead to higher-quality models. However,
random seed is often not included in hyper-parameter optimization, perhaps
because the relationship between seed and model quality is hard to describe.
It is also unknown whether the potential marginal quality improvement of
training an additional model with a different random seed is worth the
increased computational cost. This work attempts to describe the relationship
between deep net models trained with different random seeds and the behavior
of the expected model. We adopt robust hypothesis testing to propose a novel
summary statistic for network similarity, referred to as the $\alpha$-trimming
level. Our contributions are two-fold. First, we use the $\alpha$-trimming
level to show that the empirical cumulative distribution function of an
ensemble model created from a collection of trained models with different
random seeds approximates the average of these functions as the number of
models in the collection grows large. This insight provides guidance for how
many random seeds should be sampled to ensure that an ensemble of these
trained models is a reliable representative. Second, we show that the
$\alpha$-trimming level is more expressive than different performance metrics
like validation accuracy, churn, or expected calibration error when taken
alone and may help with random seed selection in a more principled fashion. We
demonstrate the value of the proposed statistic in real experiments and
illustrate the advantage of fine-tuning over random seed with an experiment in
transfer learning.

##  1 Introduction

Deep learning models have achieved state-of-the-art performance on complex
tasks in healthcare, education, cyber-security, and other critical domains.
Training these models takes significant time, energy, and financial resources.
Since training algorithms use stochastic optimization for non-convex
objectives, models produced by different training runs, in general, converge
to different solutions. These trained models may correspond to different
functions, even with a similar objective value and validation/test accuracy.
Models are also often continuously retrained as new data arrives. This
necessitates algorithmic and architectural changes in state-of-the-art models
to improve their performance on new data. The run-to-run variability in
training models makes it difficult to conclude if a specific initialization or
hyper-parameter tuning made a meaningful difference in model performance or if
it just “got lucky” due to the (unavoidable) presence of randomness in the
optimization. Without this knowledge, comparing training configurations to
assess whether one is better becomes difficult.

Figure 1:  Illustration of two models with the same validation accuracy but
different decision boundaries.

Consider the following stylized model of a standard machine learning working
process: a practitioner trains $M$ models with different hyper-parameter
settings, chooses the model with the best validation performance (fine-
tuning), and evaluates it on test data. This could also involve fine-tuning a
model pre-trained on large datasets for a downstream task using a final task-
specific layer. In this work, we consider the random seed in the stochastic
optimization as a hyper-parameter itself. Gundersen et al. [13] provided a
taxonomy of sources of irreproducibility in machine learning practice and
identified different seeds controlling model initialization, batch shuffling
during mini-batch stochastic gradient descent (SGD), and data sampling as some
of the major sources of randomness present in the training procedure. The
importance of using random seeds as part of hyper-parameter tuning in deep
neural network (DNN) training has been previously highlighted in the
literature. Henderson et al. [14], Reimers et al. [30], Melis et al. [24],
Bouthillier et al. [7], and Picard et al. [28] all found the effect of random
seed on model performance to be statistically significant. Other works have
focused on understanding the random seed effect on specific sources of
randomness in training. For example, Fort et al. [11] found that
initialization has a larger effect than batch order in SGD on model
performance, while Bouthillier et al. [6] showed the opposite. Summers and
Dinneen [32], and Jordan [18] show that networks converge to vastly different
values even with the change of a single bit of network parameters during
initialization. While studying run-to-run variability in pre-trained BERT
models over different random seeds, Dodge et al. [8] noted a validation
performance gain of 7% over previously reported results, highlighting the
importance of fine-tuning over random seed. All the above-mentioned works have
assessed the impact randomness has on the _validation accuracy_ , or _churn_
[25] between models, which quantifies how two models differ in their
predictions on the same test point. These summary statistics only focus on the
_decisions_ made by predictive (classification) models and do not directly
assess differences in the functions learned by these models.

In this paper, we assess DNN training variability over random seeds in terms
of the _network outputs used to make the decision_. Figure 1 illustrates the
difference: the solid and dashed lines represent two decision boundaries
between red (circle) and blue (cross). Accuracy measures incorrect decisions,
and the churn is determined by the region between the two curves. We are
interested in the shaded regions, which show the thresholded network output to
make decisions. If we think of the training algorithm as generating a random
sample from a function space, we can use other tools to understand model
variability. As a first step, we can examine the distribution of the pre-
thresholded outputs (the _logit gap_) from functions learned by different runs
of a fixed DNN architecture. To that end, we adopt a non-parametric hypothesis
testing framework to measure the similarity of functions learned by DNN models
that vary in random seeds.

This paper is organized as follows,

  * •

In section 2, we introduce notations and empirical cumulative distribution
functions (eCDFs) of interest.

  * •

In section 3, we introduce the classical two-sample Kolmogorov–Smirnov (KS)
test to measure the closeness between two models, propose a reference
distribution used to formulate the null hypothesis and derive a threshold for
this test.

  * •

With the limitations of a classical KS test in mind, in section 4, we
formulate a robust version of the two-sample KS test, and propose a measure of
closeness/discrepancy referred to as the $\alpha$-trimming level between a
model and the reference distribution.

  * •

In section 5, we provide empirical evidence that the proposed reference
distribution is close to the eCDF of a deep ensemble (a technique used in the
machine learning literature to boost model performance and reduce
variability). We perform two case studies highlighting the usefulness of our
proposed measure of closeness for random seed selection.

##  2 Problem setup

Notation: For a positive integer $n$, let $[n]=\\{1,2,\ldots,n\\}$. Random
variables will be denoted in boldface, with realizations being non-bolded, so
$\bm{\theta}$ is a random variable and $\theta$ is a realization. We consider
training a DNN for predicting an output taking values in a set ${\mathcal{Z}}$
with input from a (feature vector) $x$ taking values in a space
${\mathcal{X}}$. A DNN with a given architecture is specified by a set of
parameters (e.g. the weights) $\theta$ taking values in parameter space,
$\varTheta$.

We will restrict attention to _binary classification_ using a DNN with
parameters $\theta$. Given a label space ${\mathcal{Y}}$, a training algorithm
takes a training set,
$\mathcal{D}_{\mathrm{train}}=\\{(x_{i},y_{i})\in{\mathcal{X}}\times{\mathcal{Y}}\colon
i\in[N_{\mathrm{train}}]\\}$, estimates the weights $\theta$ that
(approximately) minimize the empirical loss
$\hat{R}(\theta;\mathcal{D}_{\mathrm{train}})$ over the training data, and
assigns a label $y\in{\mathcal{Y}}$ using a _softmax_ operation. We model this
by assuming the DNN takes an input $x$ and computes functions
$m^{+}(x\mid\theta)$ and $m^{-}(x\mid\theta)$, with the predicted label
$\hat{y}(x;\theta)$ being $+1$ if $m^{+}(x\mid\theta)\geq m^{-}(x\mid\theta)$
and $-1$ otherwise. A Bayesian interpretation of this rule assumes the data is
generated according to an (unknown) distribution $\pi(x,y)$, with a likelihood
function $\bar{\pi}(x|y)$, and a uniform prior $\mathsf{Unif}[{\mathcal{Y}}]$.
The functions $m^{+}(x\mid\theta)$ and $m^{-}(x\mid\theta)$ for the positive
and negative classes can be converted into posterior probability estimates:

$\displaystyle\hat{\pi}(y=1\mid x)$
$\displaystyle=\frac{\exp(m^{+}(x\mid\theta))}{\exp(m^{+}(x\mid\theta))+\exp(m^{-}(x\mid\theta))}.$
(1) $\displaystyle\hat{\pi}(y=-1\mid x)$
$\displaystyle=\frac{\exp(m^{-}(x\mid\theta))}{\exp(m^{+}(x\mid\theta))+\exp(m^{-}(x\mid\theta))}.$
(2)

The prediction function is then the maximum a posteriori (MAP) estimate:

$\displaystyle\hat{y}(x;\theta)$
$\displaystyle=\operatorname{sgn}\left\lparen\hat{\pi}(y=1\mid
x)-\hat{\pi}(y=-1\mid x)\right\rparen=\operatorname{sgn}\left\lparen
m^{+}(x\mid\theta)-m^{-}(x\mid\theta)\right\rparen.$ (3)

We can then assume the learned function is
$m(x;\theta)=m^{+}(x\mid\theta)-m^{-}(x\mid\theta)$ which is an approximation
to the (unknown) log-likelihood ratio
$\log\frac{\bar{\pi}(x|y=1)}{\bar{\pi}(x|y=-1)}$. We therefore refer to
$m^{+}(x\mid\theta)$ and $m^{-}(x\mid\theta)$ as _logits_ and $m(x;\theta)$ as
the _logit gap_. In our setting, for a given $\theta\in\varTheta$, the network
computes the function $m(x;\theta)$ (logit gap function), and the DNN defines
a family of functions
$\mathcal{M}=\\{m(x;\theta):{\mathcal{X}}\rightarrow{\mathcal{Z}}\colon\theta\in\varTheta\\}$.
The goal of _training_ a DNN is to find a “good” setting for the parameters
$\theta$ or, equivalently, to find a “good” function $m\in\mathcal{M}$.

DNN training algorithms almost always use _stochastic optimization_ for
training. The optimization algorithms used in DNN training are approximate in
two ways. First, they will generally converge to a local minimum because the
risk minimization problem is non-convex. Second, they are usually stochastic,
meaning the estimated parameters, $\bm{\theta},$ are random variables. Thus,
two runs of the same training algorithm, on the same training set
$\mathcal{D}_{\mathrm{train}}$, can produce different functions. One natural
question is to ask how different these functions are. We can try to answer
this using a _test set_ , $\mathcal{D}_{\mathrm{test}}=\\{(x_{j},y_{j})\colon
j\in[N_{\mathrm{test}}]\\}$. If we assume $\mathcal{D}_{\mathrm{test}}$ is
drawn i.i.d. from the data distribution $\pi$, then given a trained model
$m(x;\theta)$, the set of values $\\{m(x_{j};\theta)\colon
j\in[N_{\mathrm{test}}]\\}$ is also an i.i.d. sample from a distribution on
$\mathbb{R}$ induced by $\pi$.

Let $\tau$ be the distribution of $\bm{\theta}$ corresponding to the
randomness in the training algorithm run on data set
$\mathcal{D}_{\mathrm{train}}$. Then
$\mathcal{D}_{\mathrm{param}}=\\{\theta_{1},\theta_{2},\ldots,\theta_{M}\\}$
is sampled i.i.d. $\sim\tau$. Let $\mathbbm{1}\left\lparen.\right\rparen$
denote the indicator function, such that $\mathbbm{1}\left\lparen x\leq
t\right\rparen=1$ if $x\leq t$, and $0$ otherwise. Let $\mathcal{P}$ be the
set of probability measures on $\mathbb{R}$ and a CDF an element of
$\mathcal{P}$. Given a fixed DNN architecture, the training set
$\mathcal{D}_{\mathrm{train}}$, and a stochastic training algorithm, the
training algorithm takes the training set $\mathcal{D}_{\mathrm{train}}$ and
uses randomization in several different ways to produce a trained model. Thus,
given the set $\mathcal{D}_{\mathrm{param}}$ sampled i.i.d. from $\tau$, these
parameters correspond to $M$ i.i.d. samples, $\\{m(x;\theta_{k}):k\in[M]\\},$
taking values in $\mathcal{M}$. When the models are applied to
$\mathcal{D}_{\mathrm{test}}$, they produce $M$ individual CDFs taking values
in $\mathcal{P}$, which can be written as,

$\displaystyle
G_{k}(t)=\mathbb{E}_{\pi|\theta_{k}}\left[\mathbbm{1}\left\lparen
m(\mathbf{x};\theta_{k})\leq t\right\rparen\right]\>\>\text{for }k\in[M].$ (4)

Without access to the underlying data distribution $\pi$, we work with the
eCDFs instead, which can be written as,

$\displaystyle\hat{G}_{k}(t)=\frac{1}{N_{\mathrm{test}}}\sum_{j=1}^{N_{\mathrm{test}}}\mathbbm{1}\left\lparen
m(x_{j};\theta_{k})\leq t\right\rparen\>\>\text{for }k\in[M].$ (5)

The distribution on $\mathcal{P}$ of the CDFs from the test set is induced by
$\tau\times\pi$, and the expected CDF is,

$\displaystyle F_{\tau\times\pi}(t)$
$\displaystyle=\mathbb{E}_{\tau\times\pi}\left[\mathbbm{1}\left\lparen
m(\mathbf{x};\bm{\theta})\leq t\right\rparen\right].$ (6)

Without access to $\tau\times\pi$, we can consider an (approximate) empirical
version of $F_{\tau\times\pi}$ if we condition Equation (6) on
$\mathcal{D}_{\mathrm{param}}$, and get the following conditional expected
CDF,

$\displaystyle\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}(t)=\mathbb{E}_{\pi|\mathcal{D}_{\mathrm{param}}}\left[\frac{1}{M}\sum_{k=1}^{M}(\mathbbm{1}\left\lparen
m(\mathbf{x};\theta_{k})\leq t)\right\rparen\right]$
$\displaystyle=\frac{1}{M}\sum_{k=1}^{M}\mathbb{E}_{\pi|\mathcal{D}_{\mathrm{param}}}\left[(\mathbbm{1}\left\lparen
m(\mathbf{x};\theta_{k})\leq t)\right\rparen\right],$ (7)

by the linearity of expectation.
To remove the dependence on the underlying data distribution $\pi$, we further
look at an empirical version of $\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$,
by taking the average of eCDFs of $M$ individually trained models,

$\hat{G}(t)=\frac{1}{M}\sum_{k=1}^{M}\left\lparen\frac{1}{N_{\mathrm{test}}}\sum_{j=1}^{N_{\mathrm{test}}}\mathbbm{1}\left\lparen
m(x_{j},\theta_{k})\leq t\right\rparen\right\rparen.$ (8)

In contrast with $\hat{G}$, which is the average of individual eCDFs, we can
also sample $P$ models from our collection of $M$ models and construct an
ensemble model, $\bar{m}(x_{j})=\frac{1}{P}\sum_{k=1}^{P}m(x_{j},\theta_{k})$,
by averaging their logit gaps. The eCDF of this ensemble model can then be
written as,

$\hat{G}_{e}(t)=\frac{1}{N_{\mathrm{test}}}\sum_{j=1}^{N_{\mathrm{test}}}\mathbbm{1}\left\lparen\frac{1}{P}\sum_{k=1}^{P}m(x_{j},\theta_{k})\leq
t)\right\rparen,k\in M.$ (9)

Given the test set, our central object of inquiry is the empirical
distribution of logit gaps as defined in (5). We can control which sources of
randomness are used in the training by switching on or off these different
sources. For example, we can use deterministic initialization or fixed batch
ordering. Under these scenarios, we can generate $M$ models and analyze the
variability of these models with respect to (w.r.t.) a reference function,
representative of what the group has learned in consensus.

##  3 Kolmogorov-Smirnov test to analyze model variability

We use a framework from non-parametric hypothesis testing to measure how
different the models are. Ideally, we would want to analyze variability in the
training procedure by measuring the closeness of the eCDF of a trained model
$\hat{G}_{k}$ defined in Equation (5) to their expected CDF
$F_{\tau\times\pi}$ defined in Equation (6). However, without access to
$\tau\times\pi$, we consider the expected CDF
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$ conditioned over
$\mathcal{D}_{\mathrm{param}}$ as defined in Equation (2). If we had access to
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$, we could assess whether or not a
given model $m(x;\theta_{k})$ is a representative sample from the training
algorithm by testing if the model’s eCDF $\hat{G}_{k}$ is close to
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$. This can be achieved through the
following hypothesis test,

$\displaystyle\mathcal{H}_{0}$
$\displaystyle\colon\\{m(x_{j};\theta_{k})\colon
j\in[N_{\mathrm{test}}]\\}\sim\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$
$\displaystyle\mathcal{H}_{1}$
$\displaystyle\colon\\{m(x_{j};\theta_{k})\colon
j\in[N_{\mathrm{test}}]\\}\not\sim\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}},$
(10)

where $\mathcal{{H}}_{0}$, and $\mathcal{{H}}_{1}$ are the null, and alternate
hypothesis. However, in our application, we also do not have access to
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$. Instead, we can use the average
of eCDFs from $M$ models, $\hat{G}$, defined in Equation (8), as a proxy to
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$. Thus, our new goal is to test
the closeness of $\hat{G}_{k}$ to $\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$
through a hypothesis test using $\hat{G}$. In particular, we want to show if
$\hat{G}_{k}$ is close to $\hat{G}$ holds with high probability (w.h.p.), then
$\hat{G}_{k}$ will also be “close” to
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$, provided w.h.p.
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$ is close to $\hat{G}$. Thus,
without direct access to $\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$, the KS-
test to detect whether $m(x;\theta_{k})$ has been drawn from
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$ requires three steps,

  * •

Step 1: We show that the distance between
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$ and $\hat{G}$ can be bounded
w.h.p. using the DKW inequality and a union bound over the collection of
parameters.

  * •

Step 2: We choose $\hat{G}$ as the reference function and apply a two-sample
KS-test to determine if samples from $\hat{G}$ and $\hat{G}_{k}$ have been
drawn from the same distribution.

  * •

Step 3: Finally, we conclude that if samples from $\hat{G}$ and $\hat{G}_{k}$
have indeed been drawn from the same distribution, i.e., if $\hat{G}_{k}$ is
close to $\hat{G}$, then $\hat{G}_{k}$ is also close to
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$, by combining the thresholds in
Step 1 and Step 2.

To show that Step 1 holds, in section 7.1 of the Appendix, we develop an
analysis that results in the following theorem,

######  Theorem 1.

Let $\hat{G}$ and $\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$ be given by
Equation (8) and (2). Then for any $\delta_{b}>0$,

$\displaystyle\mathbb{P}_{\pi|\mathcal{D}_{\mathrm{param}}}\left(\left\lVert\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}-\hat{G}\right\rVert_{\infty}>\delta_{b}\right)\leq\epsilon_{b},$
(11)

,where $\epsilon_{b}=2M\exp(-2N_{\mathrm{test}}\delta_{b}^{2})$.

After proving Step 1, we proceed to Step 2, where we would like to test the
binary hypothesis whether samples from $\hat{G}_{k}$ and $\hat{G}$ have been
drawn from the same distribution through the following two-sample hypothesis
test,

$\displaystyle\mathcal{\tilde{H}}_{0}$
$\displaystyle\colon\hat{G},\hat{G}_{k}\text{\ are the same distribution}$
$\displaystyle\mathcal{\tilde{H}}_{1}$
$\displaystyle\colon\hat{G},\hat{G}_{k}\text{\ are not the same
distribution},$ (12)

where $\mathcal{\tilde{H}}_{0}$, and $\mathcal{\tilde{H}}_{1}$ are the null,
and alternate hypothesis. A two-sample KS test can solve the above problem,
with the two-sample DKW inequality setting a threshold for the test [33],

$\displaystyle\prob_{\pi|\mathcal{D}_{\mathrm{param}}}*{\left\lVert\hat{G}-\hat{G}_{k}\right\rVert_{\infty}>\delta_{a}}\leq\epsilon_{a},$
(13)

where $\epsilon_{a}=C\exp({-2(\frac{mn}{m+n})\delta_{a}^{2})}$, $m,n$ are the
number of i.i.d. samples from $\hat{G}$ and $\hat{G}_{k}$ respectively, and
$m=n=\frac{N_{\mathrm{test}}}{2}$. For $m=n\geq 458$, the leading constant $C$
can be replaced with $2$, as shown by Wood et. al. [33, Theorem c]. By fixing
$\epsilon_{a}$ we can set the threshold as,

$\displaystyle\delta_{a}=\sqrt{\frac{2}{N_{\mathrm{test}}}\ln{\frac{2}{\epsilon_{a}}}}.$
(14)

For models not rejected under $\mathcal{\tilde{H}}_{0}$, we can use the
triangle inequality to find an upper bound on the $L_{\infty}$-distance
between $\hat{G}_{k}$ and $\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$
resulting in,

$\displaystyle\left\lVert\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}-\hat{G}_{k}\right\rVert_{\infty}\leq\delta_{a}+\delta_{b}.\>\>\>\>$
(15)

Thus, to test our hypothesis in Equation $\eqref{eq: hypothesis test}$, we
carry out the following two-sample KS test,

$\displaystyle\lVert\hat{G}-\hat{G}_{k}\rVert_{\infty}\underset{\mathcal{\tilde{H}}_{0}}{\overset{\tilde{H}_{1}}{\gtrless}}\delta_{a}.$
(16)

which implies the following threshold for our original hypothesis,

$\displaystyle\lVert\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}-\hat{G}_{k}\rVert_{\infty}\underset{\mathcal{H}_{0}}{\overset{\mathcal{H}_{1}}{\gtrless}}\delta_{a}+\delta_{b}.$
(17)

##  4 Two-sample robust KS test to estimate model variability

In large sample settings, the KS test often rejects the null because even
small changes in the sample can result in a significant shift in the
$L_{\infty}$-norm [12, p.245]. Note that due to a lack of knowledge about our
null distribution or possible ambiguity regarding what the true
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$ is, we also do not expect our
empirical samples to look exactly like they were drawn from
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$. In this scenario, a classical KS
test will likely reject the null even for some samples drawn from
$\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$. To reduce the number of false
rejections, we relax the test to allow for a controlled number of outliers
using ideas from _robust statistics_ [17].

Consider the following scenario. Let $\hat{G}_{k}$ be an independently trained
model (candidate model) not rejected under the null hypothesis of our
classical KS test in Equation (16) i.e. w.h.p.
$\lVert\hat{G}-\hat{G}_{k}\rVert_{\infty}\leq\delta_{a}$. An upper bound on
the $L_{1}$-distance between $\hat{G}$ and $\hat{G}_{k}$ can be established
as,

$\displaystyle\lVert\hat{G}-\hat{G}_{k}\rVert_{1}\leq
N_{\mathrm{test}}\delta_{a}.$ (18)

Thus, the null hypothesis of our classical KS test corresponds to an
$L_{1}$-ball of radius $N_{\mathrm{test}}\delta_{a}$ around $\hat{G}$ as shown
in Figure 2. Consider another candidate model $\hat{G}_{j}:,j\neq k$, rejected
by the KS test in Equation (16). For a small $\delta$, the
$L_{\infty}$-distance between $\hat{G}$ and $\hat{G}_{j}$ can be bounded by:

$\displaystyle\delta_{a}\leq\lVert\hat{G}-\hat{G}_{j}\rVert_{\infty}\leq\delta_{a}+\delta.$
(19)

The $L_{1}$-distance between $\hat{G}$ and $\hat{G}_{j}$ can be similarly
bounded by:

$\displaystyle
N_{\mathrm{test}}\delta_{a}\leq\lVert\hat{G}-\hat{G}_{j}\rVert_{1}\leq
N_{\mathrm{test}}(\delta_{a}+\delta).$ (20)

Thus, with a robust KS test, we aim to broaden our null hypothesis by
including nearby models like $\hat{G}_{j}$ that are rejected under the null
but within an $L_{1}$-ball slightly greater than
$N_{\mathrm{test}}\delta_{a}$. The concept of nearby models was quantified by
Huber through contamination neighborhoods [16]. Given a distribution
$P_{0}\in\mathcal{P}$ and a scalar, $\alpha\in[0,1]$, the
$\alpha$-contaminated neighborhood around $P_{0}$ is
$\mathcal{V}_{\alpha}(P_{0})=\\{(1-\alpha)P_{0}+\alpha Q\colon
Q\in\mathcal{P}\\}$. Alvarez-Esteban et al. [1] introduce a related concept of
$\alpha$-contaminated neighborhoods as the set of so-called $\alpha$-trimmings
of $P\in\mathcal{P}$ defined by,

$\displaystyle\mathcal{R}_{\alpha}(P)=\left\\{P^{*}\in\mathcal{P}\colon
P^{*}\ll P,\frac{dP}{dP^{*}}\leq\frac{1}{(1-\alpha)}P-\text{ a.s}.\right\\}.$
(21)

Trimmings are related to contamination neighborhoods as [2],

$\displaystyle P\in\mathcal{V}_{\alpha}(P_{0})\Leftrightarrow
P_{0}\in\mathcal{R}_{\alpha}(P).$ (22)

The advantage of using trimming sets to define nearby models, as shown by del
Barrio et al. [4] is that we can compute, for a given $\hat{G}_{k}$, the
closest $L_{\infty}$ approximation to $\hat{G}$, in the set of
$\alpha$-trimmings of $\hat{G}_{k}$:

$\displaystyle
d(\hat{G},\mathcal{R}_{\alpha}(\hat{G}_{k}))=\min_{\tilde{F}\in\mathcal{R}_{\alpha}(\hat{G}_{k})}\lVert\hat{G}-\tilde{F}\rVert_{\infty}.$
(23)

The core idea behind trimming is to ask if removing a small fraction of
samples (from the test set) would allow a KS test to not reject the null
hypothesis. We point the readers to Section 7.2, and 7.3 in the Appendix for a
detailed description.

Thus, we want to know how much trimming is needed for nearby models like
$\hat{G}_{j}$, rejected under a classical KS test, i.e. not included within
$L_{1}$-ball of radius $N_{\mathrm{test}}\delta_{a}$ around $\hat{G}$, to not
be rejected under the null hypothesis. We can think of the trimming level
$\alpha$ as the radius of another $L_{1}$-ball around $\hat{G}_{j}$.
Increasing the trimming level is the same as increasing the radius of this
ball. We want to know how big the $L_{1}$-ball around $\hat{G}_{j}$ needs to
be to intersect the $L_{1}$-ball of radius $N_{\mathrm{test}}\delta_{a}$
around $\hat{G}$. This is demonstrated in a toy diagram in Figure 2. The pink
ball is an $L_{1}$-ball (trimming set) of radius $\alpha^{\prime}$ around
$\hat{G}_{j}$. A trimming level of $\alpha=\alpha^{\prime}$ wasn’t enough to
include the trimmed function $\tilde{F}$ inside $N_{\mathrm{test}}\delta_{a}$
contamination ball of $\hat{G}$. But if we increase the radius to create an
$L_{1}$-ball of radius $\alpha=\hat{\alpha}$ around $\hat{G}_{j}$, we see that
the optimal trimmed function $\tilde{F}$ is included within the
$N_{\mathrm{test}}\delta_{a}$ contamination ball of $\hat{G}$. Thus our goal
is to increment $\alpha$ and at each increment run the following hypothesis
test,

$\displaystyle
d(\hat{G},\mathcal{R}_{\alpha}(\hat{G}_{k}))\underset{\mathcal{\tilde{H}}_{0}}{\overset{\tilde{H}_{1}}{\gtrless}}\delta_{a}.$
(24)

i.e. if $\tilde{F}$ is within $N_{\mathrm{test}}\delta_{a}$ ball of $\hat{G}$.

Figure 2: Illustration of $L_{1}$-contamination balls and trimming sets

Proposed two-sample robust KS test: To that end, we develop a two-sample
version of the robust KS test discussed in Section 4. We use bootstrap
sampling to resample $\frac{N_{\mathrm{test}}}{2}$ i.i.d. test points to
compute $\hat{G}$, and $\hat{G}_{k}$. Since we only have access to $\hat{G}$
and not samples from it, let $\frac{N_{\mathrm{test}}}{2}$ be the number of
i.i.d. samples generated from $\hat{G}$ using inverse transform sampling, and
$z=z_{1},z_{2},....,z_{N_{\mathrm{test}}}$ be the combined ordered arrangement
of $N_{\mathrm{test}}$ samples from $\hat{G}$, and $\hat{G}_{k}$. For a fixed
$\alpha$, we calculate the $\alpha$-trimming of $\hat{G}_{k}$ that minimizes
the $L_{\infty}$-distance to $\hat{G}$. Our robust hypothesis test is then,

$\displaystyle\max_{z}\min_{\tilde{F}\in\mathcal{R}_{\alpha}(\hat{G}_{k})}|\hat{G}(z)-\tilde{F}(z)|\underset{\mathcal{\tilde{H}}_{0}}{\overset{\tilde{H}_{1}}{\gtrless}}\delta_{a},$
(25)

where $\delta_{a}$ can be set as discussed in Equation (13). We use this test
in Algorithm 1 to compute the smallest $\alpha$ for which the test fails to
reject $\mathcal{\tilde{H}}_{0}$. The output $\hat{\alpha}$ is our proposed
measure of discrepancy between model $\hat{G}_{k}$ and $\hat{G}$. We expect
candidate models with small $\hat{\alpha}$ to be similar to those not rejected
under the null hypothesis, i.e., the eCDF of these models are expected to all
be close to the average of their eCDFs. In Section 5 we demonstrate how
$\hat{\alpha}$ relates to metrics like accuracy and churn, commonly adopted to
understand model performance.

Data: $\mathcal{D}_{\mathrm{test}}$ samples, trained models $\\{m_{k}\colon
k\in[M]\\}$, candidate model $m_{l}\notin[M]$, threshold $\delta_{a}$,
trimming levels $(\alpha_{1},\ldots,\alpha_{T})$

Result: trimming level estimate $\hat{\alpha}$

for  _$b=1$ to $B$_ do

Compute $\hat{G}$ using $\\{m_{k}\colon k\in[M]\\}$, and
$\frac{N_{\mathrm{test}}}{2}$ resampled from $\mathcal{D}_{\mathrm{test}}$;

Compute $\hat{G}_{k}$ using $m_{l}\notin[M]$, and
$\frac{N_{\mathrm{test}}}{2}$ resampled from $\mathcal{D}_{\mathrm{test}}$;

Generate $\frac{N_{\mathrm{test}}}{2}$ samples from $\hat{G}$ using inverse
transform sampling to form $z=z_{1},z_{2},....,z_{N_{\mathrm{test}}}$ as the
combined ordered arrangement of $N_{\mathrm{test}}$ samples from $\hat{G}$,
and $\hat{G}_{k}$.

$\mathsf{Reject}\leftarrow 1$, $\hat{\alpha}\leftarrow 0$, $t\leftarrow 1$;

while  _$\mathsf{Reject}=1$ and $t\leq T$_ do

$\alpha\leftarrow\alpha_{t}$;

$\mathcal{H}\leftarrow$ output of (25) ;

if  _$\mathcal{H}=\mathcal{\tilde{H}}_{0}$_ then

$\mathsf{Reject}\leftarrow 0$, $\hat{\alpha}_{b}\leftarrow\alpha$

else

$t=t+1$;

$\hat{\alpha}\leftarrow\frac{1}{B}\sum_{b=1}^{B}\hat{\alpha}_{b}$;

Algorithm 1 Estimate $\hat{\alpha}$ measure

###  4.1 Commonly used metrics to analyze model variability

A standard way to analyze variability is by using commonly used metrics like
accuracy and churn. The _accuracy_ of a model is,

$\displaystyle
A(\theta)=\frac{1}{N_{\mathrm{test}}}\sum_{j=1}^{N_{\mathrm{test}}}\mathbbm{1}\left\lparen\hat{y}(x_{j};\theta)=y_{j}\right\rparen.$
(26)

The _churn_ is defined by

$\displaystyle
C(\theta_{1},\theta_{2})=\frac{1}{N_{\mathrm{test}}}\sum_{j=1}^{N_{\mathrm{test}}}\mathbbm{1}\left\lparen\hat{y}(x_{j};\theta_{1})\neq\hat{y}(x_{j};\theta_{2})\right\rparen,$
(27)

which is the fraction of training points where the models disagree. Both
accuracy and churn focus on the _predictions_ made by models and do not use
information about the logit gap function $m(x;\theta_{k})$ beyond its sign.
Looking at $m(x;\theta_{k})$ directly gives us other approaches to assess
whether models are similar or not: two models may have similar accuracy and
low churn but can have very different logit gaps.

We also consider a third metric for the qualitative assessment of these
models, the Expected Calibration Error (ECE) [26]. The ECE is a summary
statistic of model calibration which measures the difference in accuracy and
expected confidence and is obtained by partitioning the predictions into $R$
equally spaced bins $B_{r}$,

$\displaystyle\text{ECE}(\theta)=\sum_{r=1}^{R}\frac{\left\lvert
B_{r}\right\rvert}{N_{\mathrm{test}}}\left\lvert
A(B_{r};\theta)-\text{CONF}(B_{r};\theta)\right\rvert,$ (28)

where $A(B_{r};\theta)=\frac{1}{\left\lvert B_{r}\right\rvert}\sum_{j\in
B_{r}}\mathbbm{1}\left\lparen\hat{y}(x_{j};\theta)=y_{j}\right\rparen$, and
$\text{CONF}(B_{r};\theta)=\frac{1}{\left\lvert B_{r}\right\rvert}\sum_{j\in
B_{r}}\hat{\pi}(y_{j}=\hat{y}(x_{j};\theta)|x_{j})$. A perfectly calibrated
model will have $A(B_{r};\theta)=\text{CONF}(B_{r};\theta).$ Although
Niculescu-Mizil et al. [27] noted that DNN are well-calibrated on binary
classification tasks, it is impossible for a model to achieve perfect
calibration in reality, making this another useful metric to assess the
quality of models produced through different random seeds in conjunction with
accuracy.

##  5 Experiments

Our proposed measure of model closeness/discrepancy offers new insight into
questions around neural network training. In this section, we perform a series
of illustrative experiments to demonstrate the utility of this measure:

  * •

In Section 5.1, we show that the eCDF of a deep ensemble model reliably
approximates the reference as the size of the pool of candidate models in the
ensemble increases. Establishing this relationship at the very onset allows us
to qualitatively compare candidate models with the deep ensemble substitute of
our reference function, while comparing candidate models to the reference
function through our robust hypothesis test.

  * •

In Section 5.1.2, we use our proposed measure of model closeness
$\hat{\alpha}$ to show the minimum number of candidate models needed to form a
reliable deep ensemble that has less variability in different performance
metrics and approximates our reference function well. We show that the
discrepancy measured by the $\alpha$-trimming levels computed by our proposed
test between an ensemble model and the average of the eCDFs of those same
models decreases as the number of models in the ensemble pool increases.

  * •

In Section 5.2, and 5.3, we connect the proposed discrepancy measure to other
metrics used to measure network variability. We show how $\hat{\alpha}$ is
more informative than validation accuracy alone through two case studies. The
first uses a small CNN to perform a binary classification task on the CIFAR-10
dataset. The second uses a ViT variant, pre-trained on ImageNet, to perform
the same task on the same dataset. The former case study allowed us to use
neural networks with relatively few parameters, so we can test the limits of
the expressivity of the proposed test and the latter to show how the same test
applies to large-scale pre-trained models as well.

For experiments in Section 5.1, and 5.2, we chose the following experimental
setup. We used a small convolutional neural network with two convolutional
layers (having $32$ and $16$ features, respectively, with a $3\times 3$ kernel
size) followed by one hidden layer of 64 units and a final layer of 2 units
that output the raw logits of the network. Details of the network architecture
and parameters are included in Section 7.4.1 of Appendix. We train this
network on a subset of the CIFAR-10 dataset [19], under the following setting,

  * •

Out of the 10 classes, we take 8 and create a binary classification problem by
merging them into two super groups, with class 1 comprising of airplane,
automobile, ship, truck, and class -1 comprising of bird, dog, frog and horse.

  * •

Training set size $N_{\mathrm{train}}=40000$ and test size
$N_{\mathrm{test}}=8000$.

  * •

We fix all hyper-parameters other than the random seeds by training all models
for 50 epochs, with a fixed learning rate of 0.001, and fix the batch size to
32.

We chose a small example to allow us to train many models so we can explore
model training in different scenarios. All training runs were performed using
the Amarel HPC clusters provided by Rutgers [31]. Code to compute the test
statistic and $\hat{\alpha}$ metric of the two-sample version of the robust KS
test has been made available
online.111https://github.com/kanupolo10/SIGNAL/experiments/Robust-KS-Tests

###  5.1 The relation between our reference function and a deep ensemble

This section demonstrates how the eCDF of a deep ensemble model $\hat{G}_{e}$
defined in Equation (9), closely approximates our reference function $\hat{G}$
defined in Equation (8). Establishing this relation will allow us to compare
candidate models with the deep ensemble substitute of the reference function
using different performance metrics like validation accuracy, churn w.r.t. a
deep ensemble, and the expected calibration error (ECE) while comparing the
eCDF of candidate models to the reference function through the robust KS test.
A lower ECE indicates a better-calibrated model, while lower churn w.r.t. an
ensemble indicates less disagreement between candidate models and an ensemble.
We point the readers to Section 4.1 for a formal definition of these metrics.

Figure 3:  (Left) Histogram of logit gaps from the ensemble model with the
upper and lower envelopes representing the maximum and minimum probability
attained in each bin among individual candidate models. (Right) A plot showing
the evolution of validation accuracy of candidate models over different
epochs. The solid red dots represent the mean of validation accuracy over 800
seeds at each epoch, the light-colored area represents one standard deviation,
and the purple area represents the maximum and minimum values at that epoch.

We train a total of $1600$ models by randomly fixing a seed for each model.
The seeds control both initialization and batch order during SGD. We use the
first $M=800$ models to create a reference function $\hat{G}$ and set the rest
$M^{\prime}=800$ models to create the eCDF of a deep ensemble $\hat{G}_{e}$.
We choose different ensemble sizes $P\in[3,5,10,20,30,70,100,150,200]$ to
observe the variability of ensembles in different metrics and to understand
how large of an ensemble we need to approximate the reference function well.
For each value of $P$, we sample $P$ models without replacement from the
remaining $M^{\prime}=800$ models and compute $\hat{G}_{e}$. We repeat this
experiment $500$ times to get $500$ ensemble models using $P$ components for
each value of $P$. Deep ensemble predictors have been widely used in the
literature to reduce variability in DNN models [20, 21]. Candidate models will
have varying degrees of certainty on individual test points. Taking the
average of model “confidences” across independent training runs makes them
closer to their expected values, lowering variability. Although ensembling
through averaging over softmax probabilities is common practice in the
literature, averaging over logits has also been investigated to address the
shortcomings of probability averaging [15, 34].

Figure 3 (Left) shows how the logit gap samples obtained from ensembling
$P=M^{\prime}$ candidate models compare with candidate models in the pool. The
ensemble model produces fewer samples with small logit gaps (samples with
higher uncertainty) and large logit gaps (overconfident samples).

Figure 3 (Right) shows the evolution of validation accuracy of candidate
models over epochs. The solid red dots in the plot correspond to the mean
accuracy over $M^{\prime}$ seeds at each epoch, the light-colored region
corresponds to one standard deviation, and the purple region corresponds to
the minimum and maximum accuracy at that epoch. As observed in the plot,
validation accuracy stops increasing from epoch 30 onwards, hence the decision
to stop training at epoch 50, which is well past this optimization
convergence. The same strategy was adopted by Picard to study deep net
variability; see [28, Section 4.1] for a discussion on this.

####  5.1.1 Closeness of a deep ensemble to the reference function

The first question we want to answer is how quickly does the eCDF of an
ensemble model approximate the chosen reference? In other words, as $P\text{
approaches }M^{\prime}$ how quickly does
$\left\lVert\hat{G}-\hat{G}_{e}\right\rVert_{\infty}$ decrease? To answer
this, we compute the $L_{\infty}$-distance between $\hat{G}$ and $\hat{G}_{e}$
and compare it to the threshold set by the test in Equation (13). To set the
threshold $\delta_{a}$ in Equation (14), we fix $\epsilon_{a}=0.01$.

Figure 4 demonstrates this convergence behavior. The eCDF of all ensemble
models $\hat{G}_{e}$, produced with $P=100,150,\text{and }200$ candidate
models have $L_{\infty}$-distances from $\hat{G}$ less than the threshold set
for our KS test. This empirical evidence indicates that a deep ensemble eCDF
more closely resembles the reference function $\hat{G}$ as $P\text{ approaches
}M^{\prime}$. Therefore, we expect any candidate model close to our reference
function $\hat{G}$, formed with $M$ candidate models, to be also close to a
deep ensemble, formed with $P$ candidate models, in terms of different
performance metrics, provided $P\text{ approaches }M^{\prime}=M$. When
analyzing the variability of candidate models, we use the robust KS test for
comparing the closeness of the eCDF of each candidate model to the reference
function $\hat{G}$, and use the deep ensemble substitute of this reference
function for qualitative comparison with candidate models.

Figure 4: A plot showing $L_{\infty}$-distance of the eCDF of ensemble models
$\hat{G}_{e}$ (formed with $P$ models) w.r.t. the reference $\hat{G}$ (formed
with $M$ models) against validation accuracy of ensemble models. For each
value of $P$, we choose $P$ candidate models to form one ensemble model and
repeat this experiment 500 times through sampling with replacement to create
500 ensemble models. Thus, each dot with a fixed color represents one out of
these 500 ensemble models and each color indicates the value of $P$ or the
number of candidate models in the pool. We chose $\epsilon=0.01$, to compute
the threshold $\delta_{a}$ in Equation (14).

Since the benefits of ensembling are observed even for pool sizes as small as
$P=5$ [20], it might be tempting to focus on a few random seeds to generate a
deep ensemble model due to computational constraints. This brings us to the
next question: Does averaging over a small number of candidate models to
create an ensemble model result in a reliable representative of the training
procedure? We try to understand this qualitatively. We investigate three
metrics of ensemble models as $P$ approaches $M^{\prime}$: validation
accuracy, churn w.r.t. a deep ensemble of pool size $P=M$, and ECE. In Figure
5, we use a 3-D scatter plot to visualize the relation among the ensemble
models w.r.t. these three metrics. We expect the best ensembles to have high
validation accuracy, low churn, and low ECE.

Figure 5: 3-D scatter plot to visualize the relationship among ensemble models
in terms of validation accuracy, ECE, and churn w.r.t. an ensemble of size
$P=M$. For each value of $P$, we choose $P$ candidate models to form one
ensemble model and repeat this experiment 500 times through sampling with
replacement to create 500 ensemble models. Thus, each dot with a fixed color
represents one out of these 500 ensemble models and each color indicates the
value of $P$ or the number of candidate models in the pool.

As observed in Figure 5, there is a cluster of ensembles that all have high
accuracy, low churn, and low ECE, which is the desired behavior. These “good”
ensembles are almost all created from a large number of models. These are also
the same models whose eCDFs have achieved a smaller $L_{\infty}$-distance from
the reference function $\hat{G}$ in Figure 4. Some ensembles created from a
small number of models may outperform those created from a large number of
models in one particular metric, but it is always at the cost of performance
w.r.t. another metric. For instance, one of the highest validation accuracy
achieved in Figure 5 is an ensemble model belonging to $P=5$. However, this
model also has a higher churn w.r.t. an ensemble of size $P=M$. This is also
reflected by the high $L_{\infty}$-distance of this model from the reference
$\hat{G}$ as shown in Figure 4. We conclude that contrary to standard
practice, exploring more than 5 seeds is important to create a reliable
ensemble that has less variability in different performance metrics. Thus,
ensembling over smaller pool sizes results in a deep ensemble that has
significant variability over different metrics used to evaluate model
performance, and increasing the size $P$ of ensemble models reduces this
variability.

####  5.1.2 Selecting the number of models to be ensembled via the robust KS
test

Section 5.1.1 showed that it is possible to create a good ensemble with enough
models that have less variability over different performance metrics and whose
eCDF approximates the reference function well. The more models we ensemble,
the better and more consistent the ensembles become. This brings us to our
final question: How many candidate models do we need to create a deep ensemble
model that is a reliable representative of the training procedure i.e. has
less variability over different performance metrics and approximates the
reference function well? Based on the threshold set by our test on the
$L_{\infty}$-distance and as demonstrated in Figure 4, a classical KS test
will indicate that we need to explore 100 random seeds to create a reliable
ensemble, which can then be used as a base model for qualitative comparisons
with candidate models. Since the classical KS test is too sensitive, setting a
threshold based on it will require too many models to create the ensemble. We
can use the proposed robust KS test to also measure the similarity between the
ensemble function $\hat{G}_{e}$ and the reference function $\hat{G}$ but with
a little wiggle room. In practice, we also may not want to rely on validation
accuracy or churn to decide on the number of models we need. So, the proposed
test offers a way to estimate this number using only the logits. To that end,
we compute $\hat{\alpha}$ through a robust KS test to measure the
closeness/discrepancy of the eCDFs of ensemble models $\hat{G}_{e}$, of
different sizes, with the reference function $\hat{G}$. Similar to Figure 4,
in Figure 6, as $P$ approaches $M^{\prime}$, $\hat{\alpha}$ approaches 0.
However, if we set a threshold on the $\alpha$-trimming level instead of the
$L_{\infty}$-distance, we get a better lower bound on the number of models we
need to use for a reliable ensemble. If we consider all ensembles that require
$\hat{\alpha}\leq 0.05$, we see that ensembles created from $P=30$ models or
more are all included in this set as shown in column 2 of Table 1, while more
than $20\%$ of ensembles created from $P=3$ models have required a higher
level of trimming. From Table 1, we notice that ensembles achieve better
values in all metrics as the size of the ensemble $P$ increases and the
variability of these three metrics reduces noticeably for ensembles created
from $P=30$ models or more. Thus, our proposed metric indicates that
ensembling over at least $P=30$ candidate models may be required for an
ensemble to be a reliable representative of the training variability.

Figure 6: A plot showing $\hat{\alpha}$ of the eCDF of ensemble models $\hat{G}_{e}$ (formed with $P$ models) w.r.t. the reference $\hat{G}$ (formed with $M$ models) against validation accuracy of ensemble models. For each value of $P$, we choose $P$ candidate models to form one ensemble model and repeat this experiment 500 times through sampling with replacement to create 500 ensemble models. Thus, each dot with a fixed color represents one out of these 500 ensemble models and each color indicates the value of $P$ or the number of candidate models in the pool. $\mathbf{P}$ |  % of models with $\mathbf{\hat{\alpha}\leq 0.05}$ |  Accuracy mean $\underset{-}{+}$ std |  Churn w.r.t. ensemble of size $\mathbf{P=M}$ mean $\underset{-}{+}$ std |  ECE mean $\underset{-}{+}$ std
---|---|---|---|---
3 |  77.20 |  91.85 $\underset{-}{+}$ 0.25 |  190.79 $\underset{-}{+}$ 33.90 |  0.0226 $\underset{-}{+}$ 0.0048
5 |  89.20 |  92.03 $\underset{-}{+}$ 0.18 |  145.70 $\underset{-}{+}$ 28.61 |  0.0200 $\underset{-}{+}$ 0.0037
10 |  96.80 |  92.11 $\underset{-}{+}$ 0.14 |  103.57 $\underset{-}{+}$ 22.50 |  0.0185 $\underset{-}{+}$ 0.0023
20 |  99.40 |  92.17 $\underset{-}{+}$ 0.11 |  76.972 $\underset{-}{+}$ 17.23 |  0.0175 $\underset{-}{+}$ 0.0016
30 |  100.00 |  92.21 $\underset{-}{+}$ 0.10 |  62.84 $\underset{-}{+}$ 12.45 |  0.0171 $\underset{-}{+}$ 0.0014
70 |  100.00 |  92.22 $\underset{-}{+}$ 0.08 |  43.45 $\underset{-}{+}$ 9.64 |  0.0168 $\underset{-}{+}$ 0.0012
100 |  100.00 |  92.23 $\underset{-}{+}$ 0.08 |  38.59 $\underset{-}{+}$ 9.08 |  0.0169 $\underset{-}{+}$ 0.0012
150 |  100.00 |  92.23 $\underset{-}{+}$ 0.07 |  32.75 $\underset{-}{+}$ 7.42 |  0.0169 $\underset{-}{+}$ 0.0011
200 |  100.00 |  92.24 $\underset{-}{+}$ 0.07 |  29.3 $\underset{-}{+}$ 6.91 |  0.0169 $\underset{-}{+}$ 0.0011
Table 1: Table showing ensemble statistics for different values of $P$

###  5.2 Evaluating our proposed metric of model closeness/discrepancy

We next move on to understand how our proposed measure of model
closeness/discrepancy relates to metrics commonly used to understand model
variability. We try to understand this for candidate models generated under
different sources of randomness ( like initialization and batch order). In
particular, we considered three scenarios: $\mathsf{S}_{\mathrm{init}}$ with
only random initialization and fixed batches, $\mathsf{S}_{\mathrm{batch}}$
with only random batch selection in SGD and fixed initialization, and
$\mathsf{S}_{\mathrm{all}}$ combining both sources of randomness. We trained
200 models in each of $\mathsf{S}_{\mathrm{init}}$,
$\mathsf{S}_{\mathrm{batch}}$, and $\mathsf{S}_{\mathrm{all}}$ with the same
hyper-parameter settings as detailed at the beginning of Section 5. For each
source of randomness, we use $M=100$ models to create the reference function
$\hat{G}$ as defined in Equation $\eqref{eq:avg_cdf}$ and choose our candidate
models with eCDF $\hat{G}_{k}$ as defined in Equation $\eqref{eq:Individual
ecdf}$ from the rest $M^{\prime}=100$ models. For each candidate model, we use
its eCDF $\hat{G}_{k}$ to compute $\hat{\alpha}$ through a robust KS test
against the reference $\hat{G}$.

####  5.2.1 Comparing our proposed discrepancy measure with accuracy

Figure 7: 2-D Scatter plot to visualize the relationship of candidate models
with their average cdf, $\hat{G}$, in terms of validation accuracy and
$\hat{\alpha}$ for models in $\mathsf{S}_{\mathrm{init}}$,
$\mathsf{S}_{\mathrm{batch}}$ and $\mathsf{S}_{\mathrm{all}}$.

Our first question in this section is: How does our proposed metric
$\hat{\alpha}$ relate to accuracy? In Figure 7, we plot the relationship
between the validation accuracy of candidate models and $\hat{\alpha}$ for
candidate models in $\mathsf{S}_{\mathrm{init}}$,
$\mathsf{S}_{\mathrm{batch}}$, and $\mathsf{S}_{\mathrm{all}}$. For each plot,
we notice significant variability over random seeds in both $\hat{\alpha}$ and
validation accuracy for the same hyper-parameter setting. This is evidence
that a fixed hyper-parameter setting that works well for one random seed can
perform poorly for another and thus produce models with very different logit
gap functions. We notice more variability in both validation accuracy and
$\hat{\alpha}$ among models in $\mathsf{S}_{\mathrm{batch}}$, than models in
$\mathsf{S}_{\mathrm{init}}$. Details of accuracy statistics for each source
of randomness are listed in Table 2.

Random Sources |  Accuracy mean $\underset{-}{+}$ std |  Accuracy range |  Ensemble Accuracy |  Accuracy range of models with $\mathbf{\hat{\alpha}\leq 0.05}$
---|---|---|---|---
$\mathsf{S}_{\mathrm{init}}$ |  91.101 $\underset{-}{+}$ 0.302 |  [90.422, 91.781] |  91.731 |  [90.422, 91.781]
$\mathsf{S}_{\mathrm{batch}}$ |  90.819 $\underset{-}{+}$ 0.554 |  [87.852, 91.843] |  91.731 |  [90.409, 91.843]
$\mathsf{S}_{\mathrm{all}}$ |  90.810 $\underset{-}{+}$ 0.560 |  [88.513, 91.769] |  92.361 |  [90.409, 91.769]
Table 2: Validation Accuracy statistics of candidate models generated from
different random sources:
$\mathsf{S}_{\mathrm{init}}$,$\mathsf{S}_{\mathrm{batch}}$, and
$\mathsf{S}_{\mathrm{all}}$ along with the ensemble accuracy for each random
source.

The mean accuracy of models in $\mathsf{S}_{\mathrm{init}}$ is slightly higher
than models in $\mathsf{S}_{\mathrm{all}}$, and $\mathsf{S}_{\mathrm{batch}}$,
while all three categories have achieved similar maximum accuracy. The higher
variability among models in $\mathsf{S}_{\mathrm{batch}}$ comes from producing
models with much poorer validation accuracy than $\mathsf{S}_{\mathrm{init}}$,
a trend also observed in $\mathsf{S}_{\mathrm{all}}$ as reported in the
minimum accuracy of the accuracy range for each category. This indicates that
random batch shuffling has a higher effect on the overall variability among
models in $\mathsf{S}_{\mathrm{all}}$, than random initialization. We also
notice that although there is more variability among models in
$\mathsf{S}_{\mathrm{batch}}$, than $\mathsf{S}_{\mathrm{init}}$, the ensemble
model has performed the same for both categories. However, the combined
variability of $\mathsf{S}_{\mathrm{init}}$, and
$\mathsf{S}_{\mathrm{batch}}$, reflected among models in
$\mathsf{S}_{\mathrm{all}}$ has led to a better ensemble performance in
$\mathsf{S}_{\mathrm{all}}$.

We also look at models admitted by small values of $\hat{\alpha}$ in each
category (e.g. $\hat{\alpha}\leq 0.05$), and observe the corresponding range
of validation accuracy, reported in the fourth column of Table 2. For
$\mathsf{S}_{\mathrm{init}}$, the range of accuracy achieved by models
admitted by $\hat{\alpha}\leq 0.05$ covers the full range of accuracy of all
models in $\mathsf{S}_{\mathrm{init}}$. This range is also close to the range
of models admitted by $\hat{\alpha}\leq 0.05$ in $\mathsf{S}_{\mathrm{all}}$,
and $\mathsf{S}_{\mathrm{batch}}$. The maximum validation accuracy for each
source is also admitted by small values of $\hat{\alpha}$. This is because
candidate models aggregate to form ensembles that result in a boost over
average performance in the group. Since models with lower $\hat{\alpha}$ are
close to the reference function in terms of their eCDF, and hence also to the
ensemble, they will also have accuracy similar to the ensemble. For each
source of randomness, if we consider the range of validation accuracy for
models admitted by small values of $\hat{\alpha}$ to be representative of the
training variability (as a result of being close to the reference function),
then any model that has performed worse than this range will end up with a
high value of $\hat{\alpha}$. This can be observed for
$\mathsf{S}_{\mathrm{batch}}$, and $\mathsf{S}_{\mathrm{all}}$ since these two
categories have resulted in more models with performance worse than the
discussed range. However, the opposite is not true, i.e. models with accuracy
within the discussed range will not always have a small $\hat{\alpha}$, and
hence will not be good representatives. This shows that validation accuracy
alone is not the right metric to assess model quality. A model can end up
within the discussed range of validation accuracy but not be a good
representative of the training variability over random seeds. We conclude that
lower $\hat{\alpha}$ does not admit poor validation accuracy, but similar
validation accuracy does not imply similar $\hat{\alpha}$.

####  5.2.2 Comparing our proposed discrepancy measure with pairwise churn

Figure 8: 3-D scatter plot to visualize the relationship of candidate models
with their average cdf, $\hat{G}$, in terms of pairwise _churn_ and
$\hat{\alpha}$, for models in $\mathsf{S}_{\mathrm{init}}$,
$\mathsf{S}_{\mathrm{batch}}$ and $\mathsf{S}_{\mathrm{all}}$

.

Random Sources |  Range of pairwise churn |  Range of pairwise churn of pairs of models with $(\mathbf{\hat{\alpha_{i}}\leq 0.05},\mathbf{\hat{\alpha_{j}}\leq 0.05}),i\neq j$
---|---|---
$\mathsf{S}_{\mathrm{init}}$ |  [270, 481 ] |  [270, 481 ]
$\mathsf{S}_{\mathrm{batch}}$ |  [218, 1174 ] |  [ 283, 467]
$\mathsf{S}_{\mathrm{all}}$ |  [327, 1184 ] |  [283, 462]
Table 3: Pairwise churn statistics of models generated from different random
sources: $\mathsf{S}_{\mathrm{init}}$, $\mathsf{S}_{\mathrm{batch}}$, and
$\mathsf{S}_{\mathrm{all}}$

How does our proposed metric $\hat{\alpha}$ relate to pairwise churn? In
Figure 8, we plot the relationship between the pairwise disagreement between
candidate models or churn, and $\hat{\alpha}$, for candidate models in
$\mathsf{S}_{\mathrm{init}}$, $\mathsf{S}_{\mathrm{batch}}$, and
$\mathsf{S}_{\mathrm{all}}$. We distinguish this from churn w.r.t. a deep
ensemble by labeling it as pairwise churn. Similar to the trend observed for
validation accuracy, there is more variability in pairwise churn among pairs
of models in $\mathsf{S}_{\mathrm{init}}$, than $\mathsf{S}_{\mathrm{batch}}$,
as observed in the pairwise churn range for each random source reported in
column 2 of Table 3. Again, most of the variability in pairwise churn in
$\mathsf{S}_{\mathrm{all}}$ is due to $\mathsf{S}_{\mathrm{batch}}$. For all
plots in Figure 8, we observe that pairs of models with lower values of
$\hat{\alpha}$ have achieved a smaller range of pairwise churn among its
group. This is less obvious for models in $\mathsf{S}_{\mathrm{init}}$,
because this group’s total range of pairwise churn is low. We see this in
column 3 of Table 3, where we focus on pairs of models in each category that
have achieved pairwise $\hat{\alpha}$ less than $0.05$ and report the range of
pairwise churn for those models. Since the range of pairwise churn of pairs of
models in $\mathsf{S}_{\mathrm{init}}$ is low, small values of pairwise
$\hat{\alpha}$ have admitted this full range. We also see that the pairs of
models admitted by small values of $\hat{\alpha}$ in
$\mathsf{S}_{\mathrm{batch}}$, and $\mathsf{S}_{\mathrm{all}}$ have a pairwise
churn range close to the pairwise churn range admitted by small values of
$\hat{\alpha}$ in $\mathsf{S}_{\mathrm{init}}$. Any pairs of models with
pairwise churn higher than this range will have at least one model in the pair
with a high $\hat{\alpha}$. A possible explanation for a pair of models with
low pairwise $\hat{\alpha}$ also having a low pairwise churn is that, models
with low $\hat{\alpha}$ are expected to consist of ”good” models that are
similar to a good quality ensemble. Since these models achieve validation
accuracy within a higher range among their group, the number of points they
will make mistakes on will be low. This will also result in a low pairwise
churn because these models will only have a limited number of test points to
disagree on. As observed for validation accuracy, the opposite is again not
true for pairwise churn, i.e., models with low pairwise churn will not
necessarily have small values of pairwise $\hat{\alpha}$. Higher values of
$\hat{\alpha}$ will include models with both high and low validation accuracy,
resulting in either high or low pairwise churn. We conclude that  lower
pairwise $\hat{\alpha}$ does not admit high pairwise churn but low pairwise
churn does not imply lower pairwise $\hat{\alpha}$.

Figure 9: 2-D scatter plot to visualize the relationship of $\hat{\alpha}$ of
candidate models in $\mathsf{S}_{\mathrm{all}}$ with their average cdf,
$\hat{G}$, in terms of (Left) Validation accuracy, (Middle) Churn w.r.t. their
ensemble and (Right) ECE.

####  5.2.3 Our proposed discrepancy measure is more informative than accuracy

Figure 10: 3-D scatter plot to visualize how $\hat{\alpha}$ compares to other
metrics for 100 candidate CNN models performing a binary classification task
on CIFAR-10.

ECE |  Churn w.r.t. ensemble |  Val Accuracy |  Average Churn |  $\hat{\alpha}$
---|---|---|---|---
0.021 |  280 |  91.519 |  458.61 |  0.005
0.023 |  318 |  91.769 |  476.95 |  0.021
0.026 |  376 |  91.569 |  $\mathbf{513.83}$ |  $\mathbf{0.134}$
0.029 |  307 |  91.556 |  476.35 |  0.027
0.030 |  371 |  91.382 |  482.45 |  0.024
0.026 |  265 |  91.581 |  459.37 |  0.005
0.030 |  309 |  91.157 |  475.91 |  0.016
0.028 |  297 |  91.108 |  472.15 |  0.002
$\mathbf{0.038}$ |  $\mathbf{425}$ |  91.182 |  $\mathbf{552}$ |  $\mathbf{0.201}$
Figure 11: Table showing $\hat{\alpha}$ values for CNN models with similar
validation accuracy and their corresponding values achieved in other metrics.

We now provide evidence that $\hat{\alpha}$ is more informative than
validation accuracy for model comparison. Figure 9 and Figure 11 show the
relationship of $\hat{\alpha}$ of candidate models in
$\mathsf{S}_{\mathrm{all}}$ with three metrics: validation accuracy, churn
w.r.t the ensemble of candidate models, and their ECE, in 2-D and 3-D
respectively. If the eCDF of a candidate model is close to the eCDF of the
reference function, i.e. for smaller values of $\hat{\alpha}$ the variability
in each of the three metrics is confined to a range that includes the best
value in each metric, with the worst value not differing too much from the
best value. For each metric, we can think of the range of values corresponding
to a low $\hat{\alpha}$ as the minimum variability observed in models that are
close to their expected distribution. Any model that has achieved a value
worse than this range will also have a high $\hat{\alpha}$. Models with low
$\hat{\alpha}$ are within a higher accuracy range among their group and have
ECE and churn w.r.t. ensemble within a low range. However, the converse is not
true, i.e. high accuracy or low churn w.r.t. ensemble or low ECE alone does
not imply low $\hat{\alpha}$. These models can achieve a value that is within
a good range of values in one metric but can perform poorly in another metric,
and this is reflected in the corresponding $\hat{\alpha}$ value. To see this,
we focus on models that have achieved a similar high validation accuracy and
look at other metrics like ECE, average churn for each model w.r.t. all other
models, and churn w.r.t. their ensemble in Table 11. $\hat{\alpha}$ is low if
all three metrics fall within the previously discussed good range of values.
An increase in $\hat{\alpha}$ (denoted in bold) in Table 11 indicates a
reduction in quality in one or more of these metrics. We can see that most
models with similar validation accuracy in Table 11 have needed only a small
trimming level ($\hat{\alpha}\leq 0.05$) to not reject the null. But in row 3
and row 9, we see that despite having similar validation accuracy,
$\hat{\alpha}$ values are large due to a reduction in quality in other metrics
(denoted in bold).

###  5.3 Application in Transfer Learning

Figure 12:  (Left) Validation loss vs epochs of ViT models, pre-trained on
ImageNet, and performing binary classification task on CIFAR10. The solid red
dots represent the mean of validation loss over 45 seeds at each epoch, the
light-colored area represents one standard deviation, and the orange area
represents the maximum and minimum values at that epoch. (Right) Validation
accuracy vs epochs of Vit models pre-trained on ImageNet and performing binary
classification task on CIFAR10. The solid red dots represent the mean of
validation accuracy over 45 seeds at each epoch, the light-colored area
represents one standard deviation, and the purple area represents the maximum
and minimum values at that epoch.

We conclude our experiments by demonstrating the usefulness of our framework
in transfer learning. We use a Vision Transformer(ViT) variant [22], pre-
trained on ImageNet, to perform a downstream binary classification task on
CIFAR-10 (as described in Section 5). The pre-trained model is available as
part of the Hugging Face Community [23]. We only train the final task-specific
classification layer for 5 epochs (till close to convergence as visualized in
Figure 12) and rely on already fine-tuned hyper-parameters and the pre-trained
weights of the upstream task to “fine-tune” the downstream binary
classification task over random seeds only. Thus, keeping the fine-tuning
regime and pre-trained weights fixed, we only vary the random seed that
controls different sources of randomness in the training procedure
(initialization and batch order during SGD), generate 90 models, and observe
the variability over random seeds for this fixed setup. As observed by Picard
et al. [28], there is less variability in validation accuracy over random
seeds when using pre-trained models.

We create our reference function $\hat{G}$ as defined in Equation
$\eqref{eq:avg_cdf}$, using the first 45 models and treat the rest 45 models
as our candidate models $\hat{G}_{k}$ as defined in Equation
$\eqref{eq:Individual ecdf}$, for comparison with the reference. We conduct
the same experiment as in Section 5.2.3, where we run our robust hypothesis
test between the reference function $\hat{G}$ and candidate models
$\hat{G}_{k}$ to compute our proposed discrepancy measure $\hat{\alpha}$.
Following similar reasoning to Section 5.2.3, in Figure 14, we observe how our
proposed metric is more informative than validation accuracy. $\hat{\alpha}$
is indicative of the quality of models in terms of other metrics like ECE,
churn w.r.t. an ensemble and average churn of candidate models w.r.t. all
other models in the pool alongside accuracy. A higher $\hat{\alpha}$ is
usually accompanied by a reduction in quality in one or more of these three
metrics, denoted in bold in Table 14. Although pre-trained large-scale models
provide us with very high validation accuracy and less variability in
different metrics used to assess model performance when fine-tuned over the
final classification layer, our experiment demonstrates how low values of
$\hat{\alpha}$ are still useful in identifying the best representatives of the
training procedure.

Figure 13:  3-D scatter plot to visualize how $\hat{\alpha}$ compares to other
metrics for 45 candidate ViT models pre-trained on ImageNet for a binary
classification task on CIFAR-10. Each dot represents one candidate model.

Average churn |  Churn w.r.t. ensemble |  Val accuracy |  ECE |  $\hat{\alpha}$
---|---|---|---|---
48.33 |  30 |  99.3 |  0.002 |  0.026
$\mathbf{49.44}$ |  28 |  99.25 |  $\mathbf{0.003}$ |  $\mathbf{0.118}$
45.88 |  24 |  99.2 |  0.001 |  0.007
44.82 |  26 |  99.2 |  0.002 |  0.004
47.97 |  $\mathbf{34}$ |  99.2 |  0.001 |  $\mathbf{0.069}$
Figure 14: Table showing $\hat{\alpha}$ values for candidate models with
similar validation accuracy and their corresponding values achieved in other
metrics.

Through these experiments, we highlight the importance of treating seeds
controlling random elements in a DNN architecture as a separate hyper-
parameter that needs “tuning” and propose a framework to reliably select seeds
that do not rely solely on commonly used performance metrics like validation
accuracy. For practitioners interested in a good seed selection to account for
the variability caused by random seeds, the recommendation is to explore at
least 30 seeds and choose the seed that generates a model with high validation
accuracy and low $\hat{\alpha}$.

##  6 Conclusion and Future Work

In this work, we propose a framework for random seed selection of a DNN with a
fixed architecture and a fixed hyper-parameter setting. Our proposed framework
is based on a robust two-sample hypothesis testing problem that uses trimming
of the eCDF of samples obtained from the logit gap function. Our test assesses
the closeness of candidate models in a pool with their expected eCDF by down-
weighting that part of the data that has a greater influence on the
dissimilarity. We also provide some evidence that our new measure, the
trimming level $\hat{\alpha}$, could be a more informative metric to assess
model performance than commonly used test/validation accuracy. This allows us
to perform random seed selection in a more principled fashion instead of
relying solely on trial-and-error methods or metrics like validation accuracy.
Since the usefulness of our methodology relies on ensembles achieving
significant performance gains over their constituents, we investigate the
behavior of the ensemble as the number of models in the ensemble pool
increases. Our proposed metric $\hat{\alpha}$ can also be used as a stopping
criterion for the number of different seeds that need to be explored for an
ensemble model to reliably approximate the expected eCDF of the logit gap
distribution of candidate models, and to have less variability in different
performance metrics. Although ensembling by model averaging in the output
space is quite straightforward and usually results in performance gains, this
comes at the cost of reduced interpretability, since such ensembles do not
learn any parameters or features. By selecting models that are close to their
ensembles in the output space of a DNN, instead of the ensemble itself, our
framework leaves room for model interpretability while at the same time
maintaining high accuracy compared to their counterparts. Our methodology can
find usefulness in critical application areas like credit risk assessment
where a single model is often preferred over the improved accuracy of an
ensemble due to interpretability concerns [10].

In scenarios where ensembling by model averaging might not lead to significant
performance gains over candidate models in the pool, future extensions of this
work include investigating when candidate models form effective ensembles. We
are also interested in extensions to multi-class classification and exploring
robust two-sample hypothesis testing based on other distance metrics like the
Wasserstein metric [3]. In this work, we focused on samples from the logit gap
function as our probe to understand deep net variability. We can look at other
functions of the trained models, such as the Jacobian or the Neural Tangent
Kernel of functions learned by these models. The eigen-distribution of these
matrices for instance can inform us how well individual candidate models can
generalize to test data.

## References

  * [1] Pedro César Alvarez-Esteban, Eustasio Barrio, Juan Antonio Cuesta-Albertos and Carlos Matran  “Trimmed comparison of distributions”  In _Journal of the American Statistical Association_ 103.482 Taylor & Francis, 2008, pp. 697–704  URL: https://www.jstor.org/stable/27640092
  * [2] Pedro C Álvarez-Esteban, Eustasio Barrio, Juan A Cuesta-Albertos and Carlos Matrán  “Uniqueness and approximate computation of optimal incomplete transportation plans”  In _Annales de l’IHP Probabilités et statistiques_ 47.2, 2011, pp. 358–375
  * [3] Pedro C. Álvarez-Esteban, Eustasio Barrio, Juan A. Cuesta-Albertos and Carlos Matrán  “Similarity of samples and trimming”  In _Bernoulli_ 18.2 Bernoulli Society for Mathematical StatisticsProbability, 2012  DOI: 10.3150/11-bej351
  * [4] Eustasio Barrio, Hristo Inouzhe and Carlos Matrán  “On approximate validation of models: a Kolmogorov–Smirnov-based approach”  In _TEST_ 29.4 Springer, 2020, pp. 938–965
  * [5] Eustasio del Barrio, Hristo Inouzhe and Carlos Matrán  “Box-Constrained Monotone Approximations to Lipschitz Regularizations, with Applications to Robust Testing”  Springer, 2020, pp. 65–87
  * [6] Xavier Bouthillier, Pierre Delaunay, Mirko Bronzi, Assya Trofimov, Brennan Nichyporuk, Justin Szeto, Nazanin Mohammadi Sepahvand, Edward Raff, Kanika Madan, Vikram Voleti, Samira Ebrahimi Kahou, Vincent Michalski, Tal Arbel, Chris Pal, Gael Varoquaux and Pascal Vincent  “Accounting for variance in machine learning benchmarks”  In _Proceedings of Machine Learning and Systems_ 3, 2021
  * [7] Xavier Bouthillier, César Laurent and Pascal Vincent  “Unreproducible Research is Reproducible”  In _Proceedings of the 36th International Conference on Machine Learning_ 97, Proceedings of Machine Learning Research  PMLR, 2019, pp. 725–734  URL: https://proceedings.mlr.press/v97/bouthillier19a.html
  * [8] Jesse Dodge, Gabriel Ilharco, Roy Schwartz, Ali Farhadi, Hannaneh Hajishirzi and Noah Smith  “Fine-Tuning Pretrained Language Models: Weight Initializations, Data Orders, and Early Stopping”, 2020  arXiv:2002.06305 [cs.CL]
  * [9] Aryeh Dvoretzky, Jack Kiefer and Jacob Wolfowitz  “Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator”  In _The Annals of Mathematical Statistics_ JSTOR, 1956, pp. 642–669
  * [10] Raquel Florez-Lopez and Juan Manuel Ramon-Jeronimo  “Enhancing accuracy and interpretability of ensemble strategies in credit risk assessment. A correlated-adjusted decision forest proposal”  In _Expert Systems with Applications_ 42.13 Elsevier, 2015, pp. 5737–5753
  * [11] Stanislav Fort, Huiyi Hu and Balaji Lakshminarayanan  “Deep Ensembles: A Loss Landscape Perspective”, 2020  DOI: 10.48550/arXiv.1912.02757
  * [12] Jean Dickinson Gibbons and Subhabrata Chakraborti  “Nonparametric Statistical Inference (Sixth Edition)”  Boca Raton, FL, USA: CRC Press, 2021
  * [13] Odd Erik Gundersen, Kevin Coakley, Christine Kirkpatrick and Yolanda Gil  “Sources of Irreproducibility in Machine Learning: A Review”, 2023  DOI: 10.48550/arXiv.2204.07610
  * [14] Peter Henderson, Riashat Islam, Philip Bachman, Joelle Pineau, Doina Precup and David Meger  “Deep Reinforcement Learning that Matters”, 2019  arXiv:1709.06560 [cs.LG]
  * [15] Geoffrey Hinton, Oriol Vinyals and Jeff Dean  “Distilling the knowledge in a neural network”  In _arXiv preprint arXiv:1503.02531_ , 2015
  * [16] Peter J Huber  “Robust estimation of a location parameter”  In _Breakthroughs in statistics: Methodology and distribution_ Springer, 1992, pp. 492–518
  * [17] Peter J. Huber and Elvezio M. Ronchetti  “Robust Statistics, Second Edition”  Hoboken, NJ, USA: John Wiley & Sons, 2009  DOI: 10.1002/9780470434697
  * [18] Keller Jordan  “Calibrated Chaos: Variance Between Runs of Neural Network Training is Harmless and Inevitable”, 2023  DOI: 10.48550/arXiv.2304.01910
  * [19] Alex Krizhevsky  “Learning Multiple Layers of Features from Tiny Images”, 2009  URL: https://www.cs.toronto.edu/~kriz/learning-features-2009-TR.pdf
  * [20] Balaji Lakshminarayanan, Alexander Pritzel and Charles Blundell  “Simple and Scalable Predictive Uncertainty Estimation Using Deep Ensembles”  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , NIPS’17  Long Beach, California, USA: Curran Associates Inc., 2017, pp. 6405–6416
  * [21] Stefan Lee, Senthil Purushwalkam, Michael Cogswell, David Crandall and Dhruv Batra  “Why M Heads are Better than One: Training a Diverse Ensemble of Deep Networks”, 2015
  * [22] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin and Baining Guo  “Swin Transformer: Hierarchical Vision Transformer using Shifted Windows”, 2021  arXiv:2103.14030 [cs.CV]
  * [23] Ze Liu, Yutong Lin, Yue Cao, Han Hu, Yixuan Wei, Zheng Zhang, Stephen Lin and Baining Guo  “Swin Transformer: Hierarchical Vision Transformer using Shifted Windows”  In _CoRR_ abs/2103.14030, 2021  arXiv: https://arxiv.org/abs/2103.14030
  * [24] Gábor Melis, Chris Dyer and Phil Blunsom  “On the State of the Art of Evaluation in Neural Language Models”, 2017  arXiv:1707.05589 [cs.CL]
  * [25] Mahdi Milani Fard, Quentin Cormier, Kevin Canini and Maya Gupta  “Launch and Iterate: Reducing Prediction Churn”  In _Advances in Neural Information Processing Systems_ 29 Curran Associates, Inc., 2016  URL: https://proceedings.neurips.cc/paper/2016/file/dc5c768b5dc76a084531934b34601977-Paper.pdf
  * [26] Mahdi Pakdaman Naeini, Gregory Cooper and Milos Hauskrecht  “Obtaining well calibrated probabilities using bayesian binning”  In _Proceedings of the AAAI conference on artificial intelligence_ 29.1, 2015
  * [27] Alexandru Niculescu-Mizil and Rich Caruana  “Predicting good probabilities with supervised learning”  In _Proceedings of the 22nd international conference on Machine learning_ , 2005, pp. 625–632
  * [28] David Picard  “Torch. manual_seed (3407) is all you need: On the influence of random seeds in deep learning architectures for computer vision”  In _arXiv preprint arXiv:2109.08203_ , 2021
  * [29] Rahul Rahaman and Alexandre H. Thiery  “Uncertainty Quantification and Deep Ensembles”, 2021  arXiv:2007.08792 [stat.ML]
  * [30] Nils Reimers and Iryna Gurevych  “Reporting Score Distributions Makes a Difference: Performance Study of LSTM-networks for Sequence Tagging”, 2017  arXiv:1707.09861 [cs.CL]
  * [31] “Rutgers Office of Advance Research Computing”, https://oarc.rutgers.edu/resources/amarel/
  * [32] Cecilia Summers and Michael J. Dinneen  “Nondeterminism and Instability in Neural Network Optimization”, 2021  DOI: 10.48550/arXiv.2103.04514
  * [33] Fan Wei and Richard M. Dudley  “Two-sample Dvoretzky–Kiefer–Wolfowitz inequalities”  In _Statistics & Probability Letters_ 82.3, 2012, pp. 636–644  DOI: https://doi.org/10.1016/j.spl.2011.11.012
  * [34] Danny Wood, Tingting Mu, Andrew M Webb, Henry WJ Reeve, Mikel Lujan and Gavin Brown  “A unified theory of diversity in ensemble learning”  In _Journal of Machine Learning Research_ 24.359, 2023, pp. 1–49

##  7 Appendix

###  7.1 Analysis

Theorem 1. Let $\hat{G}$ and $\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}$ be
given by Equation (8) and (2). Then for any $\delta_{b}>0$,

$\displaystyle\mathbb{P}_{\pi|\mathcal{D}_{\mathrm{param}}}\left(\left\lVert\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}-\hat{G}\right\rVert_{\infty}>\delta_{b}\right)\leq\epsilon_{b},$

where $\epsilon_{b}=2M\exp(-2N_{\mathrm{test}}\delta_{b}^{2})$.
Proof of Theorem 1:

$\displaystyle\left\lVert\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}-\hat{G}\right\rVert_{\infty}$
$\displaystyle=\left\lVert\frac{1}{M}\sum_{k=1}^{M}\left\lparen\mathbb{E}_{\pi|\theta_{k}}\left[\mathbbm{1}\left\lparen
m(\mathbf{x},\theta_{k})\leq
t\right\rparen\right]-\frac{1}{N_{\mathrm{test}}}\sum_{j=1}^{N_{\mathrm{test}}}\mathbbm{1}\left\lparen
m(x_{j};\theta_{k})\leq t\right\rparen\right\rparen\right\rVert_{\infty}$ (29)
$\displaystyle\leq\frac{1}{M}\sum_{k=1}^{M}\left\lVert\mathbb{E}_{\pi|\theta_{k}}\left[\mathbbm{1}\left\lparen
m(\mathbf{x},\theta_{k})\leq
t\right\rparen\right]-\frac{1}{N_{\mathrm{test}}}\sum_{j=1}^{N_{\mathrm{test}}}\mathbbm{1}\left\lparen
m(x_{j};\theta_{k})\leq t\right\rparen\right\rVert_{\infty}.$ (30)

For ease of notation, replace the two terms in the difference of Equation (30)
by,

$\displaystyle A_{k}(t)$
$\displaystyle\doteq\mathbb{E}_{\pi|\theta_{k}}\left[\mathbbm{1}\left\lparen{m(\mathbf{x},\theta_{k}})\leq
t\right\rparen\right],\textrm{ and}$ (31) $\displaystyle B_{k}(t)$
$\displaystyle\doteq\frac{1}{N_{\mathrm{test}}}\sum_{j=1}^{N_{\mathrm{test}}}\mathbbm{1}\left\lparen
m(x_{j};\theta_{k})\leq t\right\rparen$ (32)

so that,

$\displaystyle\left\lVert\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}-\hat{G}\right\rVert_{\infty}\leq\frac{1}{M}\sum_{k=1}^{M}\left\lVert
A_{k}-B_{k}\right\rVert_{\infty}.$ (33)

For each $\theta_{k}$, $B_{k}$ is the eCDF of samples drawn from a
distribution with CDF $A_{k}$. By the Dvoretzky-Kiefer-Wolfwowitz (DKW)
inequality [9],

$\displaystyle\mathbb{P}_{\pi|\theta_{k}}\left\lparen\left\lVert
A_{k}-B_{k}\right\rVert_{\infty}>\delta_{b}\right\rparen\leq
2\exp(-2N_{\mathrm{test}}\delta_{b}^{2}).$ (34)

Taking a union bound over
$\mathcal{D}_{\mathrm{param}}=\\{\theta_{k}:k\in[M]\\}$ we get,

$\displaystyle\mathbb{P}_{\pi|\mathcal{D}_{\mathrm{param}}}\left\lparen\exists
k\in[M]\text{\ s.t.\ }\left\lVert
A_{k}-B_{k}\right\rVert_{\infty}>\delta_{b}\right\rparen\leq
2M\exp(-2N_{\mathrm{test}}\delta_{b}^{2}).$ (35)

Therefore, each term on the right-hand side of (33) is smaller than
$\delta_{b}$ w.h.p. and,

$\displaystyle\mathbb{P}_{\pi|\mathcal{D}_{\mathrm{param}}}\left\lparen\left\lVert\hat{F}_{\pi|\mathcal{D}_{\mathrm{param}}}-\hat{G}\right\rVert_{\infty}\leq\delta_{b}\right\rparen\geq
1-2M\exp(-2N_{\mathrm{test}}\delta_{b}^{2}),$ (36)

as desired.

###  7.2 Impartial trimming

In section 4, we introduce a robust version of the KS test, where we allow for
some outliers in samples from the null hypothesis using $\alpha$-trimming of
distributions. For completeness purposes, we present the concepts on
$\alpha$-trimming of a distribution w.r.t. a reference distribution as
introduced by Alvarez-Esteban et al. [1]. Trimming allows us to assume some
outliers in samples that did come from the null distribution and helps
quantify the fraction of these outliers. To elaborate on the concept of
trimming, let $k$ be the number of trimmed observations, and let $\alpha$ be
the trimming fraction, which implies $k\leq n\alpha$. Given n i.i.d samples
$\\{x_{i}\\}_{i=1}^{n}$ with probability distribution $P$, the empirical
measure can be defined as
$\frac{1}{n}\sum_{i=1}^{n}\delta^{\text{Dirac}}({x_{i}})$, where
$\delta^{\text{Dirac}}({x})$ is the Dirac measure. To not reject the null, we
can remove outliers by assigning a weight of 0 to the bad observations in the
sample and adjusting the weight of good observations by $\frac{1}{n-k}$.
However, we may not want to completely get rid of samples from the feasible
set and instead downplay the importance of bad observations by modifying the
empirical measure to be $\frac{1}{n}\sum_{i=1}^{n}w_{i}\delta_{x_{i}}$, where
$0\leq w_{i}\leq\frac{1}{(1-\alpha)}$, and $\frac{1}{n}\sum_{i=1}^{n}w_{i}=1$.
This is called an impartial trimming of the probability measure $P$. Using
impartial trimming interchangeably with trimming we look at the definition of
$\alpha$-trimming, denoted by $R_{\alpha}(P)$, of the distribution $P$ [1,
Definition 1],

$\mathcal{R}_{\alpha}(P)=\\{P^{*}\in\mathcal{P}:P^{*}\ll
P,\frac{dP}{dP^{*}}\leq\frac{1}{(1-\alpha)}P\text{ a.s}.\\}.$ (37)

where $P,P^{*}$ are probability measures on $\mathbb{R}$, $0\leq\alpha\leq 1$,
and we say $P^{*}$ is an $\alpha$-trimming of P if $P^{*}$ is absolutely
continuous w.r.t. $P$ and satisfies the above definition. The set of
$R_{\alpha}(P)$, the $\alpha$-trimmings of P, can be characterized in terms of
the trimming function h. The function h determines which zones in the
distribution P are downplayed or removed.

Let $\mathcal{C}_{\alpha}$ be the class of absolutely continuous functions
$h:[0,1]\rightarrow[0,1]$, such that $h(0)=0$, and $h(1)=1$, with derivative
$h^{\prime}$, such that $0\leq h^{\prime}\leq\frac{1}{1-\alpha}$. For any real
probability measure P, the following holds [1, Proposition 1 a.],

$\mathcal{R}_{\alpha}(P)=\\{P^{*}\in\mathcal{P}:P^{*}(-\infty,t]=h(P(-\infty,t]),h\in
C_{\alpha}\\}.$ (38)

###  7.3 Trimmed Kolmogorov-Smirnov distance

Given two distribution functions, $F_{\tau\times\pi}$ and $G_{k}$, the KS
distance is given by,

$d(F_{\tau\times\pi},G_{k})=\sup_{x\in\mathbb{R}}\lvert
F_{\tau\times\pi}(x)-G_{k}(x)\rvert.$ (39)

We analogously define the $\alpha$-trimmed KS distance functional as,

$d(F_{\tau\times\pi},R_{\alpha}(G_{k}))=\min_{\tilde{F}\in\mathcal{R}_{\alpha}(G_{k})}d(F_{\tau\times\pi},\tilde{F}).$
(40)

The plug-in estimator for $d(F_{\tau\times\pi},R_{\alpha}(G_{k}))$ is
$d(F_{\tau\times\pi},R_{\alpha}(\hat{G}_{k}))$, where $\hat{G}_{k}$ is the
empirical distribution function based on a sample of $N_{\mathrm{test}}$
independent random variables with common distribution function $G_{k}$, and
$F_{\tau\times\pi}$ is our null distribution. A practical computation for
$d(F_{\tau\times\pi},\mathcal{R}_{\alpha}(\hat{G}_{k}))$ uses function
$F_{\tau\times\pi}\circ G_{k}^{-1}$ to express
$d(F_{\tau\times\pi},R_{\alpha}(G_{k}))$. If $F_{\tau\times\pi}$ is continuous
and $x$ (generated samples) has distribution function $G_{k}$,
$F_{\tau\times\pi}\circ G_{k}^{-1}$ is the quantile function associated with
the random variable $Y=F_{\tau\times\pi}(x)$. This gives rise to the following
lemma and theorem for computing the trimmed KS distance between a theoretical
distribution and an eCDF [4],

######  Lemma 1 ([5], Lemma 2.4).

If $G_{k},F_{\tau\times\pi}$ are continuous distribution functions and $G_{k}$
is strictly increasing then,

$\displaystyle d(F_{\tau\times\pi},\mathcal{R}_{\alpha}(G_{k}))$
$\displaystyle=\min_{h\in C_{\alpha}}\lvert\lvert h-F_{\tau\times\pi}\circ
G_{k}^{-1}\rvert\rvert,\hskip
14.22636ptd(F_{\tau\times\pi},\mathcal{R}_{\alpha}(\hat{G}_{k}))$
$\displaystyle=\min_{h\in C_{\alpha}}\lvert\lvert
h-F_{\tau\times\pi}\circ\hat{G}_{k}^{-1}\rvert\rvert.$ (41)

######  Theorem 2 ([5], Theorem 2.5).

Suppose $\Gamma:[0,1]\rightarrow[0,1]$ is a continuous non-decreasing function
and let

$\displaystyle B(t)$ $\displaystyle=\Gamma(t)-\frac{t}{1-\alpha},$ (42)
$\displaystyle U(t)$ $\displaystyle=\sup_{t\leq s\leq 1}B(s),$ (43)
$\displaystyle L(t)$ $\displaystyle=\inf_{0\leq s\leq t}B(s),\textrm{ and}$
(44) $\displaystyle\tilde{h}_{\alpha}(t)$
$\displaystyle=\max\bigg{(}\bigg{(}\min\bigg{(}\frac{U(t)+L(t)}{2}\bigg{)},0\bigg{)}\frac{-\alpha}{(1-\alpha)}\bigg{)}.$
(45)

Then

$\displaystyle
h_{\alpha}\triangleq\tilde{h}_{\alpha}+\frac{(\cdot)}{(1-\alpha)}$ (46)

is an element of $C_{\alpha}$ and $\min_{h\in C_{\alpha}}\lvert\lvert
h-\Gamma\rvert\rvert=\lvert\lvert
h_{\alpha}-\Gamma\rvert\rvert=\lvert\lvert\tilde{h}_{\alpha}-\Gamma\rvert\rvert$,
with the assumptions on $\Gamma$ holding for $\Gamma=F_{\tau\times\pi}\circ
G_{k}^{-1}$.

In our application we don’t have access to $F_{\tau\times\pi}$, and as
discussed in Section 4, our robust hypothesis test is against the reference
function $\hat{G}$. Since this function is an average of eCDFs, which is not a
continuous function, for practical computation of
$d(\hat{G},\mathcal{R}_{\alpha}(\hat{G}_{k}))$, we consider the linearly
interpolated cdf of $\hat{G}$, which we denote as $\hat{G^{l}}$. Thus running
the hypothesis test in Equation $\eqref{eq: two samp robust threshold}$ is
equivalent to running the following test,

$\displaystyle
d(\hat{G^{l}},\mathcal{R}_{\alpha}(\hat{G}_{k}))\underset{\mathcal{\tilde{H}}_{0}}{\overset{\tilde{H}_{1}}{\gtrless}}\delta_{a}+\frac{1}{N_{\mathrm{test}}}$
(47)

or equivalently, from Equation (25),

$\displaystyle\max_{z}\min_{\tilde{F}\in\mathcal{R}_{\alpha}(\hat{G}_{k})}|\hat{G^{l}}(z)-\tilde{F}(z)|\underset{\mathcal{\tilde{H}}_{0}}{\overset{\tilde{H}_{1}}{\gtrless}}\delta_{a}+\frac{1}{N_{\mathrm{test}}},$
(48)

where we arrive at the threshold on the right-hand side by considering the
inequality
$\left\lVert\hat{G}-\hat{G^{l}}\right\rVert_{\infty}\leq\frac{1}{N_{\mathrm{test}}}$
and then applying the triangle inequality.

###  7.4 Additional Experimental Details

####  7.4.1 CNN architecture

Figure 15 shows the CNN architecture and parameters used in Section 5.

Figure 15: CNN Architecture

####  7.4.2 Observing the evolution of ensemble models over the training
procedure

Figure 16: 3-D scatter plot with two different views to visualize the
relationship of candidate models in $\mathsf{S}_{\mathrm{all}}$ at different
epochs with their ensemble at that epoch in terms of validation accuracy, ECE,
and churn w.r.t. ensemble (Z,Y,X axis) for View 1 (Left) and (Z,X,Y) axis for
View 2 (Right).

We once again focus on models in $\mathsf{S}_{\mathrm{all}}$, observing the
evolution of the ensemble of all candidate models in $M^{\prime}=100$ w.r.t.
its constituents at fixed training epochs. We track this evolution w.r.t. our
previously used metrics, validation accuracy, ECE, and churn w.r.t. the
ensemble. Figure 16 presents four snapshots of this evolution for candidate
models and their ensemble in $\mathsf{S}_{\mathrm{all}}$ at different epochs.
As candidate models evolve to become more accurate but also over-confident in
their predictions, their ensemble also evolves to become a better model: At
smaller epochs, poorly trained candidate models result in ensembles of poor
quality. This is evident for epoch 2, and epoch 10 in the top row of Figure 16
where the ensemble of candidate models is of poorer or comparable quality
w.r.t. to its constituents in terms of validation accuracy and ECE. As
candidate models are trained for longer durations, their ensembles are not
only better calibrated but also achieve the highest validation accuracy among
their constituents. We also note an overall reduction in the churn of
candidate models w.r.t. their ensemble as observed along the x-axis, with
increasing epochs, of View 1. This reduction can be attributed to the
improvement in the validation accuracy of candidate models over time. However,
churn between the ensemble model and its constituents that have achieved the
highest range of validation accuracy among its group increases from roughly a
range of 180-250 in epoch 10, to 250-400 in epoch 50. This can be explained by
the significant improvement in the validation accuracy of the ensemble
compared to its constituents at epoch 50. We also notice an increase in the
ECE values of candidate models that have achieved a high range of accuracy in
their group, as observed along the x-axis across different epochs in View 2.
This is observed as we move from epoch 10 to epoch 50, where the ECE values of
the aforementioned group of high-accuracy candidate models decreased roughly
from 0.01-0.02 in epoch 10 to 0.02-0.03 in epoch 50. This overconfidence in
candidate models that have achieved high accuracy in their group results in a
significant improvement in the quality of their ensemble. This non-inherent
trait of deep ensembles to be a better-calibrated counterpart of its
constituents under the commonly encountered scenario of the pool of candidate
models being overconfident was reported by Rahaman et al. [29].