Uncertainty quantification in neural network classifiers -- a local
  linear approach
Classifiers based on neural networks (NN) often lack a measure of uncertainty
in the predicted class. We propose a method to estimate the probability mass
function (PMF) of the different classes, as well as the covariance of the
estimated PMF. First, a local linear approach is used during the training phase
to recursively compute the covariance of the parameters in the NN. Secondly, in
the classification phase another local linear approach is used to propagate the
covariance of the learned NN parameters to the uncertainty in the output of the
last layer of the NN. This allows for an efficient Monte Carlo (MC) approach
for: (i) estimating the PMF; (ii) calculating the covariance of the estimated
PMF; and (iii) proper risk assessment and fusion of multiple classifiers. Two
classical image classification tasks, i.e., MNIST, and CFAR10, are used to
demonstrate the efficiency the proposed method.
Magnus Malmström, Isaac Skog, Daniel Axehill, Fredrik Gustafsson
# Uncertainty quantification in neural network classifiers - a local linear approach

###### Abstract

Classifiers based on neural networks (nn) often lack a measure of uncertainty in the predicted class. We propose a method to estimate the probability mass function (pmf) of the different classes, as well as the covariance of the estimated pmf. First, a local linear approach is used during the training phase to recursively compute the covariance of the parameters in the nn. Secondly, in the classification phase another local linear approach is used to propagate the covariance of the learned nn parameters to the uncertainty in the output of the last layer of the nn. This allows for an efficient Monte Carlo (mc) approach for: (i) estimating the pmf; (ii) calculating the covariance of the estimated pmf; and (iii) proper risk assessment and fusion of multiple classifiers. Two classical image classification tasks, i.e., mnist, and cfar10, are used to demonstrate the efficiency the proposed method.

Neural networks; Uncertainty descriptions; Information and sensor fusion; Identification and model reduction; Intelligent driver aids; Nonlinear system identification; +
Footnote †: footnote]This paper was not presented at any IFAC meeting. Corresponding author M. Malmström.

## 1 Introduction

In this paper, the problem of quantifying the uncertainty in the predictions from a neural network (nn) is studied. The uncertainty in the prediction stems from three different sources: errors caused by the optimization algorithm that is used to train the nn, errors in the data (aleatoric uncertainty), and errors in the model (epistemic uncertainty). In this paper, the focus is on uncertainty from the two latter sources.

In numerous applications, e.g., image recognition [1], learning properties in atoms [2], and various control tasks [3, 4], nns have shown high performance. Despite their high performance, the use of nns in safety-critical applications is limited [5, 6, 7]. It is partly a consequence of the fact that their predictions usually do not come with any measure of certainty of the prediction, which is crucial to have in a decision-making process in order to know to what degree the prediction can be trusted. Moreover, the quantified measure of uncertainty can be used to detect and remove outliers in the data. Furthermore, it is not possible to fuse the prediction from the nn with information from other sensors without knowledge about the uncertainty.

Autonomous driving is an example of a safety-critical application in which it is relevant to be able to perform reliable classifications of, e.g., surrounding objects. In particular, this need was highlighted in the fatal Uber accident in 2018 where the lack of reliable classifications of surrounding objects played a role in the development of events that eventually led to the accident [8].

The problem to quantify the uncertainty in the prediction of nns has lately gained increasing attention, and numerous methods to calculate the uncertainty have been suggested [9, 10, 11, 12, 13]. For a survey of methods see [14]. The methods suggested in the literature can broadly be divided into one out of two categories. One category is based on creating an ensemble of predictions from which the uncertainty in the prediction is computed [15, 16, 17, 18, 19, 20, 21, 22]. In the other category, the nn structure is extended and the nn is trained to learn its own uncertainty [23, 24, 25, 26, 27, 28].

Concerning the first category, it has for example been suggested to create an ensemble by training multiple
nns, from whose predictions the uncertainty is computed by [15]. Since training a single nn is often computationally expensive, this method has high computational complexity. In practice, it is only feasible from a computational perspective to create small ensembles. To decrease the computational complexity, it was in [16, 17] suggested to use already existing regularization techniques (dropout and batch norm) to sample values of the parameters of the nn from which these ensembles can be created. Another method to create ensembles is by sampling values of the parameters during the last part of the training phase [18, 19]. So-called test-time data augmentation methods have also been suggested to do perturbation on the test data to create an ensemble of predictions [20]. Even though the methods in [16, 17, 18, 19, 20] do not need multiple models to be trained they require multiple forward passes. Furthermore, they require specially tailored training algorithms and carefully constructed structures of the nn.

Another limitation of methods relying on creating ensembles is that they have trouble representing the uncertainty caused by the bias in the prediction from a model mismatch. The bias can be caused by an insufficiently flexible model, which could be a result of too high regularization or too low model order.

The problem can be solved by nns from the second category, i.e., where the structure of the nn is extended such that it learns its own uncertainty in the prediction. However, this requires a more intricate nn structure with tailored loss functions [23, 24, 25, 26]. As a consequence, the training becomes more complex and computationally expensive. It also makes the methods sensitive to errors caused by the training algorithm, which are not possible to learn. Furthermore, there is also a need for more data to train complex model structures.

In this paper, we address the two limitations of the aforementioned methods using classical local approximations from the area of system identification [29], which is sometimes referred to as the _delta method_[30, 31, 32]. For regression tasks, the delta method has previously been used to quantify the uncertainty in the prediction of nns, see e.g., [31, 32, 33, 34, 35, 36], and extended to classification tasks in [30].

## 2 Problem formulation and contributions

Consider the problem of learning a classifier from the training data set

\[\mathcal{T}\triangleq\{y_{n},x_{n}\}_{n=1}^{N} \tag{1}\]

Here \(y_{n}\in\{1,\ldots,M\}\) is the class labels and \(x_{n}\!\in\!\mathbb{R}^{n_{x}}\) is the input data of size \(n_{x}\), e.g., pixels in an image. From a statistical point of view, the learning of the classifier can be seen as a system identification problem where a model \(f(x;\theta)\) that predicts the conditional probability mass function (pmf) \(p(y|x)\) of a categorical distribution, are to be identified. That is, the probability for \(y=m\) given the input \(x\) is modeled as

\[p(y=m|x;\theta)=f_{m}(x;\theta),\quad m=1,\ldots,M \tag{2}\]

Here \(\theta\!\in\!\mathbb{R}^{n_{\theta}}\) denote the \(n_{\theta}\)-dimensional parameter vector that parameterize the model. Further, the subscript \(m\) denotes the \(m\):th element of the vector-valued output of the function.

To ensure that the model \(f(x;\theta)\) fulfills the properties associated with pmf, i.e., \(f_{m}(x;\theta)\geq 0\)\(\forall m\) and \(\sum_{m}f_{m}(x;\theta)=1\), it is typically structured as

\[f(x;\theta)=\text{softmax}\left(g(x;\theta)\right) \tag{3}\]

where

\[\text{softmax}(z)\triangleq\frac{1}{\sum_{m=1}^{M}e^{z_{m}}}\begin{bmatrix}e^ {z_{1}}\\ \vdots\\ e^{z_{M}}\end{bmatrix} \tag{4}\]

and \(g(x;\theta)\) describes the underlying model of the classifier.

In the case \(g_{m}(x;\theta)=\theta^{\top}\phi_{m}(x)\), where \(\phi_{m}(x)\) denotes, a possible nonlinear, transformation of the input \(x\), then the model in (3) becomes a standard multinomial logistic regression model [37]. Furthermore, if the transformation \(\phi_{m}(x)\) is chosen randomly, the model becomes similar to the one used in extreme learning machine classifiers [38].

If a nn is used for classification, then the model is given by

\[h^{(0)}=x, \tag{5a}\] \[a^{(l+1)}=\left(h^{(l)}\ 1\right)^{\top}W^{(l)},\quad l=0, \ldots,L-1,\] (5b) \[h^{(l)}=\sigma\big{(}a^{(l)}),\quad l=1,\ldots,L-1,\] (5c) \[g(x;\theta)=a^{(L)}. \tag{5d}\]

Here \(\sigma(\cdot)\) denotes the activation function, where the ReLu function \(\sigma(z)=\max(0,z)\) is often used. The latent variable \(a^{(l)}\) denotes the value of all the nodes in the \(l\)'th layer of the nn, and \(h^{(l)}\) denotes the transformation using the activation function of the values in all the nodes in the \(l\)'th layer of the nn. The parameters of the nn model consist of all the weights and biases included in the matrices \(W^{(L)},\ldots,W^{(0)}\), i.e.,

\[\theta=\left[\text{Vec}(W^{(L)})^{\top}\ \ldots\ \text{Vec}(W^{(0)})^{\top} \right]^{\top}. \tag{5e}\]

Here \(\text{Vec}(\cdot)\) denotes the vectorization operator.
### Parameter estimation

For most nn the number of model parameters \(n_{\theta}>N\) and the model parameters \(\theta\) cannot be uniquely identified from the training data \(\mathcal{T}\) without some regularization or prior information regarding the parameters. Let \(p(\theta)\) denote the prior for the model parameters. The maximum a posteriori estimate of the model parameters is then given by

\[\hat{\theta}_{N}=\arg\max_{\theta}p(\theta|\mathcal{T})=\arg\max_{\theta}L_{N} (\theta)+\ln p(\theta), \tag{6}\]

where \(p(\theta|\mathcal{T})\) denotes the a posteriori distribution of the parameters and

\[L_{N}(\theta)=\sum_{n=1}^{N}\ln f_{y_{n}}(x_{n};\theta) \tag{7}\]

denotes the cross-entropy likelihood function [37]. Here \(y_{n}\) is used as an index operator for the subscript \(m\) of \(f_{m}(x;\theta)\).

### Prediction and classification

Once the classifier has been learned, i.e., a parameter estimate \(\hat{\theta}_{N}\) has been computed, then for a new input data point \(x^{\star}\) the probability mass function can be predicted as

\[\hat{p}(y^{\star}=m|x^{\star};\hat{\theta}_{N})=f_{m}(x^{\star};\hat{\theta}_{ N}),\quad m=1,\ldots,M \tag{8}\]

and the most likely class can be found as

\[\hat{y}^{\star}=\operatorname*{arg\,max}_{m}f_{m}(x^{\star};\hat{\theta}_{N}). \tag{9}\]

Note that, the full pmf estimate \(f(x;\hat{\theta}_{N})\) is needed both for temporal fusion using several inputs from the same class and fusion over different classifiers. Furthermore, even small probabilities can pose a large risk, e.g., there might be a pedestrian in front of a car even if another harmless object is more likely according to the classifier. Hence, it is important that the prediction \(\hat{p}(y^{\star}=m|x^{\star};\hat{\theta}_{N})\) is accurate. However, it is well known that due to, among other things, uncertainties in the parameter estimates \(\hat{\theta}_{N}\) the disagreement between true and estimated pmf may be significant. Therefore, methods to calibrate the prediction \(\hat{p}(y^{\star}|x^{\star};\hat{\theta}_{N})\) such that it better matches \(p(y^{\star}|x^{\star})\) has been developed.

### Temperature scaling

One of the most commonly used methods to calibrate the predicted pmf is called temperature scaling [39]. In temperature scaling \(g(x;\theta)\) is scaled by a scalar quantity \(T\) before the normalization by the softmax operator. With a slight abuse of notation, introduce

\[f(x^{\star};\hat{\theta}_{N},T)=\operatorname*{softmax}\left(g(x^{\star};\hat{ \theta}_{N})/T\right). \tag{10}\]

Via the temperature scaling parameter \(T\) the variations between the components (classes) in the predicted pmf can be enhanced or reduced. When \(T\to 0\), then \(f(x^{\star};\hat{\theta}_{N},T)\to\vec{e}_{i}\), where \(\vec{e}_{i}\) denotes the \(i\)th standard basis vector, thereby indicating that input \(x^{\star}_{n}\) with total certainty belongs to class \(i\). Similarly, when \(T\to\infty\), then \(f_{m}(x^{\star};\hat{\theta}_{N},T)\to 1/M\;\forall m\), thereby indicating that input \(x^{\star}_{n}\) is equally probable to belong to any of the classes.

Noteworthy is that the temperature scaling is typically done after the parameters \(\theta\) have been estimated. For notational brevity, the dependency on the temperature scaling parameter \(T\) will only be explicitly stated when temperature scaling is considered.

### Marginalization of parameter uncertainties

A more theoretically sound approach to take the uncertainties in the parameter estimate into account is via marginalization of the pmf with respect to the parameter distribution. That is, an estimate of the pmf and its covariance are calculated as

\[f(x^{\star}|\mathcal{T})\triangleq\int_{\theta}f(x^{\star}; \theta)p\big{(}\theta|\mathcal{T}\big{)}d\theta \tag{11a}\] \[P^{f}\triangleq\int_{\theta}\!\big{(}f(x^{\star};\theta)-f(x^{ \star}|\mathcal{T})\big{)}\big{(}\big{)}^{\top}\!p\big{(}\theta|\mathcal{T} \big{)}d\theta \tag{11b}\]

From hereon \((x)(\cdot)^{\top}\) is used as shorthand notation for \(xx^{\top}\). The integral in (11a) is generally intractable, but can be approximated by Monte Carlo (mc) sampling as follows

\[\theta^{(k)}\sim p\big{(}\theta|\mathcal{T}\big{)},\quad k=1,2, \ldots,K, \tag{12a}\] \[\hat{f}(x^{\star}|\mathcal{T})=\frac{1}{K}\sum_{k=1}^{K}f(x^{ \star};\theta^{(k)})\] (12b) \[\hat{P}^{f}=\frac{1}{K}\sum_{k=1}^{K}\!\big{(}f(x^{\star};\theta^ {(k)})-\hat{f}(x^{\star}|\mathcal{T})\big{)}\big{(}\big{)}^{\top}\!. \tag{12c}\]

Here \(K\) denotes the number of samples used in the mc sampling.

### Challenges and contributions

To realize the mc scheme in (12) the posterior parameter distribution \(p\big{(}\theta|\mathcal{T}\big{)}\) must be computed and samples drawn from this high-dimensional distribution. Our
contributions are: (i) a local linearization approach that leads to a recursive algorithm of low complexity to compute an approximation of the posterior parameter distribution \(p\big{(}\theta|\mathcal{T}\big{)}\) during the training phase; (ii) a second local linearization approach to reduce the sampling space from \(n_{\theta}\) to \(M\)-dimensional space in the prediction phase; and as a by-product (iii) a low-complexity method for risk assessment and information fusion.

## 3 Posterior parameter distribution

Next, a local linearization approach that leads to a recursive algorithm of low complexity to compute an approximation of the posterior parameter distribution \(p\big{(}\theta|\mathcal{T}\big{)}\) during the training phase is presented.

### Laplace approximation

Assume the prior distribution for the model parameters to be normal distributed as \(p(\theta)=\mathcal{N}(\theta;0,P_{0})\), i.e., \(l^{2}\) regularization is used. Then a Laplace approximation of the posterior distribution \(p(\theta|\mathcal{T})\) yields that [40]

\[p(\theta|\mathcal{T})\approx\mathcal{N}(\theta;\hat{\theta}_{N},P_{N}^{\theta}), \tag{13}\]

where

\[P_{N}^{\theta}=\left(-\frac{\partial^{2}L_{N}(\theta)}{\partial\theta^{2}} \Bigg{|}_{\theta=\hat{\theta}_{N}}+P_{0}^{-1}\right)^{-1}. \tag{14}\]

That is, the prior distribution is approximated by a normal distribution with a mean located at the maximum a posteriori estimate and a covariance dependent upon the shape of the likelihood function in the vicinity of the estimate. The accuracy of the approximation will depend upon the amount of information in the training data \(\mathcal{T}\).

### Asymptotic distribution

According to Bernstein-von Mises theorem [41], if the true model belongs to the considered model set, the maximum a posteriori estimate \(\hat{\theta}\) converge in distribution to

\[\hat{\theta}_{N}\stackrel{{ d}}{{\longrightarrow}}\mathcal{N}( \hat{\theta}_{N};\theta_{*},\mathcal{I}_{\theta}^{-1}), \tag{15}\]

when the information in the training data \(\mathcal{T}\) tends to infinity. Here, \(\theta_{*}\) denotes the true parameters and

\[\mathcal{I}_{\theta}\triangleq-\mathrm{E}\bigg{\{}\frac{\partial^{2}L_{N}( \theta)}{\partial\theta^{2}}\bigg{\}}, \tag{16}\]

is the Fisher information matrix. Given the likelihood function in (7) the Fisher matrix becomes

\[\mathcal{I}_{\theta}\simeq\sum_{n=1}^{N}\sum_{m=1}^{M}\eta_{m,n}\frac{ \partial g_{m}(x_{n};\theta)}{\partial\theta}\bigg{(}\frac{\partial g_{m}(x_{ n};\theta)}{\partial\theta}\bigg{)}^{\!\!\top}\] (17a) where \[\eta_{m,n}\triangleq f_{m}(x_{n};\theta)(1-f_{m}(x_{n};\theta)). \tag{17b}\]

See derivations in Appendix A.

### Recursive computation of covariance

To compute the parameter covariance \(P_{N}^{\theta}\) defined by (14), the Hessian matrix of the log-likelihood (ll) must be calculated and then inverted. This has a complexity of \(\mathcal{O}(NMn_{\theta}^{2}+n_{\theta}^{3})\), which for large \(n_{\theta}\) and \(N\) can become intractable. However, by approximating the Hessian matrix of the ll with the Fisher information matrix as follows

\[P_{N}^{\theta}\approx\left(\mathcal{I}_{\hat{\theta}_{N}}+P_{0}^{-1}\right)^{ -1}, \tag{18}\]

the computation can be done recursively and with a complexity of \(\mathcal{O}\big{(}NMn_{\theta}^{2}+NM^{3}\big{)}\). To do so, note that the \(\mathcal{I}_{\theta}\) in (17) can be written in a quadratic form by defining

\[u_{m,n}\triangleq\sqrt{\eta_{m,n}}\frac{\partial g_{m}(x_{n};\theta)}{ \partial\theta}\Big{|}_{\theta=\hat{\theta}_{N}}. \tag{19}\]

To compute \(u_{mn}\in\mathbb{R}^{n_{\theta}}\) only the gradient of the ll in (7) is required, which is nevertheless needed for the estimation of \(\theta\). Since \(\mathcal{I}_{\theta}\), and so also the covariance \(P_{N}^{\theta}\), can be written in a quadratic form, it is possible to update it recursively as [30]

\[K_{n} =P_{n}^{\theta}U_{n}\big{(}I_{M}+U_{n}^{\top}P_{n}^{\theta}U_{n} \big{)}^{-1} \tag{20a}\] \[P_{n+1}^{\theta} =P_{n}^{\theta}-K_{n}U_{n}^{\top}P_{n}^{\theta}, \tag{20b}\]

where \(I_{r}\) denotes the identity matrix of size \(r\). Here \(P_{n}^{\theta}\) is the parameter covariance for \(n\) measurements, and \(U_{n}\) is defined as

\[U_{n}=\left[u_{1,n}\ \dots\ u_{M,n}\right]\in\mathbb{R}^{n_{\theta}\times M}. \tag{21}\]

The recursion is initialized with \(P_{0}^{\theta}=P_{0}\).

### Approximating the covariance

An nn often has millions of parameters which might result in the amount of data needed to store \(P_{N}^{\theta}\) being larger than the available memory capacity. A common approach to handle this is to approximate \(P_{N}^{\theta}\) as a block-diagonal matrix [42]. Another common approach is to use the approximation

\[P_{N}^{\theta}\approx\begin{bmatrix}P_{N}^{\theta_{r}}&0\\ 0&0\end{bmatrix}, \tag{22}\]

where \(P_{N}^{\theta_{r}}\) denotes the covariance of the estimated parameters \(\theta_{r}\) corresponding to the weights and biases of the \(r\) last layers in the nn[30, 43]. Depending of the
number of included layers, this approximation might be more or less accurate. To compensate for the approximation error when doing the marginalization in (11), a scaling of \(P_{N}^{\theta}\) with factor \(T_{c}\geq 1\) can be introduced. The scaling can be estimated from validation data in a similar manner to the temperature scaling \(T\) in Sec. 2.3.

## 4 Efficient MC sampling

With access to the parameter covariance, one can propagate the uncertainty in the parameters to uncertainty in the prediction with the delta method using the principle of marginalization. Plugging in the approximate Gaussian distribution (15) into (11a) gives

\[f(x^{\star}|\mathcal{T})=\int_{\theta}f(x^{\star};\theta)\mathcal{ N}\big{(}\theta;\hat{\theta}_{N},P_{N}^{\theta}\big{)}d\theta \tag{23a}\] \[P^{f}=\int_{\theta}\!\!\big{(}f(x^{\star};\theta)-f(x^{\star}| \mathcal{T})\big{)}\big{(}\cdot\big{)}^{\top}\mathcal{N}\big{(}\theta;\hat{ \theta}_{N},P_{N}^{\theta}\big{)}d\theta \tag{23b}\]

from which mc approximation can be performed

\[\theta^{(k)}\sim\mathcal{N}\big{(}\theta;\hat{\theta}_{N},P_{N} ^{\theta}\big{)},\quad k=1,2,\ldots,K, \tag{24a}\] \[\hat{f}(x^{\star}|\mathcal{T})=\frac{1}{K}\sum_{k=1}^{K}f\big{(} x^{\star};\theta^{(k)}\big{)}\] (24b) \[\hat{P}^{f}=\frac{1}{K}\sum_{k=1}^{K}\!\big{(}f(x^{\star};\theta^ {(k)})-\hat{f}(x^{\star}|\mathcal{T})\big{)}\big{(}\cdot\big{)}^{\top}. \tag{24c}\]

This is a feasible solution to the problem, but it comes with a high computational cost since it requires drawing mc samples from a high-dimensional Gaussian distribution and evaluating the whole network.

### Marginalization using the delta method

The delta method, see e.g., [31, 32], relies on linearization of the nonlinear model \(g(x,\theta)\) and provides a remedy to the problem of sampling from the high-dimensional Gaussian distribution. The idea is to project the uncertainty in the parameters to uncertainty in the prediction before the softmax normalization (4), thereby drastically reducing the dimension of the distribution that must be sampled. Using the delta method, the uncertainty in the parameters can be propagated to the prediction before the softmax normalization as

\[p(g(x^{\star};\theta)|\mathcal{T})\approx\mathcal{N}\big{(}g(x^{\star};\theta );\hat{g}_{N},P_{N}^{g}\big{)}\] (25a) where \[\hat{g}_{N}=\mathrm{E}\{g(x^{\star};\theta)\}\simeq g(x^{\star};\hat{\theta}_ {N}) \tag{25b}\]

and

\[P_{N}^{g} =\mathrm{Cov}\{g(x^{\star};\theta)\}\] \[\simeq\bigg{(}\frac{\partial}{\partial\theta}g(x^{\star};\theta) \big{|}_{\theta=\hat{\theta}_{N}}\bigg{)}^{\top}P_{N}^{\theta}\frac{\partial} {\partial\theta}g(x^{\star};\theta)\big{|}_{\theta=\hat{\theta}_{N}}. \tag{25c}\]

Using this Gaussian approximation of the parameter distribution, the mc approximation of the marginalization integral becomes

\[g^{(k)}(x^{\star}) \sim\mathcal{N}\big{(}g(x^{\star},\theta);\hat{g}_{N},P_{N}^{g} ),\quad k{=}1,2,\ldots K \tag{26a}\] \[f^{(k)}(x^{\star}) =\mathrm{softmax}\big{(}g^{(k)}(x^{\star})\big{)},\] (26b) \[\hat{f}(x^{\star}|\mathcal{T}) =\frac{1}{K}\sum_{k=1}^{K}f^{(k)}(x^{\star}),\] (26c) \[\hat{P}^{f} =\frac{1}{K}\sum_{k=1}^{K}\!\big{(}f^{(k)}(x^{\star})-\hat{f}(x^{ \star}|\mathcal{T})\big{)}\big{(}\cdot\big{)}^{\top}. \tag{26d}\]

To summarize, the main idea of the delta method is linearization performed in two steps. First, the parameter uncertainty is computed using (15), and second, the uncertainty is propagated to the output of the model by (25). Hence, the delta method is a local linear approach that gives a linear approximation of a nonlinear model.

### Fusion

Suppose there are a set of independent classifiers, each one providing a marginal distribution \(\mathcal{N}\big{(}g_{N,c};\hat{g}_{N,c},P_{N,c}^{\theta}\big{)}\), \(c=1,\ldots,C\). Then the predictions (before the softmax normalization) from these classifiers can be fused as follows [44]

\[P_{N}^{g} =\left(\sum_{c=1}^{C}\big{(}P_{N,c}^{g}\big{)}^{-1}\right)^{-1}, \tag{27a}\] \[\hat{g}_{N} =P_{N}^{g}\sum_{c=1}^{C}\big{(}P_{N,c}^{g}\big{)}^{-1}\hat{g}_{N, c}. \tag{27b}\]

If a single classifier is used to classify a set of inputs \(x_{c}^{\star}\), \(c=1,\ldots,C\), known to belong to the same class \(y^{\star}\), then these predictions can be fused as follows

\[P_{N}^{g} =(H^{\top}R^{-1}H)^{-1}, \tag{28a}\] \[\hat{g}_{N} =P_{N}^{g}H^{\top}R^{-1}z \tag{28b}\]

where

\[z=\begin{bmatrix}\hat{g}_{N,1}\\ \vdots\\ \hat{g}_{N,C}\end{bmatrix}\in\mathbb{R}^{CM}\quad H=\begin{bmatrix}I_{M}\\ \vdots\\ I_{M}\end{bmatrix}\in\mathbb{R}^{CM,M} \tag{28c}\]
and the block \([R]_{i,j}\in\mathbb{R}^{M,M}\), \(i,j=1,\ldots,C\), of the covariance matrix is given by

\[[R]_{i,j}=\frac{\partial}{\partial\theta}g(x_{i}^{\star};\theta)^{\top}\big{|}_{ \theta=\hat{\theta}_{N}}P_{N}^{\theta}\frac{\partial}{\partial\theta}g(x_{j}^ {\star};\theta)\big{|}_{\theta=\hat{\theta}_{N}}. \tag{28d}\]

After fusion, the mc sampling in (26) can be applied as before to compute the pmf estimate.

### Risk assessment

Risk assessment can be defined as the probability \(r_{m}\) that \(p(y_{n}^{\star}=m|x_{n}^{\star})>\gamma_{m}\). The probability \(r_{m}\) can be estimated from the identified model \(f_{m}(x_{n}^{\star}|\mathcal{T})\) as follows

\[\begin{split}\hat{r}_{m}&=\Pr\{f_{m}(x_{n}^{\star} |\mathcal{T})>\gamma_{m}\}\\ &\simeq\frac{1}{K}\sum_{k=1}^{K}\mathbb{1}\big{(}f_{m}^{(k)}(x_{ n}^{\star})>\gamma_{m}\big{)}.\end{split} \tag{29}\]

Here \(\mathbb{1}(a>b)\) denotes the indicator function which is one if \(a>b\) and zero otherwise.

## 5 Validation

Suppose now we have a validation data set \(\mathcal{V}=\{y_{n}^{\circ},x_{n}^{\circ}\}_{n=1}^{N_{\circ}}\). How can we validate the estimated pmf \(\hat{f}(x_{n}^{\circ}|\mathcal{T})\) obtained from (26)? The inherent difficulty is that the validation data, just as the training data, consists of inputs and class labels, not the actual pmf. Indeed, there is a lack of unified qualitative evaluation metrics [14]. That being said, some of the most commonly used metrics are classification accuracy, ll, Brier score, and expected calibration error (ece).

Both the negative ll and the Brier score are proper scoring rules, meaning that they emphasize careful and honest assessment of the uncertainty, and are minimized for the true probability vector [45]. However, neither of them is a measure of the calibration, i.e., reliability of the estimated pmf. Out of these metrics, only ece considers the calibration. Hence, here ece is the most important metric when evaluating a method used to measure the uncertainty [39, 46]. The calculation of the Brier score and ece, together with reliability diagrams are described next. They all can be used to tune the temperature scaling \(T\) described in Section 2.3.

### Brier score

The Brier score [45, 47] corresponds to the least squares fit

\[\frac{1}{N_{\circ}}\sum_{n=1}^{N_{\circ}}\sum_{m=1}^{M}\big{(}\delta_{m,y_{n} }-\hat{p}(y_{n}^{\circ}=m|x_{n}^{\circ})\big{)}^{2}, \tag{30}\]

where \(\delta_{i,j}\) denotes the Kronecker delta function. Furthermore, \(\hat{p}(y_{n}^{\circ}=m|x_{n}^{\circ})\) denotes a generic pmf estimate.

### Accuracy and reliability diagram

Accuracy and reliability diagrams are calculated as follows. Calculate the \(J\) bin histogram defined as

\[B_{j}=\left\{n:\frac{j-1}{J}\leq\max_{m}\hat{p}(y_{n}^{\circ}=m|x_{n}^{\circ} )<\frac{j}{J}\right\} \tag{31}\]

from the validation data. For a perfect classifier \(B_{j}=\emptyset\) for \(j<J\). For a classifier that is just guessing, all sets are of equal size, i.e., \(|B_{j}|=|B_{i}|\ \forall i,j\). Note that \(\max_{m}\hat{p}(y_{n}^{\circ}=m|x_{n}^{\circ})\geq 1/M\), so the first bins will be empty if \(J>M\).

The accuracy of the classifier is calculated by comparing the size of each set with the actual classification performance within the set. That is,

\[\text{acc}(B_{j})=\frac{1}{|B_{j}|}\sum_{n\in B_{j}}\mathbb{1}\big{(}\hat{y}_ {n}^{\circ}=y_{n}^{\circ}\big{)}\] (32a) where \[\hat{y}_{n}^{\circ}=\operatorname*{arg\,max}_{m}\hat{p}(y_{n}^{\circ}=m|x_{n} ^{\circ}) \tag{32b}\]

A reliability diagram is a plot of the accuracy versus the confidence, i.e., the predicted probability frequency. A classifier is said to be calibrated if the slope of the bins is close to one, i.e., when \(\text{acc}(B_{j})=(j-0.5)/J\).

### Confidence and expected calibration error

Instead of certainty, from hereon the standard, and equivalent, notion of confidence will be used [39, 46]. The mean confidence in a set is denoted \(\text{conf}(B_{j})\) and is defined as

\[\text{conf}(B_{j})=\frac{1}{|B_{j}|}\sum_{n\in B_{j}}\max_{m}\hat{p}(y_{n}^{ \circ}=m|x_{n}^{\circ}), \tag{33}\]

This is a measure of how much the classifier trusts its estimated class labels. In contrast to the accuracy it does not depend on the annotated class labels \(y_{n}\). Comparing accuracy to confidence gives the ece, defined as

\[\textsc{ece}=\sum_{j}^{J}\frac{1}{|B_{j}|}|\text{acc}(B_{j})-\text{conf}(B_{j} )|. \tag{34}\]

A small value indicates that the weight is a good measure of the actual performance.
## 6 Experiment study

To illustrate the application of the proposed method to quantify uncertainty in the prediction, two datasets were used. First, an nn was trained using the mnist dataset [48] to classify images of handwritten digits. Second, an nn was trained on the cfar10 dataset [49] to classify images of ten different objects including e.g., cars, cats, and aircraft.

### Classification setup

For the two datasets, the structure of the nn was chosen differently. For the mnist dataset, a five-layer nn with fully-connected nodes were used. For the cfar10 dataset, a LeNet5-inspired structure was used with six convolutional layers followed by four fully connected layers. However, for both datasets the three last layers were chosen to have the same structure, i.e., they were fully connected with \(n_{W,L-2}=100\), \(n_{W,L-1}=40\), and \(n_{W,L}=10\). To decrease the size of the parameter covariance used by the delta method, as described in Sec. 3.4 the first part of the nn was assumed fixed and used to create high-level features. Since the structure of the later layers was chosen identically, the parametric models trained on the two datasets had \(n_{\boldsymbol{\theta}}=4450\) parameters.

To estimate the model parameters \(\theta\) of the nn, the adam optimizer [50] was used. The standard adam optimizer settings, together with an initial learning rate of \(10^{-4}\) and \(l^{2}\) regularization of \(10^{-4}\), were used. Three and ten epochs were used with the mnist dataset and cfar10 dataset, respectively.

### Illustration of the uncertainties in the predictions

The low-dimensional space of the output from \(g(x_{n}^{\circ};\hat{\theta}_{N})\) is particularly interesting to study when trying to understand how the uncertainty in the parameter estimate \(\hat{\theta}_{N}\) affects the classification. Even if the parameter covariance \(P_{N}^{\theta}\) is constant and only depends on the training data, the covariance \(P_{N}^{g}\) depends on the input \(x_{n}\). Fig. 1 illustrates this via an example where we concentrate our study on the decision between just a subset of the number of classes in the mnist dataset, even though the final decision is over all classes. More generally, for some inputs \(x_{n}^{\circ}\) that are located in a dense region in the space of the training data, the covariance \(P_{N}^{g}\) is small, but for an input \(x_{n}^{\circ}\) that is very far from the training data in some norm, the covariance \(P_{N}^{g}\) can be quite large. This indicates that the parameter estimate is quite sensitive in some directions. That means that the output can also be quite sensitive, and a small change in the parameters can give a completely different output. This can be seen in the two examples on the bottom part of Fig. 1. Even though the estimate of the pmf looks similar (especially for the two classes under consideration), by studying the unnormalized prediction \(g(x_{n}^{\circ};\hat{\theta}_{N})\) it is clear that the prediction in the middle has a higher uncertainty compared to the bottom one.

### Results on quantifying the uncertainty

Six different methods to quantify the uncertainty in the classification, i.e., to estimate \(p(y_{n}^{\circ}=m|x_{n}^{\circ})\), were evaluated. These are:

1. Standard method, i.e., \(\hat{p}(y_{n}^{\circ}=m|x_{n}^{\circ})=f_{m}(x_{n}^{\circ};\hat{\theta}_{N})\).
2. Temp. scaling, i.e., \(\hat{p}(y_{n}^{\circ}=m|x_{n}^{\circ})=f_{m}(x_{n}^{\circ};\hat{\theta}_{N},T)\).
3. Deep ensemble, i.e., \(\hat{p}(y_{n}^{\circ}=m|x_{n}^{\circ})\) is estimated using the ensemble method in [15]; number of trained nns are 50 for mnist and 10 for cfar10.
4. mc-dropout, i.e., \(\hat{p}(y_{n}^{\circ}=m|x_{n}^{\circ})\) is estimated using the ensemble method in [16]; 50 samples of the parameters are used to create the ensemble.
5. Proposed method, i.e., \(\hat{p}(y_{n}^{\circ}=m|x_{n}^{\circ})=\hat{f}_{m}(x_{n}^{\circ}|\mathcal{T})\).
6. Proposed method with scaled covariance, i.e., \(\hat{p}(y_{n}^{\circ}=m|x_{n}^{\circ})=\hat{f}_{m}(x_{n}^{\circ}|\mathcal{T},T _{c})\), but with the covariance \(P_{N}^{g}\) in (25) scaled with a factor \(T_{c}\).

In Fig. 2, the reliability diagram for the six different methods to quantify the uncertainty in the prediction of the nn described in Sec. 6.3 is shown. Neither computing the uncertainty in the prediction using the softmax (i), deep ensembles (iii), mc-dropout (iv), or the proposed method without scaled covariance (v) gives calibrated estimates of the uncertainty. To get well-calibrated estimates of the uncertainty either the proposed method with scaled covariance (vi) or temperature scaling (ii) should be used. Finding \(T\) and \(T_{c}\) is commonly done by minimizing the ece. However, increasing the scaling
factor decreases the ll. Hence, there is a trade-off between high ll and low ece. In Table 1, the accuracy, ll, Brier score, and ece are shown for six different methods to quantify the uncertainty in the prediction of the nn. The methods are evaluated both using the mnist and cfar10 datasets. Table 1 shows that the proposed method attains the lowest ece for both datasets. This while still having reasonably good performance in terms of accuracy, ll, and Brier score.

## 7 Summary and Conclusion

A method to estimate the uncertainty in classification performed by a neural network has been proposed. The method also enables information fusion in applications where multiple independent neural networks are used for classification, or when a single neural network is used to classify a sequence of inputs known to belong to the same class. The method can also be used for statistical risk assessment.

The proposed method is based on a local linear approach and consists of two steps. In the first step, an approximation of the posterior distribution of the estimated neural network parameters is calculated. This is done using a Laplacian approximation where the covariance of the parameters is calculated recursively using the structure of the Fisher information matrix. In the second step, an estimate of the PMF is calculated where the effect of the uncertainty in the estimated parameters is considered using marginalization over the posterior distribution of the parameter estimate. This is done by propagating the uncertainty in the estimated parameters to the uncertainty in the output of the last layer in the neural network using a second local linear approach. The uncertainty in the output of the last layer is approximated as a Gaussian distribution of the same dimension as the number of classes. The PMF and its covariance are then calculated via MC sampling, where samples are drawn from this low-dimensional distribution.

The proposed method has been evaluated on two classical classification datasets; MNIST and CFAR10. Neural networks with standard architectures were used. To handle a large number of parameters in these network architectures, only the parameters of the last layer were considered in the uncertainty computations. The results, in terms of ECE, show that the proposed method in its standard form yielded a similar performance as standard methods which do not take the uncertainty in the estimated parameters into account. However, when using a rescaled parameter covariance matrix, used to compensate for the fact that only the uncertainty from the parameters in the last layers was considered, a significant reduction in the ECE was observed. This indicates that

\begin{table}
\begin{tabular}{l|c c c c|c c c c} \hline \hline  & \multicolumn{4}{c|}{mnist} & \multicolumn{4}{c}{cfar10} \\ Method & acc. \(\uparrow\) & ll (\(10^{3}\)) \(\uparrow\) & Brier score \(\downarrow\) & ece \(\downarrow\) & acc. \(\uparrow\) & ll (\(10^{3}\)) \(\uparrow\) & Brier score \(\downarrow\) & ece \(\downarrow\) \\ \hline Standard \(f(x_{n}^{\circ};\hat{\theta}_{N})\) & 91\% & 7.886 & 0.134 & 1.078 & 83\% & 7.904 & 0.291 & 1.328 \\ Temp. sc. \(f(x_{n}^{\circ};\hat{\theta}_{N},T)\) & 91\% & 7.818 & 0.133 & 0.951 & 83\% & 7.740 & 0.269 & 0.612 \\ Deep ensemble & 96\% & 7.856 & 0.080 & 2.868 & 87\% & 7.834 & 0.191 & 1.479 \\ mc-dropout & 93\% & 7.424 & 0.123 & 2.424 & 81\% & 9.935 & 0.301 & 2.829 \\ Prop. met. \(\hat{f}(x_{n}^{\circ}|\mathcal{T})\) & 91\% & 7.845 & 0.151 & 1.242 & 83\% & 8.176 & 0.243 & 2.140 \\ Prop. met. \(\hat{f}(x_{n}^{\circ}|\mathcal{T},T_{\mathrm{c}})\) & 91\% & 7.763 & 0.151 & 0.821 & 82\% & 7.545 & 0.239 & 0.540 \\ \hline \hline \end{tabular}
\end{table}
Table 1: Computed performance measure for the two datasets. The arrows indicate whether a high or low value is preferable.

Figure 2: Reliability diagrams for prediction on the mnist dataset. The diagrams illustrate the six different methods to measure the confidence in the prediction described in Sec. 6.3. A calibration line is also shown in black.
the proposed method works, but that more advanced low-rank methods to approximate the parameter covariance are needed. This is a direction for future research.

## Acknowledgements

This work is supported by Sweden's innovation agency, Vinnova, through project iQDeep (project number 2018-02700).
