Robust Nonparametric Hypothesis Testing to Understand Variability in
  Training Neural Networks
Training a deep neural network (DNN) often involves stochastic optimization,
which means each run will produce a different model. Several works suggest this
variability is negligible when models have the same performance, which in the
case of classification is test accuracy. However, models with similar test
accuracy may not be computing the same function. We propose a new measure of
closeness between classification models based on the output of the network
before thresholding. Our measure is based on a robust hypothesis-testing
framework and can be adapted to other quantities derived from trained models.
Sinjini Banerjee, Reilly Cannon, Tim Marrinan, Tony Chiang, Anand D. Sarwate
# Robust Nonparametric Hypothesis Testing to Understand Variability in Training Neural Networks

###### Abstract

Training a deep neural network (DNN) often involves stochastic optimization, which means each run will produce a different model. Several works suggest this variability is negligible when models have the same performance, which in the case of classification is test accuracy. However, models with similar test accuracy may not be computing the same function. We propose a new measure of closeness between classification models based on the output of the network before thresholding. Our measure is based on a robust hypothesis-testing framework and can be adapted to other quantities derived from trained models.

Sinjini Banerjee\({}^{\star}\)  Reilly Cannon\({}^{\dagger}\)   Tim Marrinan\({}^{\dagger}\)   Tony Chiang\({}^{\dagger}\)   Anand D. Sarwate\({}^{\star}\)\({}^{\star}\)Rutgers University  \({}^{\dagger}\)Pacific Northwest National Lab

+
Footnote †: The work of S.B. and A.D.S. were supported in part by a Pacific Northwest National Laboratory Program (PNNL) under contract DR00022921. R.C., T.M., and T.C. were partially supported by the Mathematics for Artificial Reasoning in Science (MARS) initiative via the Laboratory Directed Research and Development (LDRD) at PNNL. T.C. and A.D.S. were also partially supported by the Statistical Inference Generates kNowledge for Artificial Learners (SIGNAL) program at PNNL.

DNN variability, Nonparametric hypothesis testing, Robust Kolmogorov-Smirnov test

## 1 Introduction

State-of-the-art deep learning models have been remarkably successful in achieving state-of-the-art performance on complex tasks in the fields of healthcare, education, cyber-security, and other important domains. Training these models take significant time, energy, and financial resources. Because training is done with stochastic algorithms for nonconvex optimization, models produced by different training runs in general converge to different solutions. These solutions can be very different, even if they have a similar objective value and test loss. Two models may differ in their predictions on the same test point: this phenomenon is known as churn [1] and is a cause for model irreproducibility.

Deep learning models are 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 certain initialization or hyperparameter 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 if one is better than another becomes a difficult task.

Reproducibility in training is an active area of recent work. Gundersen et al. [2] identified model initialization, random batch shuffling during mini-batch stochastic gradient descent (SGD), data sampling, and parallel execution over GPUs as some of the major sources of randomness present in the training procedure. Fort et al. [3], and Bouthillier et al. [4] showed random data ordering has a smaller effect than random initialization on model performance. Sompealli et al. [5] used decision boundaries to characterize the reproducibility of models. Summers and Dinneen [6], and Jordan [7] argue that variance matters only during initial conditions of the training procedure. Most work assesses the impact randomness has on the _test accuracy_ of models or _churn_ between models, both of which focus on the _decisions_ made by predictive (classification) models. Closeness of decisions does not imply closeness of the trained models and high test accuracy does not imply the stability of the learned features or their meaningful contribution to class differentiation [8].

In this paper we try to quantify the variability caused by training in terms of the _network outputs used to make the decision_. Figure 1 illustrates the difference: the sold and dashed lines represent two decision boundaries between red (circle) and blue (cross). Test accuracy measures incorrect decisions and the churn is given by the region between the two curves. We are interested in the shading, which shows the network output which is thresholded 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, instead of looking at whether models make the same number of correct decisions, we can examine the distribution of "confidence" values (the _logit gap_) from functions learned by different runs of a fixed deep neural network architecture. In particular, we use measures derived from nonparametric hypothesis testing [9]. We can frame model similarity as asking whether two models generate similar distributions of logit gaps or, with an appro

Figure 1: Illustration of two models that have the same test accuracy but different decision boundaries.
priately defined null, use a goodness-of-fit test to assess whether a particular model is close to the null. Nonparametric goodness-of-fit tests are often sensitive in the large sample regime. To remedy this we use concepts derived from robust statistics [10].

## 2 Problem Setup

We focus on a binary classification problem1 in which data is generated from an (unknown) probability distribution \(\pi\) on the space of feature-label pairs \(\mathcal{X}\times\mathcal{Y}=\mathbb{R}^{d}\times\{0,1\}\). We interpret a predictive model (e.g. a neural network) as a function \(m\colon\mathcal{X}\to\mathbb{R}\), where we interpret \(m(x)\) as a log posterior probability ratio \(m(x)=\log\frac{\mathbb{P}(y=1|x)}{\mathbb{P}(y=0)x}\).The _prediction_ produced by the model is \(\hat{y}(x)=\mathbb{1}(m(x)\geq 0)\). We can also think of the model as computing a pair of outputs \(m^{+}(x)\) and \(m^{-}(x)\) where the posteriors are given by the softmax function: \(\mathbb{P}(y=1\mid x)=\frac{\exp(m^{+}(x))}{\exp(m^{+}(x))+\exp(m^{-}(x))}\). This makes \(m(x)=m^{+}(x)-m^{-}(x)\), so we refer to \(m^{+}(x)\) and \(m^{-}(x)\) as _logits_ and \(m(x)\) as the _logit gap_.

Footnote 1: While we focus here on binary classification, our framework can be extended to more general prediction problems.

A neural network with a given architecture is parameterized by a space of parameters (e.g. the weights) \(\Theta\). Thus for every parameter setting \(\theta\in\Theta\) the network computes a function \(m(x;\theta)\) and the neural network defines a family of functions \(\mathcal{M}=\{m_{\theta}:\mathcal{X}\to\mathcal{Y}\colon\theta\in\Theta\}\). Let \([n]=\{1,2,\ldots,n\}\). A training algorithm takes a _training set_\(\mathcal{D}_{\mathrm{train}}=\{(x_{i},y_{i}):i\in[N]\}\stackrel{{ \mathrm{i.i.d.}}}{{\sim}}\pi\) and "learns" a parameter setting \(\theta\) by approximately minimizing an empirical risk \(\hat{R}(\theta;\mathcal{D}_{\mathrm{train}})\) computed over the training data. The output of the training algorithm is a model \(m(x;\theta)\) with predictions \(\hat{y}(x;\theta)=\mathbb{1}(m(x;\theta)\geq 0)\).

The optimization algorithms used in NN training are approximate in two ways. First, the risk minimization problem is in general nonconvex so in general they will converge to a local minimum. Second, they are usually stochastic, as in SGD, which means the estimated parameters \(\theta\) are themselves random variables. We can think of an NN training algorithm as sampling from a distribution on \(\Theta\) and therefore sampling from the space of functions \(\mathcal{M}\).

Since two runs of the training algorithm can produce different functions, it is natural 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^{\prime}]\}\) which we assume is also sampled i.i.d. from \(\pi\). The _test accuracy_ of a model is

\[A(\theta)=\frac{1}{N^{\prime}}\sum_{i=1}^{N^{\prime}}\mathbb{1}(\hat{y}(x_{j} ;\theta)=y_{j}). \tag{1}\]

Two models \(m(\cdot;\theta_{1})\) and \(m(\cdot;\theta_{2})\) which have similar test accuracy make the same number of mistakes. The _churn_ is defined by

\[C(\theta_{1},\theta_{2})=\frac{1}{N^{\prime}}\sum_{i=1}^{N^{\prime}}\mathbb{1} \left(\hat{y}(x_{j};\theta_{1})\neq\hat{y}(x_{j};\theta_{2})\right), \tag{2}\]

which is the fraction of training points where the models disagree. Two models with low churn make almost the same decisions (regardless of whether they are correct or not).

Both test accuracy and churn focus on the _predictions_ made by models and do not use information about the logit gap function \(m(x)\) beyond its sign. Looking at \(m(x)\) 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. Given the test set, our central object of inquiry is the empirical distribution of logit gaps, which we can write as a cumulative distribution function (CDF):

\[F(x;\theta)=\frac{1}{N^{\prime}}\sum_{j=1}^{N^{\prime}}\mathbb{1}(m(x_{j}; \theta)\leq x). \tag{3}\]

With some abuse of notation, we consider \(\mathcal{P}\) to be the set of probability measures on \(\mathbb{R}\) and a \(F\) as an element of \(\mathcal{P}\).

We consider the following statistical setup. The training algorithm takes the training set \(\mathcal{D}_{\mathrm{train}}\) and uses randomization in several different ways to produce a model. Each run of the training algorithm uses independent randomness: training is therefore sampling parameters \(\theta_{1},\theta_{2},\ldots,\theta_{M}\) i.i.d. from an (unknown) distribution on \(\Theta\) induced by the training algorithm. These parameters correspond to \(M\) i.i.d. samples \(\{m_{k}(x)=m(x;\theta_{k}):k\in[M]\}\) taking values in \(\mathcal{M}\). When applied to the test set, these produce \(M\) CDFs \(\{F_{k}(x)=m(x;\theta_{k}):k\in[M]\}\).

We can control which sources of randomness are used in the training by altering the training procedure. For example, we can use deterministic initialization or fixed batch ordering. Under these scenarios, we can generate \(M\) models and ask if the models are more different from each other when we turn on or off the sources of randomness. To measure how different the models are, we use a framework from nonparametric hypothesis testing.

## 3 Goodness-of-fit Testing for Logit gaps

We can define an expected logit gap function \(\bar{m}(x)\) by integrating over the distribution on \(\mathcal{M}\) induced by the training algorithm. Given a test set \(\mathcal{D}_{\mathrm{test}}\), this produces an expected CDF which we will call \(F_{0}\). If we knew what \(F_{0}\) was, we could assess whether or not a given model \(m(x)\) is a representative sample from the training algorithm by testing if the model's CDF \(F(x)\) is close to the expected \(F_{0}(x)\). This corresponds to the following one-sided hypothesis test:

\[\mathcal{H}_{0}\colon\{m(x_{j})\colon j\in[N^{\prime}]\}\sim F _{0} \tag{4}\] \[\mathcal{H}_{1}\colon\{m(x_{j})\colon j\in[N^{\prime}]\}\not<F_{0} \tag{5}\]

This is a classical nonparametric goodness-of-fit testing problem which can be solved using the Kolmogorov-Smirnov (KS) test [11, 12, 13]: \(d_{k}(F_{0},F)=\|F_{0}(x)-F(x)\|_{\infty}\stackrel{{\mathcal{H}_{ 1}}}{{\sim}}\tau\), where \(F(x)\) is the empirical CDF from \(m(x)\) and we can set the threshold \(\tau\) to achieve the desired error tradeoff.

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. Note that we do not expect our empirical samples to look exactly like they were drawn from \(F_{0}\). Ideally, we want a test which allows for some outliers. Furthermore, we do not know what the null hypothesis \(F_{0}\) is, so even computing the test statistic is not possible. To address these issues, we use ideas from _robust statistics_[10]. Given a distribution \(P\in\mathcal{P}\) and \(\alpha\in[0,1]\), we can define the set of \(\alpha\)-contaminated distributions as \(\mathcal{R}_{\alpha}(P)=\{(1-\alpha)P+\alpha Q\colon Q\in\mathcal{P}\}\). Alvarez-Esteban et al. [14] interpret \(\mathcal{R}_{\alpha}(P)\) as the set of so-called \(\alpha\)-trimmings of \(P\):

\[\mathcal{R}_{\alpha}(P)=\bigg{\{}Q\in\mathcal{P}\colon Q\ll P,\frac{dP}{dQ} \leq\frac{1}{(1-\alpha)}P-\text{ a.s.}\bigg{\}}. \tag{6}\]

The advantage of this definition, as shown by del Barrio et al. [15], is that we can compute, for a given \(F\), the closest \(L_{\infty}\) approximation
to \(F_{0}\) in the set of \(\alpha\)-trimmings of \(F\).

\[d_{k}(F_{0},\mathcal{R}_{\alpha}(F))=\min_{\tilde{F}\in\mathcal{R}_{\alpha}(F)} \lVert F_{0}-\tilde{F}\rVert_{\infty}. \tag{7}\]

The core idea behind trimming is to ask if trimming a small fraction of samples (from the test set) would allow a KS test to accept the null hypothesis. The set \(\mathcal{R}_{\alpha}(F)\) of \(\alpha\) trimmings of P can be characterized in terms of a "trimming function" \(h\) and we can efficiently find the optimizing \(\tilde{F}\in\mathcal{R}_{\alpha}(F)\) in (7) by optimizing over these trimming functions. This technique led to a recent work of del Barrio et al. [15] that proposes a _robust KS test_ based on trimming.

## 4 Using Trimming to Estimate Variability

Our main contribution is a new way to measure the dissimilarity of trained models. To do this we first develop a two-sample version of the robust KS test discussed in the previous section.

**Using a deep ensemble as a base model.** While we do not know the expected CDF \(F_{0}\), note that for a given test point \((x_{j},y_{j})\), we can use the empirical mean \(\frac{1}{M}\sum_{k=1}^{M}m_{k}(x)\) as an estimate of expected logit gap function \(\bar{m}(x)=\mathbb{E}[m(x)]\) at test point \(x\). To evaluate whether a particular model with CDF \(F_{\ell}\) is close to \(F_{0}\) we can instead measure its closeness to the _leave-one-out ensemble_ (LOOE):

\[m_{-\ell}(x) =\frac{1}{M-1}\sum_{k\neq\ell}m_{k}(x) \tag{8}\] \[F_{-\ell}(x) =\frac{1}{N^{\prime}}\sum_{j=1}^{N^{\prime}}\mathbb{1}(m_{-\ell }(x_{j})\leq x). \tag{9}\]

The LOOE corresponds to a deep ensemble predictor [16, 17, 18], which has been used to reduce variability in DNN models. Taking the average of model "confidences" across independent training runs make them closer to their expected values, lowering variability. Figure 2 shows how the logit gap samples obtained from averaging over all candidate models compare with candidate models in the pool. The ensemble model produces a lower number of samples with small logit gaps (samples with higher uncertainty), and large logit gaps (overconfident samples).

**Proposed algorithm.** In our robust two-sample KS test, we are comparing two empirical CDFs, a candidate \(F_{\ell}\) and the LOOE \(F_{-\ell}\), and asking if they were generated by a common underlying distribution:

\[\mathcal{H}_{0}\colon F_{\ell},F_{-\ell}\text{ are from the same distribution} \tag{10}\] \[\mathcal{H}_{1}\colon F_{\ell},F_{-\ell}\text{ are not from the same distribution} \tag{11}\]

We use bootstrap sampling to resample \(N^{\prime}\) independent test points to compute \(F_{\ell}\) and \(F_{-\ell}\). For a fixed \(\alpha\), we calculate the \(\alpha\)-trimming of \(F_{\ell}\) that minimizes the \(L_{\infty}\) distance to \(F_{-\ell}\):

\[d_{k}(F_{-\ell},\mathcal{R}_{\alpha}(F_{\ell}))=\min_{F\in\mathcal{R}_{\alpha }(F_{\ell})}\lVert F_{-\ell}-F\rVert_{\infty}. \tag{12}\]

Our hypothesis test is then

\[d_{k}(F_{-\ell},\mathcal{R}_{\alpha}(F_{\ell}))\underset{\mathcal{R}_{0}}{ \overset{\mathcal{H}_{1}}{\gtrless}}\tau. \tag{13}\]

To set the threshold \(\tau\), we observe that if \(\tilde{F}\) is the minimizer in (12), the two-sample version of the Dvoretzky-Kiefer-Wolfowitz inequality due to Wei and Dudley [19] implies that

\[\mathbb{P}\bigg{(}\sup_{x}\Bigl{|}\tilde{F}(x)-F_{-\ell}(x)\Bigr{|}>\tau\bigg{)} \leq Ce^{-2N^{\prime}\tau^{2}}. \tag{14}\]

For \(N^{\prime}>458\) the leading constant \(C\) can be replaced with \(2\). In that case, to set the probability of falsely rejecting \(\mathcal{H}_{0}\) to \(\delta\), we can set the right hand side equal to \(\delta\) and solve for \(\tau=\sqrt{N^{\prime}\ln\frac{2}{\delta}}\). We use this test in Algorithm 1 to compute the largest \(\alpha\) for which the test fails to reject \(\mathcal{H}_{0}\). The output \(\hat{\alpha}\) is our proposed measure of discrepancy between model \(m_{\ell}\) and the LOOE model.

```
Data:\(\mathcal{D}_{\mathrm{test}}\) with \(N^{\prime}>458\) samples, trained models \(\{m_{k}\colon k\in[M]\}\), model index \(\ell\in[M]\), threshold \(\tau\), resampling number \(B\), trimming levels \((\alpha_{1},\ldots,\alpha_{T})\) Result: trimming level estimate \(\hat{\alpha}\) for\(b=1\) to \(B\)do  Compute \(F_{\ell}\) using \(N^{\prime}\) resampled from \(\mathcal{D}_{\mathrm{test}}\);  Compute \(F_{-\ell}\) in (9) using \(N^{\prime}\) resampled from \(\mathcal{D}_{\mathrm{test}}\); \(\mathrm{Reject}\gets 1\), \(\hat{\alpha}_{b}\gets 0\), \(t\gets 1\); while\(\mathrm{Reject}=1\)and \(t\leq T\)do \(\alpha\leftarrow\alpha_{t}\); \(\mathcal{H}\leftarrow\) output of (13) ; if\(\mathcal{H}=\mathcal{H}_{0}\)then  Reject \(\gets 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

## 5 Experiments

To illustrate our approach, we used a small convolutional neural network on two classes from the CIFAR-10 dataset [20]: \(N=12000\) and \(N^{\prime}=4000\). Our network had 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. The small size allows us to train many more models to explore model training in different scenarios: \(\mathsf{S}_{\mathrm{init}}\) with only random initialization, \(\mathsf{S}_{\mathrm{batch}}\) with only random batch selection in SGD, \(\mathsf{S}_{\mathrm{train}}\) with only randomly resampled training data, and \(\mathsf{S}_{\mathrm{all}}\) with all sources of randomness. For each scenario we trained \(M=100\) models for \(50\) epochs.

**Similar accuracy/churn do not imply small \(\hat{\alpha}\).** To show that our measure \(\hat{\alpha}\) provides a different understanding of model closeness, Figure 2 shows the logit gap distribution for two models with respect to their LOOE. The churn is computed w.r.t to the LOOE. In general, if models need a small trimming level to be accepted under the null hypothesis they are likely also to have a low churn and high test accuracy. However, the opposite is not necessarily always true. As seen in Figure 2, Model 2 has produced more uncertain samples with a higher probability when compared to another candidate model with a similar test accuracy and churn. This could make Model 2 more prone to adversarial attacks. Based on this we observe that a high test accuracy or a low churn is not enough to conclude that one model is a better model than all other candidate models belonging to the same pool. Thus, the recommendation is to choose models that not only have achieved good test accuracy and low churn but are also better representatives of the LOOE. Relying on \(\hat{\alpha}\) avoids models that somehow "got lucky" and lets us choose better models that have performed meaningfully well.

**Comparing the evolution of \(\hat{\alpha}\) and test accuracy.** To understand the difference between \(\hat{\alpha}\), test accuracy, and churn, we considered the baseline scenario \(\mathsf{S}_{\mathrm{all}}\). We trained \(M=100\) models for \(50\)
epochs to examine the evolution of each model over those 50 epochs compared to the LOOE at 50 epochs. We have only included test accuracy plots due to space restrictions but the churn of these models w.r.t to their LOOE also follows a similar behavior. Figure 3 shows four snapshots of the relationship between the test accuracy of a model and \(\hat{\alpha}\).

At smaller epochs, although test accuracy improves, \(\hat{\alpha}\) remains large for most models. After training the models long enough, we stop seeing a large improvement in test accuracy but candidate models move closer in distribution toward their LOOE model at a higher epoch. This suggests that candidate models evolve to become better representatives of the training algorithm long after reducing test accuracy variability, and \(\hat{\alpha}\) is indicative of this evolution.

**Using \(\hat{\alpha}\) to examine the impact of randomization.** Our focus is on two epochs, one very early on in the training (epoch 2) and another after training for a long period of time (epoch 50). At epoch 2, models differing according to \(\mathsf{S}_{\mathrm{init}}\) had a higher percentage of models that do not reject the null hypothesis, with a low \(\hat{\alpha}\), than \(\mathsf{S}_{\mathrm{train}}\), or \(\mathsf{S}_{\mathrm{batch}}\). This would suggest that \(\mathsf{S}_{\mathrm{init}}\), contributes the least among individual sources of variance to \(\mathsf{S}_{\mathrm{all}}\), by a significantly large margin. However, both \(\mathsf{S}_{\mathrm{init}}\), and \(\mathsf{S}_{\mathrm{batch}}\) seem to contribute less than \(\mathsf{S}_{\mathrm{train}}\) in test accuracy variability at epoch 2. This sensitivity of \(\hat{\alpha}\) to model similarity/dissimilarity is lost at higher epochs. At epoch 50, we see a reduction in variability both in test accuracy and \(\hat{\alpha}\). This points us to the work done by [7], that individual sources stop contributing independently to the total variance observed from changing all sources of randomness if you train the models long enough.

## 6 Conclusion and Future Work

In this work, we highlighted a different approach to analyzing variability in deep neural networks. Our proposed framework is based on a robust two-sample hypothesis testing problem that uses impartial trimming of the empirical CDF of samples obtained from the logit gap function. The purpose of this test is to assess the similarity/dissimilarity between candidate models in a pool with their consensus (leave-one-out ensemble) model 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 over the commonly used test accuracy.

While in this paper we describe the methodology and an example on a small model, future extensions of this work include applications to very large deep net models and extensions to multi-class classification. Another direction is to explore two-sample hypothesis testing based on other distance metrics like the Wasserstein metric. In this work, we chose to focus 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 eigen-distribution of the Jacobian or the Neural Tangent Kernel of functions learned by these models. The eigen-distribution of these matrices can provide us information on how well individual candidate models can generalize to test data.

Figure 4: Bar plots showing variation in test accuracy and alpha, for Epoch 2 (top row), and Epoch 50 (bottom row).

Figure 3: Scatterplot of \(\hat{\alpha}\) vs. test accuracy for different epochs under \(\mathsf{S}_{\mathrm{all}}\). Each model is compared to the LOOE model at epoch 50.

Figure 2: (Left) Histogram of logit gap from the ensemble model with the upper and lower envelope representing the maximum and minimum probability attained in each bin among individual candidate models. (Right) Histogram plots of logit gap samples from two candidate models and their LOOE models at a fixed epoch.