Variational Linearized Laplace Approximation for Bayesian Deep Learning
Luis A. Ortega, Sim\'on Rodr\'iguez Santana, Daniel Hern\'andez-Lobato
TL;DR
This paper introduces a variational sparse Gaussian Process approach to efficiently approximate the Linearized Laplace Approximation for Bayesian deep learning, significantly reducing computational costs while improving predictive uncertainty estimates.
Contribution
It proposes a novel variational sparse GP method for LLA that is computationally efficient and maintains the original DNN output as the predictive mean.
Findings
Outperforms existing LLA variants in predictive quality.
Reduces training time to be independent of dataset size.
Demonstrates effectiveness on regression and classification tasks.
Abstract
The Linearized Laplace Approximation (LLA) has been recently used to perform uncertainty estimation on the predictions of pre-trained deep neural networks (DNNs). However, its widespread application is hindered by significant computational costs, particularly in scenarios with a large number of training points or DNN parameters. Consequently, additional approximations of LLA, such as Kronecker-factored or diagonal approximate GGN matrices, are utilized, potentially compromising the model's performance. To address these challenges, we propose a new method for approximating LLA using a variational sparse Gaussian Process (GP). Our method is based on the dual RKHS formulation of GPs and retains, as the predictive mean, the output of the original DNN. Furthermore, it allows for efficient stochastic optimization, which results in sub-linear training time in the size of the training dataset.…
| ResNet-20 | ResNet-32 | ResNet-44 | ResNet-56 | Mean Rank | |||||||||
| Method | ACC | NLL | ECE | ACC | NLL | ECE | ACC | NLL | ECE | ACC | NLL | ECE | |
| MAP | |||||||||||||
| MF-VI | |||||||||||||
| LLA Diag | |||||||||||||
| LLA KFAC | |||||||||||||
| LLA∗ | |||||||||||||
| LLA∗ KFAC | |||||||||||||
| ELLA | |||||||||||||
| Sampled LLA | |||||||||||||
| VaLLA | |||||||||||||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsGaussian Processes and Bayesian Inference
MethodsGreedy Policy Search · Gaussian Process
Variational Linearized Laplace Approximation for Bayesian Deep Learning
Luis A. Ortega
Simón Rodríguez Santana
Daniel Hernández-Lobato
Abstract
The Linearized Laplace Approximation (LLA) has been recently used to perform uncertainty estimation on the predictions of pre-trained deep neural networks (DNNs). However, its widespread application is hindered by significant computational costs, particularly in scenarios with a large number of training points or DNN parameters. Consequently, additional approximations of LLA, such as Kronecker-factored or diagonal approximate GGN matrices, are utilized, potentially compromising the model’s performance. To address these challenges, we propose a new method for approximating LLA using a variational sparse Gaussian Process (GP). Our method is based on the dual RKHS formulation of GPs and retains as the predictive mean the output of the original DNN. Furthermore, it allows for efficient stochastic optimization, which results in sub-linear training time in the size of the training dataset. Specifically, its training cost is independent of the number of training points. We compare our proposed method against accelerated LLA (ELLA), which relies on the Nyström approximation, as well as other LLA variants employing the sample-then-optimize principle. Experimental results, both on regression and classification datasets, show that our method outperforms these already existing efficient variants of LLA, both in terms of the quality of the predictive distribution and in terms of total computational time.
Machine Learning, ICML
1 Introduction
Deep neural networks (DNNs) have gained widespread popularity for addressing pattern recognition problems due to their state-of-the-art performance in predicting target values from a set of input attributes (He et al., 2016; Vaswani et al., 2017). Despite this great success, DNNs exhibit limitations when computing a predictive distribution that accounts for the confidence in the predictions. Specifically, DNNs result in weak calibration (Guo et al., 2017) and in poor reasoning regarding model uncertainty (Blundell et al., 2015). These issues become particularly critical in risk-sensitive situations like autonomous driving (Kendall & Gal, 2017) and healthcare systems (Leibig et al., 2017) among others.
A Bayesian approach, where probabilities describe degrees of belief in the potential values of the DNN parameters, has proven useful for treating these pathologies (MacKay, 1992b; Neal, 2012; Graves, 2011). Bayes’ rule is used here to get a posterior distribution in the high-dimensional space of DNN parameters. Nevertheless, due to the intractability of the calculations, the exact posterior is often approximated with a simpler distribution. Different techniques can be used for this, including (but not limited to) variational inference (VI) (Blundell et al., 2015), Markov chain Monte Carlo (MCMC) (Chen et al., 2014) and the Laplace approximation (LA) (Mackay, 1992; Ritter et al., 2018).
LA offers several advantages over alternative methods. It leverages the maximum a posteriori (MAP) solution, attainable via back-propagation, along with the inverse of the Hessian of the DNN parameters there. The negative inverse Hessian provides the covariances of a Gaussian posterior approximation and the MAP estimate furnishes the mean. Consequently, LA yields a Gaussian posterior approximation following standard DNN training with back-propagation. This resembles a pre-training step followed by fine-tuning, a common practice in deep learning (Daxberger et al., 2021a).
LA’s counterpart is the necessity of computing the Hessian at the MAP, which becomes prohibitive for large DNNs. To simplify this, the Hessian is often approximated using the generalized Gauss-Newton (GGN) matrix (Martens & Grosse, 2015). While this is known to produce underfitting in LA (Lawrence, 2001), using a linearized DNN for prediction alleviates it (Foong et al., 2019). This approach aligns with the fact that the GGN Hessian estimate is exact when applying LA to a linearized DNN. This method is referred to as Linearized LA (LLA) (Immer et al., 2021).
LLA applies LA on a first-order Taylor approximation of the DNN output. Consequently, in regression problems, the mean of the predictive distribution obtained by LLA coincides with the DNN’s prediction at the MAP estimate. This means that LLA simply provides error bars on the pre-trained DNN’s predictions. Additionally, the Hessian estimate is always positive definite in LLA. This enables the usage of LLA at any point, not limited to the MAP estimate. This enhances the post-hoc nature of LLA w.r.t. that of LA.
LLA has shown competitiveness with other approaches on a variety of uncertainty quantification tasks (Immer et al., 2021; Foong et al., 2019). Nevertheless, for the practical usage of LLA in the context of large DNNs, further approximations on top of GGN are required. For instance, diagonal or Kronecker-factored (KFAC) approximations of the complete GGN matrix are used. This is motivated by the cubic cost of computing posterior variances in LLA w.r.t. the training set size or the number of parameters in the DNN.
Here, we broaden LLA’s usage for DNNs to cases with a substantial number of parameters and training instances. We reinterpret the LLA predictive distribution as that of a Gaussian Process (GP) (Khan et al., 2019), and use a sparse variational GP based on inducing points to approximate it (Titsias, 2009). Standard sparse GPs often alter the predictive mean and hence may deteriorate the accuracy of the pre-trained DNN, which is assumed to have minimal generalization error. To avoid this, we use a dual representation of the GP in the Reproducing Kernel Hilbert Space (RKHS). This enables the sparse approximation to be used only in the computation of the predictive variances of the GP. The result is a sparse GP approximation where the predictive mean coincides with the predictions of the pre-trained DNN, without introducing additional prediction error.
We call our method Variational LLA (VaLLA) and conduct a series of regression and classification experiments to compare its performance and computational cost with that of related methods from the literature. Our comparisons include (i) a method that uses Kronecker and diagonal approximations of the GGN matrix, (ii) a method that uses a Nyström approximation of that matrix (ELLA) (Deng et al., 2022), and (iii) a method that relies on generating samples from LLA’s posterior distribution using the sample-then-optimize principle (Antorán et al., 2023). VaLLA’s performance is comparable or better than that of these other methods and is obtained at a smaller cost. Specifically, VaLLA adopts mini-batch training, resulting in sub-linear cost w.r.t. the training set size. Furthermore, VaLLA often generates predictive distributions that closely resemble those of LLA. Figure 1 shows this in a simple 1-D toy regression problem.
2 Background
Consider the task of inferring an unknown function based on noisy observations at corresponding locations . Deep learning (DL) defines a Neural Network with parameters so that s.t. . Thus, is fully specified by . DL often estimates optimal parameters via back-propagation. Nonetheless, despite the remarkable performance achieved (Vaswani et al., 2017), DL methods lack proper output uncertainty estimation, potentially resulting in over-confident predictions in regions without training data.
In the context of Bayesian inference, the observations are related to via the likelihood, . In regression problems, and the likelihood is typically Gaussian. By contrast, in classification problems, , with the number of classes and the likelihood is Categorical with probabilities given by a soft-max link function. In this case, is a multi-output function with outputs, one per class label.
Bayesian neural networks (BNNs) use a probabilistic framework (Mackay, 1992), establishing a prior over the network parameters and computing the Bayesian posterior, , for predictions111We ignore the dependence on to simplify the notation.. The analytic computation of the posterior is often intractable due to the strong non-linearities of the NN. Most methods resort to an approximate posterior , later used for prediction with and Monte Carlo samples. This approximate distribution captures prediction uncertainty (Bishop, 2006).
The Laplace approximation (LA) builds a Gaussian approximate posterior (Mackay, 1992). This takes the form of , where denotes the MAP solution, i.e., , and is the inverse of the negative Hessian of the log posterior, i.e.,
[TABLE]
Often, an isotropic Gaussian prior is considered so that . Hence
[TABLE]
Due to the intractability of the Hessian in large DNNs, it is common to approximate it with the generalized Gauss-Newton (GGN) matrix (Immer et al., 2021):
[TABLE]
where and . The GGN matrix is guaranteed to be positive definite, which means that need not be a maximum of the log posterior. It can be, e.g., any solution found by early-stopping back-propagation.
The earlier formulation of LA suffers from underfitting (Lawrence, 2001). This is attributed to the fact that the GGN approximation is the true Hessian matrix of the linearized DNN (Immer et al., 2021). This implies a shift between posterior inference and predictions that can be mitigated by also predicting using the linearized model:
[TABLE]
This method is known as the linearized LA (LLA). Despite these approximations, LLA still requires the inversion of which scales cubically with the number of parameters of the DNN. A dual formulation of LLA as a Gaussian Process (GP), described in the next section, scales with cubic cost in the number of training points.
2.1 Gaussian Process (GP) Interpretation of LLA
A linear model with a Gaussian distribution on the model parameters generates a GP with specific mean and covariance functions (Williams & Rasmussen, 2006). Consider the prior . The linearized BNN generates random functions that follow a GP with mean and covariance functions defined as
[TABLE]
Using the approximate posterior from LLA, i.e., , we also obtain a GP for prediction. The mean and covariance functions, providing the predictive mean and variances of the LLA approximation, are in this case:
[TABLE]
The Woodbury formula on , as defined in (3), gives
[TABLE]
Defining as the (scaled) Neural Tangent Kernel (i.e., prior covariance function) of the GP (Immer et al., 2021), and , the covariance function of the GP in (6) takes the expression
[TABLE]
This allows us to interpret the LLA approximate predictive distribution as a posterior GP in function space, with prior covariance function given by . The bottleneck of this interpretation is the evaluation of , with cost. In the case of classification problems with classes, and have size , thus the cost is also cubic in the number of classes or DNN outputs.
2.2 Dual formulation of Gaussian Processes in RKHS
A Reproducing Kernel Hilbert Space (RKHS) is a Hilbert space of functions with the reproducing property: such that . In general, can be infinite-dimensional and uniformly approximate continuous functions on a compact set.
A zero-mean prior GP with posterior mean and covariance functions and , respectively, has a dual representation in a RKHS (Cheng & Boots, 2016). There exists and a linear semi-definite positive operator s.t. , verifying
[TABLE]
As a shorthand and an abuse of notation, we write that , where this refers to a Gaussian measure in an infinite dimensional space, not a Gaussian density.
3 Variational LLA (VaLLA)
Here, we present our proposed method, Variational LLA (VaLLA). VaLLA leverages the use of variational sparse GPs. Thus, we proceed to describe these methods first.
3.1 Variational Sparse GPs
Variational sparse GPs approximate the GP posterior using a GP parameterized by inducing points , each in , and associated process values (Titsias, 2009),
[TABLE]
where , and is fixed. The approximate distribution is obtained by minimizing the KL-divergence . In practice, the minimization problem is transformed into the maximization of the lower bound of the log-marginal likelihood
[TABLE]
which has cost due to the cancellation of the factor since .
Theorem 1**.**
(Cheng & Boots, 2016)* Using a sparse GP approximation with is equivalent to restricting the mean and covariance functions of the dual representation in the RKHS to*
[TABLE]
where the functional defines a linear combination of basis functions as , with and the functional , defines a quadratic expression where such that .
Theorem 1 indicates that the algorithm of Titsias (2009) optimizes a variational Gaussian measure where and are parameterized by a function basis . Cheng & Boots (2017) propose to generalize this so that each of the linear operators is optimized using different bases (sets of inducing points). Let and be two sets of inducing points for the mean and the variance, respectively. The parameterization of Cheng & Boots (2017) is:
[TABLE]
where and are defined as using and , respectively. Now, there are two sets of inducing points, for the mean and for the covariances, respectively. This parameterization is a generalization and cannot be obtained using the approach of Titsias (2009). Here, must be found by optimizing Gaussian measures (Cheng & Boots, 2016):
[TABLE]
where the KL term is:
[TABLE]
with and matrices with the prior covariances among and , respectively.
3.2 Using Decoupled SGP and LLA
We utilize the decoupled reparameterization of sparse GPs to establish a model where the mean of the approximated posterior distribution is anchored to a pre-trained MAP solution. We denote this method as variational LLA (VaLLA).
Proposition 1**.**
If , then exists a set of inducing points and a collection of scalar values such that the dual representation of the sparse Gaussian process defined by
[TABLE]
corresponds to a GP posterior approximation with mean and covariance functions defined as
[TABLE]
where is a set of inducing points, , is a vector with the covariances between and , and verifies , with the distance in the RKHS (see proof in Appendix A).
Proposition 1 implies that if we can find values for and inducing points for the mean s.t. can be made as small as desired. For sufficiently small , , and can be used for prediction instead of . Thus, there is no need to optimize and in (14), and the posterior distribution of VaLLA uses as its mean function. The optimal parameters and can be found by optimizing (14) with and held constant. From the following proposition, computing the optimal value of has cost .
Proposition 2**.**
The value of in Proposition 1 that minimizes Equation (14) is
[TABLE]
where is the noise variance and is a matrix with the prior covariances between and . If , the covariance function of the predictive distribution in (17) is equal to that of the full GP (see proof in Appendix A).
MAP solution and Hilbert space.
Proposition 1 presupposes that . In practice, this need not be the case. Covariance functions such as squared exponential or Matérn are recognized for encompassing the entire space of continuous functions. However, whether holds in general remains unknown. For further discussions on this matter, please refer to Appendix B. From here onwards, we assume that if , then is sufficiently expressive to include a close approximation to . Consequently, can be used as the sparse GP posterior mean.
3.3 Hyper-parameter Tuning and -divergences
The VaLLA predictive mean is anchored to the DNN output, making the maximization of the ELBO in (14) unsuitable for tuning hyper-parameters like the prior variance . Specifically, in a regression scenario with Gaussian noise with variance the first term in the r.h.s. of (14) becomes:
[TABLE]
where is constant. Maximizing (19) w.r.t. results in the prior covariances, , tending to [math]. This makes posterior covariances also tend to [math], effectively cancelling the last term in (19). The KL term in (14) is also optimal and [math] if . The reasoning is that, in sparse GPs, tuning hyper-parameters involves a trade-off between fitting the mean to the training data and reducing the predictive variance of the model. Therefore, in the VaLLA setting, where the predictive mean is fixed, the optimal predictive variance tends to be zero.
To address these issues, we propose an alternative objective to (14) that facilitates hyper-parameter optimization:
[TABLE]
Here is a parameter. Instead of minimizing , this objective minimizes, in an approximate way, the -divergence between and (Li & Gal, 2017). Remarkably, this can be achieved by simply changing the data-dependent term in the objective of (14).
The use of -divergences for approximate inference has been extensively studied (Bui et al., 2017; Villacampa-Calvo & Hernández-Lobato, 2020; Santana & Hernández-Lobato, 2022), with observations suggesting that values of result in better predictive mean estimation. Conversely, values of provide superior predictive distributions in terms of the log-likelihood. Thus, in this work we opt for . In this case, (20) does not promote , unlike (14), as the data-dependent term is the log-likelihood of the training data. An unexpected behavior, however, is that (20) may lead to overfitting. To alleviate this, we employ an early-stopping strategy using a validation set (see Appendix C for further details). Early stopping is also used in other LLA approximations s.a. ELLA (Deng et al., 2022).
Mini-batch Optimization.
The objective in (20) supports mini-batch optimization with cost . For ,
[TABLE]
Here, denotes a mini-batch, and the expectation can be computed in closed form in regression problems. In classification, an approximation is available via using the soft-max approximation of Daxberger et al. (2021b). This sub-linear cost of VaLLA enables its usage in very large datasets.
Prediction.
Predictions for test points are computed using (17) with the DNN output as the mean. is evaluated as in training.
Inducing Points.
The locations of the inducing points are found by optimizing (21) with K-means initialization.
3.4 Limitations of VaLLA
VaLLA is limited by three factors: (i) Computing the predictive distribution at each training iteration involves inverting in (17), with cubic cost in the number of inducing points . Therefore, VaLLA cannot accommodate a very large number of inducing points. (ii) The objective in (20) requires a validation set and the use of early-stopping to effectively optimize the prior variance , thus further increasing training time. (iii) Mini-batch optimization in (21) involves evaluating and . Hence, we require efficient evaluation of the (scaled) Neural Tangent Kernel, and its gradients to find . While there are libraries that use structure in the derivatives for the efficient computation of , these are limited to a few DNN models (Novak et al., 2022). A simple but inefficient approach to evaluate involves computing and storing all full Jacobians in memory, for each mini-batch instance and inducing point. This is tractable in our problems, but makes VaLLA infeasible in very large problems, e.g., ImageNet. Appendix E.3 shows a very efficient layer-by-layer method to get and . However, this requires computing each layer’s contribution to the Jacobian at hand, which is difficult for large and complex DNNs.
4 Related Work
LA for DNNs was originally introduced by Mackay (1992), applying it to small networks using the full Hessian. MacKay (1992a) also proposed an approximation similar to the generalized Gauss-Newton (GGN). The combination of scalable factorizations or diagonal Hessian approximations (Martens & Grosse, 2015; Botev et al., 2017) with the GGN approximation (Martens, 2020) played a crucial role in the resurgence of LA for modern DNNs (Ritter et al., 2018; Khan et al., 2019). Recent works aim to relax the Gaussian assumption of LLA adopting a Riemannian-Laplace approximation, where samples naturally fall into weight regions with low negative log-posterior (Bergamin et al., 2023).
To address the underfitting issue associated with LA (Lawrence, 2001), particularly when combined with the GGN approximation, Ritter et al. (2018) proposed a Kronecker factored (KFAC) LLA approximation. This approach outperforms LA with a diagonal Hessian matrix.
The GP interpretation of LLA (Khan et al., 2019) allows to use GP approximate methods to speed-up the computations. Deng et al. (2022) proposed a Nyström approximation of the true GP covariance matrix by using points chosen at random from the training set. The method, called ELLA, has cost . ELLA also requires computing the costly Jacobian vectors required in VaLLA, but does not need their gradients. Unlike VaLLA, the Nyström approximation needs to visit each instance in the training set. However, as stated by Deng et al. (2022), ELLA suffers from over-fitting. An early-stopping strategy, using a validation set, is proposed to alleviate it. In this case, ELLA only considers a subset of the training data. ELLA does not allow for hyper-parameter optimization, unlike VaLLA. The prior variance must be tuned using grid search and a validation set, which increases training time significantly.
Samples from a GP posterior can be efficiently computed using stochastic optimization, elluding the explicit inversion of the kernel matrix (Lin et al., 2023). This approach can be extended to LLA to generate samples from the GP posterior, avoiding the cost (Antorán et al., 2023). However, this method cannot provide an estimate of the log-marginal likelihood for hyper-parameter optimization. To address this limitation, Antorán et al. (2023) propose using the EM-algorithm, where samples are generated (E-step) and hyper-parameters are optimized afterwards (M-step) iteratively. The EM-algorithm significantly increases computational cost, as generating a single sample is as expensive as training the original DNN on the full data. Finally, the method of Antorán et al. (2023) only considers classification problems.
5 Experiments
We compare VaLLA with other methods using the LLA implementation by Daxberger et al. (2021a). VaLLA utilizes a batch size of . For regression, MNIST and FMNIST problems, we train our own DNN (standard multi-layer perceptron), which is stored for reproducibility. In the CIFAR10 experiments with ResNet, we extract the results for other methods from (Deng et al., 2022) and use the same DNN for VaLLA. Hyper-parameters in all LLA variants (diagonal, KFAC, last-layer LLA) are optimized based on the marginal log-likelihood estimate. Additional experimental details are given in Appendix E. VaLLA’s code is given in the supplementary material and is also available here.
5.1 Synthetic Regression
We compare the predictive distribution of VaLLA with that of LLA (which is considered the optimal method), other LLA variants and ELLA, on the 1-D regression problem of Izmailov et al. (2020). In ELLA and VaLLA, we use the optimal hyper-parameters from LLA. The results in Figure 1 illustrate that VaLLA’s predictive distribution closely aligns with that of LLA. Figure 7 (see Appendix D) depicts the predictive distributions of VaLLA and ELLA for varying numbers of inducing points and points in the Nyström approximation, respectively. It shows that VaLLA converges to the true posterior faster than ELLA, with VaLLA tending to over-estimate the predictive variance while ELLA under-estimates it. In Figure 6 (see Appendix D) we observe the effect of tuning the prior variance in VaLLA in another toy 1-D problem, with and without early-stopping. Notably, early stopping, using a validation set, prevents overly small predictive variances in VaLLA. Finally, we observe that when VaLLA estimates the prior variance by maximizing (20), it tends to underestimate LLA’s predictive variance.
5.2 Airline, Year and Taxi Regression Problems
We carry out experiments on large regression datasets. (i) The Year dataset (UCI) with instances and features. We use the original train/test splits. (ii) The US flight delay (Airline) dataset (Dutordoir et al., 2020). Following Ortega et al. (2023) we use the first instances for training and the next for testing. features are considered: month, day of month, day of week, plane age, air time, distance, arrival time and departure time. (iii) The Taxi dataset, with data recorded on January, 2023 (Salimbeni & Deisenroth, 2017). attributes are considered: time of day, day of week, day of month, month, PULocationID, DOLocationID, distance and duration. We filter trips shorter than 10 seconds and larger than 5 hours, resulting in million instances. The first is used as train data, the next as validation data, and the last (10%) as test data. In all experiments, a 3-layer DNN with 200 units, tanh activations and L2 regularization is considered. VaLLA and ELLA use 100 inducing points and 100 random points, respectively.
The table displayed on the l.h.s. of Figure 2 presents the averaged results over 5 random seeds. LLA is not considered here due to intractability. Negative log likelihood (NLL), continuous ranked probability score (CRPS) (Gneiting & Raftery, 2007) and a centered quantile metric (CQM), described below, are reported. We observe that VaLLA performs best according to NLL and CQM, while it gives worse results in terms of CRPS compared to the other methods.
Centered Quantile Metric (CQM).
In regression, CQM assesses the calibration of the predictions, extending the Expected Calibration Error (ECE) to regression problems with Gaussian predictions with the same mean but different variance. CQM calculates the centered interval around the mean with probability mass . For Gaussian predictions , the open interval is defined as , where with the CDF of a Gaussian with mean and variance . The fraction of test points falling inside the interval is then computed. If the predictive distribution is well calibrated, . Formally,
[TABLE]
All methods utilize the pre-trained DNN solution as predictive mean. Thus, the differences in stem from the predictive variance. Evaluating the integrand in (22) on a grid of values allows us to visually interpret the uncertainty estimation of each method. The r.h.s. of Figure 2 shows for several models on the Taxi dataset. In general, all methods tend to over-estimate the actual predictive variance, as evidenced by the values above the diagonal. The l.h.s. of Figure 2 shows CQM estimated using trapezoid integration with points. We refer to Appendix F for more details on the CQM metric.
5.3 Image Classification Problems
MNIST and FMNIST.
We employ a 2-layer fully connected DNN with units in each layer and tanh activations. In VaLLA we considered 100 and 200 inducing points, while in ELLA, 2000 random points are used. The Out-of-distribution (OOD) detection ability of each method is evaluated using the entropy of the predictive distribution as a score. We compute the area under the ROC curve (AUC) of the binary problem that distinguishes between instances from pairs of datasets MNIST/FMNIST and FMNIST/MNIST (Immer et al., 2021). Moreover, in FMNIST we also assess the robustness of the predictive distribution by rotating the test images up to degrees and computing the ECE and NLL on rotated images (Ovadia et al., 2019).
The left table in Figure 3 shows the results on MNIST. VaLLA gives better uncertainty estimates in terms of NLL and the Brier score, but performs less effectively in terms of ECE. Remarkably, VaLLA improves prediction accuracy (ACC) due to the approximation of Daxberger et al. (2021b) to compute class probabilities in multi-class problems. In terms of OOD-AUC VaLLA outperforms the MAP solution but lags behind other methods s.a. Sampled-LLA or LLA with Kronecker approximations. Figure 3 (right) illustrates the training times for each method, with VaLLA being faster than ELLA, Sampled-LLA or Last-Layer LLA.
Finally, the left table in Figure 4 displays the results on FMNIST. Here, VaLLA excels in prediction accuracy and provides the best uncertainty estimates in terms of NLL and the Brier score. Although it does not perform as well in terms of ECE, the differences are small. VaLLA also achieves the best results in OOD-AUC. Figure 4 (right) shows VaLLA holds better performance in terms of ECE and NLL as the test images’ corruption increases (rotation level), indicating the greater robustness of VaLLA’s predictive distribution.
CIFAR10 and ResNet.
Various ResNets architectures are used, and the corresponding pre-trained models are those of Deng et al. (2022) (accessible here). Table 1 shows ACC, NLL and ECE for each method, including LLA variants and a mean-field VI approach (Deng & Zhu, 2023). Since we use the same pre-trained models, the results of all other methods are consistent with those reported by Deng et al. (2022). VaLLA with outperforms other methods in most cases, always being either the best or second-best method. Figure 5 shows the NLL on the perturbed test set with 5 increasing levels of image corruptions (Deng et al., 2022). Each box-plot summarizes the test NLL for each intensity level across all corruptions. The results again highlight VaLLA’s robust predictive distribution, achieving also lower NLL compared to the other methods.
6 Conclusions
We introduced VaLLA, a method derived from the formulation of a generalized sparse GP that offers the flexibility to fix the predictive mean to any desired function in the RKHS. VaLLA excels in computing error bars for pre-trained DNNs with a vast number of parameters on extensive datasets, handling even millions of training instances. VaLLA’s applicability spans both regression and classification problems, showcasing costs independent of the number of training points . In comparison, the Nyström approximation by Deng et al. (2022) incurs a linear cost in , unless early-stopping is employed. Furthermore, VaLLA surpasses the sample-then-optimize method of Antorán et al. (2023) in terms of speed, while also providing predictive distributions robust to input corruptions. In essence, VaLLA stands out by delivering robust predictive distributions akin to LLA, all while maintaining noteworthy computational efficiency.
Acknowledgements
Authors gratefully acknowledge the use of the facilities of Centro de Computacion Cientifica (CCC) at Universidad Autónoma de Madrid. The authors acknowledge financial support from project PID2022-139856NB-I00 funded by MCIN/ AEI / 10.13039/501100011033 / FEDER, UE and project PID2019-106827GB-I00 / AEI / 10.13039/501100011033 and from the Autonomous Community of Madrid (ELLIS Unit Madrid). The authors also acknowledge financial support from project TED2021-131530B-I00, funded by MCIN/AEI /10.13039/501100011033 and by the European Union NextGenerationEU PRTR.
Societal Impact
As machine learning models play an ever-growing role in influencing decisions with substantial consequences for society, industry, and individuals —such as ensuring the safety of autonomous vehicles (mcallister2017concrete) and improving disease detection (sajda2006machine; singh2021better)— it becomes imperative to possess a comprehensive understanding of the employed methodologies and be capable of offering robust performance assurances. Our diligent examination of VaLLA’s performance across diverse datasets and tasks as part of our empirical assessment showcases its adaptability to various domain-specific datasets.
Appendix A Proofs
See 1
Proof.
First of all, notice that , leading to
[TABLE]
Furthermore, defined as can be seen as a vector of , considering the image of any orthonormal basis of . In fact, let be the usual basis of :
[TABLE]
Assume a variational distribution with and . Then, defining the correspondent dual vectors as
[TABLE]
the mean and covariance functions in the dual formulation of the sparse GP are
[TABLE]
and
[TABLE]
which is the same approximate GP posterior found in Equation () of Titsias (2009). ∎
Proposition 1**.**
If , then exists a set of inducing points and a collection of scalar values such that the dual representation of the sparse Gaussian process defined by
[TABLE]
corresponds to a GP posterior approximation with mean and covariance functions defined as
[TABLE]
where is a set of inducing points, , is a vector with the covariances between and , and verifies , with the distance in the RKHS.
Proof.
First of all, if , the reproducing property of the RKHS verifies that there exists , with the input space, and such that verifies
[TABLE]
As a result, the mean function of the approximate posterior is
[TABLE]
On the other hand, using that
[TABLE]
the covariance function is
[TABLE]
Where the characterization of and as elements of and of Equations (23) and (24) are used. ∎
Proposition 2**.**
The value of in Proposition 1 that minimizes Equation (14) is
[TABLE]
where is the noise variance and is a matrix with the prior covariances between and . If , the covariance function of the predictive distribution in (17) is equal to that of the full GP.
Proof.
For simplicity, assume a non-inverse reparameterization where
[TABLE]
We will find the optimal value for and compute the corresponding value for (using the inverse reparameterization). First, we will show that the true GP posterior can be written in the dual formulation as the following
[TABLE]
where it verifies that the GP posterior mean function is
[TABLE]
Where the characterization of and as elements of and is used, similarly to Equations (23) and (24). Furthermore, using Woodbury matrix identity,
[TABLE]
where again, correspondence between operators and matrices is used. This leads to the covariance function:
[TABLE]
This mean and covariance function are exactly the ones obtained from the original GP formulation (Titsias, 2009). Using the ELBO is equivalent to the KL divergence between the true posterior and the variational approximation, we got
[TABLE]
where
[TABLE]
Naming , it is important to notice that given , despite being an operator, can be seen as a matrix whose entries are the application of to the usual basis of matrices. In short
[TABLE]
The partial derivative of w.r.t. a single position in the matrix is , where denotes a zero-vector with a single at position . Then, using the chain rule for matrices,
[TABLE]
we can compute the optimum in the two terms in the KL. The optimum for the logarithm term can be computed as
[TABLE]
leading to
[TABLE]
On the other hand, the optimum for the trace term is
[TABLE]
where we used that . As a result,
[TABLE]
Using all the derivations
[TABLE]
Using the Woodbury matrix identity
[TABLE]
Thus, the value of where verifies
[TABLE]
Using that , we can take further derivations on the expression of as
[TABLE]
Applying again Woodbury matrix identity:
[TABLE]
If we substitute this value on the predictive distribution
[TABLE]
This expression coincides with the optimal sparse GP solution described by Titsias (2009), with optimal variational covariance . As a result, the optimal solution for VaLLA coincides with the optimal solution for standard sparse GPs. Let us now compute the optimal value of given the optimal value of . First, notice that using Woodbury Matrix identity
[TABLE]
Therefore, the relation between and is . Meaning that
[TABLE]
Global optimum
To complete the proof of the solution in Eq. (55) being not only optimal but also a maximum of the ELBO, we must test the behavior of the second derivative w.r.t. . Let us reuse the previous results, where we found that
[TABLE]
Taking a second derivative w.r.t. another location yields
[TABLE]
Considering that does not depend on , the first term drops from the derivative. Thus
[TABLE]
Here we aim to use the chain rule for matrices,
[TABLE]
Then, consider that
[TABLE]
Then,
[TABLE]
We can work on this trace to simplify the expression as
[TABLE]
Naming , we got
[TABLE]
where in the last equality we used that is symmetric. Using the definition of we know that
[TABLE]
This shows that is the posterior covariance of a GP with prior covariances and noise covariances .
[TABLE]
with denoting the Kronecker product and a definite positive matrix. As the Kronecker product of two definite positive matrices is definite positive, the optimal is a minimum of the KL. Given that this second derivative is positive-definite, the value found is the global minimum of the ELBO objective.
∎
Appendix B MAP solution in Hilbert Space
Whether the map solution is in the Hilbert space might be difficult (if not impossible) to know. However, there are cases where it can be theoretically shown; for example in a linear model. Let the map solution be a linear model as
[TABLE]
Then, the features (Jacobians) are
[TABLE]
With a single inducing point and scalar value , the mean function would be
[TABLE]
This recovers the MAP solution if and .
However, if the model is not linear but has a linear last layer as
[TABLE]
where is a non-linear function that depends on parameters . Then, the features (Jacobians) are
[TABLE]
Here it might be difficult to check if there exists a combination that yields the map solution as the mean function.
Appendix C Over-fitting and Early Stopping
In Section 3.3 we discussed the fact that the standard maximization of the ELBO does not allow the optimization of the prior variance. To summarize that section, the optimal value for the prior variance is infinite as a result of the mean being fixed to the optimal MAP solution. As discussed, we circumvent this by applying -divergences, which are not ill-defined in this learning setup; allowing the optimization of the prior. However, the use of this optimization objective is not perfect and we faced the fact that it tends to over-fit the prior variance to the training data. The middle column of Figure 6 (middle) shows the obtained predictive distribution (two times standard deviation) learned from VaLLA using the black points as training data. The MAP solution is obtained using a 2 hidden layer MLP with hidden units and tanh activation, optimized to minimize the RMSE of the training data for iterations with Adam and learning rate . VaLLA on the other hand is trained for iterations. As one may see in the image, the prior variance is fitted to the data to the point where the uncertainty does not increase in the middle gap of the data.
In this experiment, VaLLA optimizes hyper-parameters along with the variational objective. LLA optimizes the prior variance and likelihood variance by maximizing the marginal log likelihood and ELLA uses LLA’s optimal hyper-parameters.
In this situation there are two simple courses of action: we could return to the original ELBO and choose the prior variance by cross-validation; or we could perform early-stopping with the -divergences objective, using a validation set to stop training before over-fitting the prior variance parameter. This last approach may not work as it assumes that there is a point during training where the prior variance truly explains the underlying data without over-fitting. However, as the prior variance is set to a relatively large value compared to the optimal one (which is small and leads to over-fitting), this method resulted in great performance for VaLLA, while avoiding using a costly cross-validation approach. The left column of Figure 6 shows the obtained predictive distribution (two times standard deviation) learned from VaLLA, in this case, using the black points as training data and the orange points as the validation set. In the experiments, we computed the NLL of the validation set every training iterations and stopped training when it worsens. This also allowed us to save computational time.
Appendix D Increasing Inducing Points
Using the optimal covariance in Proposition 2, if the set of inducing points equals the training points , the posterior distribution of VaLLA equals that of the exact LLA Gaussian Process. This suggest that increasing the number of inducing points would lead to better uncertainty estimations. In this section, we aim to show how close is the predictive distribution of VaLLA to that of LLA when we increase the number of inducing points in .
Figure 7 shows the obtained predictive distribution of VaLLA (first row) and ELLA (second row) for , and inducing points/samples. The initial and final locations of the inducing points are also shown for VaLLA. The posterior distribution obtained by LLA is shown in dotted orange. The MAP solution is obtained using a 2 hidden layer MLP with hidden units and tanh activation, optimized to minimize the RMSE of the training data for iterations with Adam and learning rate . VaLLA on the other hand is trained for iterations. For this experiment, VaLLA and ELLA use the optimal prior variance and likelihood variance obtained by optimizing LLA’s marginal log likelihood. As one may see in the image, it is clear that one of the main differences between the two methods is that VaLLA tends to over-estimate the variance whereas ELLA tends to infra-estimate it, compared to LLA. Furthermore, the value of for which the model is closer to the LLA posterior is lower for VaLLA than for ELLA. As we increase , VaLLA’s predictive distribution becomes closer and closer to that of LLA.
The initial and final position of the inducing locations is also shown in the figure. For this experiments, the initial values are computed using K-Means. It can be seen how VaLLA is capable of tuning the inducing locations and move them from one cluster of points to another as needed. This is one of the main advantages of this method compared to ELLA.
Appendix E Experimental Details
Source code for the conducted experiments can be accessed in the following anonymous repository: https://github.com/anonymous-valla/valla.
E.1 MAP solutions
For regression problems (Year, Airline and Taxi datasets), a 3-layer fully connected NN was used with units in each layer. The optimal weights are obtained by minimizing the RMSE using iterations of batch size and Adam optimizer (Kingma & Ba, 2014) with learning rate and weight decay .
For MNIST and FMNIST experiments, a 2-layer fully connected NN was used with units in each layer. The optimal weights are obtained by minimizing the NLL using iterations of batch size and Adam optimizer (Kingma & Ba, 2014) with learning rate and weight decay .
E.2 Laplace Library
The Laplace library (Daxberger et al., 2021a) was used to perform last-layer, KFAC and diagonal approximations of the LLA method and optimize the prior variance on each case. The latter is done by optimizing the log marginal likelihood of the data using the library’s log_marginal_likelihood method for iterations with the Adam optimizer and learning rate .
E.3 Efficient Kernel Computation for MLP
In this section we discuss an efficient implementation for computing the Neural Tangent Kernel . First of all, take into account that the computation of the kernel can be reduced to a summation on the number of parameters of the model:
[TABLE]
One of the limitations of computing the kernel is storing in memory, which is a 3 dimensional tensor of (batch size, number of classes, number of parameters). Computing the kernel as a sum allows to simplify the required computations significantly (we no longer have to store in memory the Jacobians). Consider now a MLP as
[TABLE]
where each function is a non-linear activation function and each function is a linear function of the form
[TABLE]
With this, is supposed to be a fully-connected neural network of layers. Each of the partial derivatives of the neural network are
[TABLE]
and the kernel is computed simply by adding the product of these derivatives. Here, is a sub-index denoting input -th to layer . Similarly, is a sub-index denoting each component of the bias vector parameter at layer , or similarly, each output of that layer.
In fact, using the structure of the model and the chain rule, the derivative of the output of the network w.r.t. the weight parameter of the layer is:
[TABLE]
where each of the two vectors in the r.h.s. has length equal to the number of units in the layer . In fact
[TABLE]
where {\color[rgb]{.75,0,.25}a}({\color[rgb]{0,0,1}h_{l-1}})_{i} corresponds to the inputs of the layer. Moreover, can also be computed using the chain rule:
[TABLE]
which can be easily computed by back-propagating the derivatives. The same derivations can be easily done for the biases of each layer . As a result, the derivatives only depend on a back-propagating term {\color[rgb]{0,.5,.5}\frac{\partial}{\partial h_{l}}g_{o}(\mathbf{x},\hat{\bm{\theta}})} for each layer, the value of the parameters {\color[rgb]{1,.5,0}\bm{W}_{l}},{\color[rgb]{1,.5,0}\bm{b}_{l}} and the propagated outputs at each layer {\color[rgb]{0,0,1}h_{1}},\dots,{\color[rgb]{0,0,1}h_{L-1}} evaluated at the non-linear activation {\color[rgb]{.75,0,.25}a}(\cdot) and its derivative {\color[rgb]{0,1,1}a^{\prime}}(\cdot). This means that, if we store the intermediate outputs of each layer ({\color[rgb]{0,0,1}h_{1}},\dots,{\color[rgb]{0,0,1}h_{L-1}}) on the forward pass of the model, by using a single backward pass, we can compute {\color[rgb]{0,.5,.5}\frac{\partial}{\partial h_{l}}g_{o}(\mathbf{x},\hat{\bm{\theta}})} for each layer.
Critically, given each {\color[rgb]{0,.5,.5}\frac{\partial}{\partial h_{l}}g_{o}(\mathbf{x},\hat{\bm{\theta}})}, we can add the contribution of each layer to the kernel, using (70). In this process, we can sped-up the computations by using structure in the derivatives. For example, in (71) we observe that the derivative has a simple form which is a vector of ones times a scalar. Furthermore, there is no dependence on , the output unit corresponding to the weight . Therefore, for two instances and , the kernel contribution (ignoring the prior variance parameter) corresponding to outputs and is:
[TABLE]
with s_{o,\mathbf{x}}^{l}=\left({\color[rgb]{0,.5,.5}\frac{\partial}{\partial h_{l}}g_{o}(\mathbf{x},\hat{\bm{\theta}})}\right)^{T}\bm{1}_{l} a scalar. Similar simplifications occur in the case of, e.g., a convolutional layer.
Summing up, by using this method, all the required kernel matrices can be easily and efficiently computed, for a mini-batch of data points and a set of inducing points, with a similar cost as that of letting the mini-batch or the inducing points go through the DNN. A disadvantage is, however, that the described computations will have to be manually coded for each different DNN architecture. This becomes tedious in the case of very big DNN with complicated layers, as described in Section 3.4.
Appendix F Further Analysis of the Quantile Metric
As stated in Section 5, we proposed a new metric for regression problems that, in a way, extends ECE to regression problems with Gaussian predictive distributions with the same mean. This kind of metric is desirable for LLA methods as all of them rely on keeping the optimal MAP solution as the predictive mean of the model. They only differ in the predictive variance. Formally, CQM computes for each the probability that points fall into the predictive centered interval of probability . The underlying reasoning is that, if the model explains the data well enough, of the points will fall inside the centered quantile interval. Thus, the metric defined as
[TABLE]
should be roughly [math] when the model predictive distribution is similar to the actual one, given by the observed data.
Figure 8 shows the evolution of w.r.t. for the best performing models in the regression problems. CQM corresponds to the area between the shown curve and (black line). This figure allows to argue that (in general) all methods are over-estimating the predictive variance as they are giving values above the diagonal. That is, for a specific value of , the reported probabilities are higher than , meaning that, on average, there are more points in than they should. That is, the predicted interval is larger than it should, which can only mean that the variance is over-estimated. From a geometrical perspective, it is clear that CQM is always greater than [math] and lower than ; independently of the model and dataset used.
In fact, this figure allow to visually study the level of over/infra-estimation of the prediction uncertainty, for each degree of confidence . For example, in the Year dataset (Figure 8) we see that VaLLA slightly over-estimates the uncertainty for while it infra-estimates it for larger values of .
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Antorán et al. (2023) Antorán, J., Padhy, S., Barbano, R., Nalisnick, E. T., Janz, D., and Hernández-Lobato, J. M. Sampling-based inference for large linear models, with application to linearised Laplace. In International Conference on Learning Representations , 2023.
- 2Bergamin et al. (2023) Bergamin, F., Moreno-Muñoz, P., Hauberg, S., and Arvanitidis, G. Riemannian laplace approximations for bayesian neural networks. ar Xiv preprint ar Xiv:2306.07158 , 2023.
- 3Bishop (2006) Bishop, C. M. Pattern Recognition and Machine Learning (Information Science and Statistics) . Springer, 2006.
- 4Blundell et al. (2015) Blundell, C., Cornebise, J., Kavukcuoglu, K., and Wierstra, D. Weight uncertainty in neural network. In International conference on machine learning , pp. 1613–1622. PMLR, 2015.
- 5Botev et al. (2017) Botev, A., Ritter, H., and Barber, D. Practical gauss-newton optimisation for deep learning. In International Conference on Machine Learning , pp. 557–565. PMLR, 2017.
- 6Bui et al. (2017) Bui, T. D., Yan, J., and Turner, R. E. A unifying framework for gaussian process pseudo-point approximations using power expectation propagation. Journal of Machine Learning Research , 18:1–72, 2017.
- 7Chen et al. (2014) Chen, T., Fox, E., and Guestrin, C. Stochastic gradient hamiltonian monte carlo. In International conference on machine learning , pp. 1683–1691. PMLR, 2014.
- 8Cheng & Boots (2016) Cheng, C.-A. and Boots, B. Incremental variational sparse Gaussian process regression. Advances in Neural Information Processing Systems , 29, 2016.