A Comparison of the Delta Method and the Bootstrap in Deep Learning Classification

In the context of deep learning classification, the classical Bootstrap method starts by creating $B$ datasets from the original dataset by sampling with replacement. Subsequently, $B$ networks are trained separately on each of the bootstrapped datasets. The epistemic uncertainty for each of the $T_{L}$ class predictions (in standard deviations) associated with prediction of $x_{0}$ is obtained by the sample standard deviation over the ensemble of $B$ predictions,

$$\widetilde{\sigma}_{\text{boot}}(x_{0})=\sqrt{\frac{1}{B-1}\sum_{b=1}^{B}(\hat{y}_{0}^{(b)}-\overline{\hat{y}_{0}})^{2}}\in\mathbb{R}^{T_{L}},$$

(1)

where the vector $\hat{y}_{0}^{(b)}$ represents the $T_{L}$ predictions for $x_{0}$ (one probability per class) obtained from the $b$th bootstrapped network, and where $\overline{\hat{y}_{0}}$ is the sample mean,

$$\overline{\hat{y}_{0}}=\frac{1}{B}\sum_{b=1}^{B}\hat{y}_{0}^{(b)}\in\mathbb{R}^{T_{L}}.$$

(2)

The method is easy to implement efficiently in practice. Training $B$ networks is an ‘embarrassingly’ parallel problem, and the space complexity for the bootstrapped datasets is just $O(BN)$ when an indexing scheme is used for the sampling with replacement. The experiments conducted in this paper is based on the example pydeepboot.py provided in the pydeepdelta provision [14].

The Delta method was adapted to the deep learning classification context by [11]. The adaption addresses several fundamental difficulties that arise when the method is applied in deep learning. In essence, it is shown that an approximation of the eigendecomposition of the Fisher information matrix utilizing only $K$ eigenpairs allows for an efficient implementation with bounded worst-case approximation errors. We briefly review the standard method here for convenience.

An approximation of the epistemic component of the uncertainty associated with the prediction of $x_{0}$ can be found by the formula

$$\widetilde{\sigma}_{\text{delta}}(x_{0})=\sqrt{\text{diag}\big{(}F\Sigma F^{T}\big{)}}\in\mathbb{R}^{T_{L}},$$

(3)

where the sensitivity matrix $F$ in (3) is defined

$$F=\begin{bmatrix}F_{ij}\end{bmatrix}\in\mathbb{R}^{T_{L}\times P},~{}F_{ij}=\frac{\partial}{\partial\omega_{j}}f_{i}(x_{0},\omega)\bigg{\rvert}_{\omega=\hat{\omega}}.$$

(4)

The covariance matrix $\Sigma$ in (3) can be estimated by several alternative estimators. In [11] it was demonstrated that the Hessian estimator, the Outer-Products of Gradients (OPG) estimator and the Sandwich estimator lead to nearly perfect correlated results for two different deep learning models. Since the models discussed in this paper are identical to those in [11], we thus focus only on one of the estimators, namely the OPG estimator defined by

$$\Sigma=\frac{1}{N}G^{-1}=\frac{1}{N}\left[\frac{1}{N}\sum_{n=1}^{N}\frac{\partial C_{n}}{\partial\omega}\frac{\partial C_{n}}{\partial\omega}^{T}\bigg{\rvert}_{\omega=\hat{\omega}}+\lambda I\right]^{-1}\in\mathbb{R}^{P\times P},$$

(5)

where the summation part of $G$ corresponds to the empirical covariance of the gradients of the cost function evaluated at $\hat{\omega}$. As discussed in [11], the term $\lambda I$ is explicitly added in order to make the OPG estimator asymptotically equal to the Hessian estimator, as is the primary motivation for the former as a plug-in replacement of the latter in the first place.

When the Delta method is implemented under the framework of [11], it has several desirable properties: a) requires only $O(PK)$ space and $O(KPN)$ time, b) fits well with deep learning software frameworks based on automatic differentiation, c) works with any $L_{2}$-regularized neural network architecture, and d) does not interfere with the training process as long as the norm of the gradient of the cost function is approximately equal to zero after training.

There are $L=6$ layers, layer $l=1$ is the input layer represented by the input vector. Layer $l=2$ is a $3\times 3\times 1\times 32$ convolutional layer followed by max pooling with stride equal to $2$ and with a ReLU activation function. Layer $l=3$ is a $3\times 3\times 32\times 64$ convolutional layer followed by max pooling with a stride equal to $2$, and with ReLU activation function. Layer $l=4$ is a $3\times 3\times 64\times 64$ convolutional layer with ReLU activation function. Layer $l=5$ is a $576\times 64$ dense layer with ReLU activation function, and the output layer $l=6$ is a $64\times T_{L}$ dense layer with softmax activation function, where the number of classes (outputs) is $T_{L}=10$. The total number of parameters is $P=93322$.

There are $L=6$ layers, layer $l=1$ is the input layer represented by the input vector. Layer $l=2$ is a $3\times 3\times 3\times 32$ convolutional layer followed by max pooling with stride equal to $2$ and with a ReLU activation function. Layer $l=3$ is a $3\times 3\times 32\times 64$ convolutional layer followed by max pooling with a stride equal to $2$, and with ReLU activation function. Layer $l=4$ is a $3\times 3\times 64\times 64$ convolutional layer with ReLU activation function. Layer $l=5$ is a $1024\times 64$ dense layer with ReLU activation function, and the output layer $l=6$ is a $64\times 10$ dense layer with softmax activation function, where the number of classes (outputs) is $T_{L}=10$. The total number of parameters is $P=122570$.

For the Bootstrap networks, we test two different weight initialization variants: dynamic random normal weight initialization (DRWI) and static random normal weight initialization (SRWI). The former uses a different (e.g. dynamic) seed across the replicates, meaning that each network in the DRWI Bootstrap ensemble will start out with different random weight values. The latter case uses the same (e.g. static) seed across the replicates, and hence all the networks in the SRWI Bootstrap ensemble receives the same random initial weight values. For all networks, we use zero bias initialization. Futhermore, to investigate the impact of random weight initialization on the Delta method, we apply the Delta method 16 times on a set of 16 networks distinguished only by DRWI.

We use the cross-entropy cost function with a $L_{2}$-regularization rate $\lambda=0.01$, and utilize the Adam [7, 1] optimizer with a batch size of $100$, and no form of randomized data shuffling. To ensure convergence (e.g. $||\nabla C(\hat{\omega})||_{2}\approx 0$), we apply two slightly different learning rate schedules given by the following (step, rate) pairs: MNIST = $\{(0,10^{-3}),(60\text{k},10^{-4}),(70\text{k},10^{-5}),(80\text{k},10^{-6})\}$ and CIFAR-10 = $\{(0,10^{-3}),(55\text{k},10^{-4}),(85\text{k},10^{-5}),(95\text{k},10^{-6},(105\text{k},10^{-7})\}$. For MNIST, we stop the trainings after $90,000$ steps, while for CIFAR-10, after $115,000$ steps – corresponding to the overall training statistics shown in Table 1.

The results from the full set of regressions ($d=1,2,\ldots,16$) holding a fixed $B=100$ are shown in Figure 3. The primary observations are as follows: The mean squared correlation coefficients $R^{2}$ are generally high for MNIST and CIFAR-10, meaning that there is a strong linear relationship between the uncertainty levels obtained by the Bootstrap and the Delta method. For the lowest $K$, the $R^{2}$ starts out at 90% for MNIST, and at 81% for CIFAR-10. As $K$ grows, an increase by only $4$% is observed for MNIST, while 8% for CIFAR-10. The major difference observed as $K$ increases lies in the absolute uncertainty levels expressed by the slope $\beta$: for MNIST, the slope stabilizes at around $K=600$ while at about $K=1000$ for CIFAR-10. The same trend is reflected in the maximum approximation errors $\epsilon$, where we respectively see them approach zero at the same values for $K$. Although not shown in the plots, the regression intercepts $\alpha$ are always zero, meaning that there is no offset in the uncertainty estimates by the two methods.

The main difference found from applying DRWI opposed to SRWI for the Bootstrap ensembles, is that the absolute level of uncertainty increases with DRWI. This is expected, since the DRWI version of the Bootstrap will be more prone to reaching different local minima, and therefore also captures this additional variance. Supporting evidence for this hypothesis is evident by CIFAR-10’s wider confidence intervals. A more pronounced geometry difference across various local minima will ultimately lead to higher variability in the $R^{2}$ and $\beta$. A slightly higher mean $R^{2}$ (+1-2%) is also observed for the DRWI version of the Bootstrap. This is reasonable given the fact that also the Delta method networks are more prone to reaching different local minima across the 16 repetitions because of DRWI.

Figure 4 shows the same type of comparison when the number of Bootstrap replicates $B$ varies, and the number of eigenpairs are fixed ($K=1500$ for MNIST and $K=2500$ for CIFAR-10). The main observation from this experiment is that there is very little to achieve by selecting a larger ensemble size $B$ than about 50, as this is the point where the mean slope and squared correlation coefficient stabilizes.

Table 2 shows the computation time for the two methods when executed on a Nvidia RTX 2080 Ti based GPU. For MNIST, the smallest $K$ leading to acceptable approximation errors and stable absolute uncertainty levels for the Delta method is at $K=600$, while for CIFAR-10 the same applies at $K=1000$. Furthermore, the smallest acceptable $B$ leading to stable correlation and absolute uncertainty levels for the Bootstrap is at $B=50$. We conclude that in these experiments the Delta method outperforms the Bootstrap in terms of computation time by a factor $4.6$ on MNIST, and a factor $5.9$ for CIFAR-10.