Paper deep dive
Scalable AI Uncertainty Quantification via Generalized Laplace Active Subspaces
Wouter N. Edeling, Peter V. Coveney
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 10/10/2026, 2:24:38 AM
Summary
This paper proposes a scalable Uncertainty Quantification (UQ) method for neural networks using a low-rank generalized Laplace approximation within active subspaces. The authors compare standard Bayesian scaling (summed loss) with generalized mean-loss scaling. They find that standard scaling causes posterior variance contraction dependent on data size, requiring more curvature directions and degrading predictive coherence. In contrast, mean-loss scaling yields a stable, lower-dimensional active subspace, producing calibrated and coherent predictive confidence intervals without needing to linearize the network.
Entities (10)
Relation Signals (8)
Wouter N. Edeling → affiliatedwith → CWI
confidence 99% · Affiliation: Scientific Computing Group, CWI
Peter V. Coveney → affiliatedwith → University College London
confidence 99% · Affiliation: Centre for Computational Science, University College London
Wouter N. Edeling → affiliatedwith → University of Twente
confidence 95% · Affiliation: Faculty of Electrical Engineering, Mathematics and Computer Science, University of Twente
Generalized Laplace Active Subspaces → proposessolutionfor → Uncertainty Quantification
confidence 95% · We propose a low-rank generalized Laplace approximation for neural-network UQ
Mean-Loss Scaling → yields → Calibrated Predictive Intervals
confidence 92% · the generalized mean-loss scaling yields a more stable, lower dimensional active subspace and produces calibrated, coherent predictive confidence intervals
Standard Bayesian Scaling → causes → Posterior Variance Contraction
confidence 90% · the standard scaling induces a data-size dependent contraction of the posterior variance
Laplace Approximation → istypeof → Posterior Approximation
confidence 90% · the Laplace approximation MacKay (1992) is a post-hoc posterior approximation.
Active Subspaces → →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Reliable uncertainty quantification (UQ) is essential for deploying neural networks in scientific and high-stakes applications, but full Bayesian inference over the network parameters is computationally infeasible. We propose a low-rank generalized Laplace approximation for neural-network UQ based on a small number of data-informed curvature directions. Starting from a generalized Bayesian posterior defined through an empirical loss, we construct a local Gaussian approximation around a pretrained set of weights in this active curvature subspace. The posterior variances in the retained subspace are available in closed form, and the prior variance is calibrated by an empirical Bayes procedure. The generalized Bayesian formulation allows us to compare two posterior scalings: the standard Bayesian scaling associated with the summed negative log likelihood, and a mean-loss scaling in which the empirical loss is normalized by the number of data. A central finding is that the standard scaling induces a data-size dependent contraction of the posterior variance in the leading active directions. In regression problems, this can force the low-rank framework to retain additional weak-curvature directions in order to achieve nominal coverage of calibration data. When posterior samples are propagated through the non-linear network, these additional directions can degrade the coherence of the predictive intervals and shift the posterior predictive mean away from the pretrained model. In contrast, the generalized mean-loss scaling yields a more stable, lower dimensional active subspace and produces calibrated, coherent predictive confidence intervals. These results indicate that generalized Laplace active subspaces provide a practical and scalable route to calibrated uncertainty quantification in neural networks.
Tags
Links
- Source: https://arxiv.org/abs/2610.11738v1
- Canonical: https://arxiv.org/abs/2610.11738v1
Trouble viewing inline? Open PDF directly →
Full Text
94,917 characters extracted from source content.
Expand or collapse full text
Scalable AI Uncertainty Quantification via Generalized Laplace Active Subspaces Wouter N. Edeling Email: edeling@cwi.nl Affiliation: Scientific Computing Group, CWI, Science park 123, Amsterdam, 1098XG, the Netherlands Affiliation: Faculty of Electrical Engineering, Mathematics and Computer Science, University of Twente, Hallenweg 15, Enschede, 7500AE, the Netherlands Peter V. Coveney Email: p.v.coveney@ucl.ac.uk Affiliation: Centre for Computational Science, University College London, 20 Gordon Street, London, WC1H 0AJ, UK Affiliation: Advanced Research Computing Centre, University College London, Gordon Street, London, WC1E 6BT, UK Abstract Reliable uncertainty quantification (UQ) is essential for deploying neural networks in scientific and high-stakes applications, but full Bayesian inference over the network parameters is computationally infeasible. We propose a low-rank generalized Laplace approximation for neural-network UQ based on a small number of data-informed curvature directions. Starting from a generalized Bayesian posterior defined through an empirical loss, we construct a local Gaussian approximation around a pretrained set of weights in this active curvature subspace. The posterior variances in the retained subspace are available in closed form, and the prior variance is calibrated by an empirical Bayes procedure. The generalized Bayesian formulation allows us to compare two posterior scalings: the standard Bayesian scaling associated with the summed negative log likelihood, and a mean-loss scaling in which the empirical loss is normalized by the number of data. A central finding is that the standard scaling induces a data-size dependent contraction of the posterior variance in the leading active directions. In regression problems, this can force the low-rank framework to retain additional weak-curvature directions in order to achieve nominal coverage of calibration data. When posterior samples are propagated through the non-linear network, these additional directions can degrade the coherence of the predictive intervals and shift the posterior predictive mean away from the pretrained model. In contrast, the generalized mean-loss scaling yields a more stable, lower dimensional active subspace and produces calibrated, coherent predictive confidence intervals. These results indicate that generalized Laplace active subspaces provide a practical and scalable route to calibrated uncertainty quantification in neural networks. We demonstrate the method on two regression problems, and a classification case involving an autoregressive molecular language model. keywordsuncertainty quantification, neural networks, Laplace posterior, active subspaces 1 Introduction One of the biggest open questions concerning artificial intelligence (AI) is the reliability and trustworthiness of the predictions which emanate from its output. After training, AI software operating in inference mode can spew out predictions concerning the behaviour of a new system in the same class as those it has been trained on in virtually instantaneous fashion. Such predictions can have enormous value in many real world scenarios. In some situations, qualitative classifications suffice while others require only a modicum of accuracy to be beneficial; still others require accurate and precise predictions or they are not actionable. This last scenario is the one which is most often necessary in scientific and engineering contexts. To use the results of a calculation inside a scientific procedure, one generally needs to know that the predictions are reliable within defined and suitably limited bounds of uncertainty. Actionable predictions must therefore not only be timely; they need to be of a sufficient level of accuracy and/or precision to be used with confidence. The problem with AI methods in general is that what they output often appears plausible and may pass muster as such for some applications; but in more quantitative scientific and mathematical domains, plausibility is not sufficient: we need to be clear about the level of certainty of its predictions. The key assessments of validation, verification and uncertainty quantification are regularly visited on physics-based models, e.g. for weather and climate prediction Palmer (2000); Murphy et al. (2004), turbulence modelling Cheung et al. (2011); Gorlé and Iaccarino (2013), tsunami hazards Liu and Guillas (2017); Selva et al. (2016) or biological applications Avitabile et al. (2026). In these instances there is an exactly defined mathematical model based on clearly defined equations that typically include a variety of parameters, whereas such systematic analysis are significantly more challenging for AI systems. One immediately noticeable difference between AI and physics-based codes resides in the number of parameters which arise in the two cases. Whereas physics-based codes typically contain around five parameters Cheung et al. (2011), possibly several hundreds Edeling et al. (2024) (extending to a few thousand in some large scale codes), each parameter has a well-defined physical meaning pertaining to the underlying model. In AI software there are many orders of magnitude more, none of which has a physical interpretation, as they are simply free parameters. The impressive fitting of input-output data can be readily performed with so many free parameters in play, particularly when one works with held-out data sets where all manner of optimisations can be conducted. The approach lends itself to the making of overconfident predictions, which can and often do unravel when the trained system is applied to data which it has not previously encountered. We are interested in the uncertainty of such predictions in terms of the uncertainty in the parameter values of the trained model; unlike a physics-based approach, lacking any physical meaning, these parameters cannot be meaningfully adjusted based on insights into the system which can be gleaned from mathematical models we have created. The purpose of the present paper is to demonstrate that, despite the lack of clarity as to the meaning of the parameters in an AI system, it is possible to obtain predictive uncertainty estimates. A natural framework for quantifying parameter uncertainty is Bayesian inference Gelman (2013). In principle, one would place a prior distribution on the neural-network weights, update this prior by conditioning on data, and propagate the resulting posterior distribution through the network to obtain predictive uncertainty. In practice, since the parameter dimension is extremely large, standard sampling strategies are computationally infeasible. Scalable approximations are therefore unavoidable. Scalable Bayesian neural-network posterior approximations include variational approaches such as Bayes by Backprop Blundell et al. (2015), which learns a parametric approximation to the posterior distribution over the network weights, typically using a Gaussian variational family. The resulting training objective is an evidence lower bound, consisting of an expected data-fit term and a Kullback-Leibler (KL) divergence to the prior. In practice, the balance between these two terms can be delicate, and may require specialized training strategies such as KL annealing or cyclical schedules Fu et al. (2019). By contrast, the Laplace approximation MacKay (1992) is a post-hoc posterior approximation. It also models the parameter posterior by a Gaussian distribution, but this Gaussian is centred at a pre-trained set of weights, and its covariance is obtained from the local curvature of the log posterior density rather than by optimizing a variational distribution during training. While the posterior is now conceptually simple, the computational bottleneck shifts to the covariance matrix of the Laplace approximation, as it is too large to form explicitly for many networks. In this work we circumvent this bottleneck via a generalized-Bayes low-rank Laplace approximation. The posterior is restricted to a low-dimensional active subspace spanned by the dominant Gauss-Newton curvature directions of the empirical mean loss. These directions are computed using matrix-free curvature-vector products, so that the full covariance matrix is never formed. Although active subspaces are originally derived from the dominant directions of an uncentered covariance matrix of the gradient of a (physics-based) model Constantine et al. (2014), the overarching idea is similar, namely to find a linear projection to a low-dimensional subspace relevant for a quantity of interest. For the original active subspace approach this quantity was the model output subject to a prior distribution, in this work it is the posterior predictive uncertainty of a neural network. Within the subspace spanned by these dominant active directions, the posterior variances are available in closed form, while the prior variance is calibrated using a generalized-evidence criterion. The active subspace dimension is selected by empirical coverage on held-out calibration data, yielding predictive intervals that are both calibrated and sharp. We derive our low-rank Laplace posterior within the framework of generalized Bayesian inference, where the data is connected to the parameters via a loss function, rather than a traditional likelihood function Bissiri et al. (2016). In particular we contrast summed with mean loss functions, where the latter is scaled by 1/N1/N, N being the number of data records in the training set. The use of the summed loss equates to the standard Bayesian likelihood function. A central observation of the paper is that the scaling of the generalized posterior is not merely a technical detail. The standard Bayesian scaling leads to an increasingly concentrated posterior variance with the number of training data. To avoid such overconfident predictions, a low-rank approximation should then retain many weak-curvature directions with relatively high posterior variance. This increases the computational cost, and reduces the statistical coherence of the predictive output distribution. In other words, when samples from a local Gaussian weight posterior are propagated through the nonlinear network, these weak directions can shift the predictive mean or distort the resulting uncertainty intervals. Existing methods Immer et al. (2021); Faller and Martin (2026) deal with this issue by replacing the full network with its linear expansion. By contrast, our generalized Bayesian mean-loss scaling avoids the purely data-size-driven contraction, results in a much lower-dimensional active subspace and can retain the full non-linear network. Numerical studies demonstrate that the proposed approach yields coherent predictive intervals with only a small number of active directions. 2 Results Let f(x,θ)∈ℝMf(x;θ) ^M represent our Quantity of Interest (QoI), namely the output of a pretrained neural network. Here x are the inputs to the neural network, whereas θ∈ℝDθ ^D is the collection of all connection weights. Furthermore, let (x,y)(x,y) denote all supervised training data, and specific independent identically-distributed training samples are indexed as (xn,yn)(x_n,y_n), n=1,⋯,Nn=1,·s,N, with yn∈ℝMy_n ^M. It is our goal to perform scalable uncertainty quantification of f by propagating a low-rank Bayesian posterior distribution on θ through the network. Before we describe our results, we give a brief overview of our UQ framework, called generalized Laplace active subspaces. Full details are found in the Methods section; the PyTorch code and data to reproduce the results below are available through Edeling (2026). 2.1 Generalized Laplace active subspaces Two choices define our Bayesian posterior. First, instead of seeking a posterior distribution p(θ∣y)p(θ y) in the original D-dimensional coordinate system of the connection weights, we seek a low-rank (d≪Dd D) coordinate system that is most relevant for the posterior predictive uncertainty. Secondly, we do not employ the standard Bayesian belief update p(θ∣y)∝p(y∣θ)p(θ)p(θ y) p(y θ)p(θ), which uses the likelihood function p(y∣θ)p(y θ) to connect the parameters θ to the observed data y. For independent Gaussian observations, the corresponding negative log-likelihood is proportional to a summed loss. In contrast, neural networks are typically trained using an average loss ℒNL_N, such as the mean squared error, which includes an additional factor 1/N1/N. To obtain a Bayesian-style update that is consistent with the average loss used in practice, we adopt the generalized Bayesian framework of Bissiri et al. (2016) instead, which gives rise to the following posterior: πβ(θ∣y)∝exp[−βℒN(θ)]p(θ), _β(θ y) \! [- _N(θ) ]p(θ), (1) see Section 4.1.1 for more details. The so-called ‘temperature’ parameter β controls the contribution of the loss relative to the prior. This formulation includes the standard Bayesian posterior as a special case, obtained by the appropriate choice of β, while also allowing alternative scalings such as the mean-loss scaling considered in this work. It therefore provides a convenient framework for studying how posterior scaling (i.e. the choice of β) affects the dimension and structure of the resulting low-rank posterior approximation. This posterior approximation is obtained by constructing a local Laplace approximation centered at the pre-trained weights θ0∈ℝD _0 ^D. With a commonly-assumed isotropic Gaussian prior p(θ)=(0,σ02ID)p(θ)=N(0,σ^2_0I_D), this results in a Gaussian posterior distribution πβ(θ∣y)=(θ0,Σ) _β(θ y)=N( _0, ), with Σ=(βHL+1σ02ID)−1. = (β H_L+ 1 _0^2I_D )^-1. (2) Here, HLH_L is defined to be the Hessian (curvature) of the mean loss: HL:=∇θ2ℒN(θ0)∈ℝD×DH_L:=∇^2_θL_N( _0) ^D× D. For large neural networks this matrix cannot be formed or inverted explicitly, and it may not be positive semi-definite. We therefore replace the full curvature by a low-rank approximation in the dominant eigenspace of the generalized Gauss-Newton matrix G (defined by (28)), G being a positive semi-definite approximation to HLH_L. The d-dimensional dominant G eigenspace with λ1≥λ2≥⋯≥λd _1≥ _2≥·s≥ _d; Gpj=λjpj,j=1,…,d, Gp_j= _jp_j, j=1,…,d, (3) defines the projection matrix P=[p1,…,pd]∈ℝD×dP=[p_1,…,p_d] ^D× d, used to obtain the posterior perturbations on top of θ0 _0 via θ=θ0+Pz,z∈ℝd. θ= _0+Pz, z ^d. (4) In this way uncertainty is represented only via the dominant curvature-informed directions, along which the data was able to inform the parameters well. In the coordinate system of this subspace the Laplace posterior has diagonal covariance, with σj2=1βλj+σ0−2. _j^2= 1β _j+σ^-2_0. (5) Thus, once the leading eigenpairs (λj,pj)( _j,p_j) of G are available, the posterior variances σj2 _j^2 in the retained active directions are obtained in closed form. The eigenpairs are computed using automatic differentiation, matrix-free vector products GvGv and the Lanczos algorithm (see section 4.1.5). This avoids explicitly forming G, and makes the construction applicable to networks with a large number of parameters. The resulting approximation is therefore both local and low-rank: local because it describes uncertainty around a single trained solution θ0 _0, and low-rank because only a small number of active directions are retained. The prior variance σ02 _0^2 sets the overall scale of the posterior, and can be calibrated by an empirical-Bayes criterion within the active subspace w0=PTθ0∈ℝdw_0=P^T _0 ^d. The resulting prior variance is found by solving ∥w0∥22=∑i=1d1α−1βλi+α=∑i=1dσ02−σi2, w_0 ^2_2= _i=1^d 1α- 1β _i+α= _i=1^dσ^2_0-σ^2_i, (6) for α=σ0−2α=σ^-2_0. Note that this relation chooses the prior scale such that the norm of the projected parameter vector is consistent with the amount of variance reduction induced by the data in the active subspace. The posterior predictive distribution is then approximated by sampling (4) and propagating these samples through the neural network f(x,θ)f(x;θ). This sampling strategy is denoted as sampled Laplace. However, the non-linear nature of the network may distort the UQ estimates. So-called linearized Laplace methods Immer et al. (2021), which replace the full network with its linear expansion (see (45)), are a means to counteract this. We will contrast these two sampling methods below. 2.2 Regression with known observational noise We start with the same scalar (M=1M=1) regression function with Gaussian noise as in Jantre et al. (2024), namely; y(x)=sin(4πx)+sin(7πx)+ϵ,ϵ∼(0,σ2),x∈[0,1], y(x)= (4π x)+ (7π x)+ε, ε (0,σ^2), x∈[0,1], (7) using a known noise of σ=0.1σ=0.1. Implementation verification. This example is used to verify the implementation, in particular the matrix-free computation of GvGv and the dominant eigenpairs of G. We train a fully-connected feed-forward neural network, using back propagation with Stochastic Gradient Descent and the Adam optimizer Kingma and Ba (2014); Aggarwal (2018). The learning rate was set at 0.001, and we used two hidden layers at 50 neurons each. The dimension of the weight space is D=2701D=2701, small enough to explicitly form the full G matrix (28). The matrix-free GvGv vector is contrasted against the corresponding explicit matrix-vector product in Figure 1. The absolute error |Gv−(Gv)explicit||Gv-(Gv)_explicit| is (10−8)O(10^-8) in general (using the torch.float32 data type). A similar verification for the matrix-free eigenpairs is displayed in Figure 2. Figure 1: The matrix-vector product Gv∈ℝDGv ^D with D=2701D=2701, computed using the matrix-free method outlined in Section 4.1.5, and by explicitly forming matrix G via (28). The absolute error |Gv−(Gv)explicit||Gv-(Gv)_explicit| is plotted on the right vertical axis. (a) Eigenvalue fit and errors. (b) Eigenvector fit and errors Figure 2: The d leading eigenvalues λj _j and leading eigenvector p1p_1 of G, computed using the matrix-free method outlined in Section 4.1.5, and by explicitly forming matrix G via (28). The absolute errors are plotted on the right vertical axis. Standard vs generalized scaling. Next we examine the scaling of the output confidence intervals (CIs) with the temperature parameter β, with the CIs computed using the sampled Laplace strategy. As described in Section 4.1.1, the value β=Nσ−2β=Nσ^-2 corresponds to the standard Bayesian setting, whereas with β=σ−2β=σ^-2 we employ the generalized Bayesian scaling where the task of standard likelihood function is taken up by the mean-squared loss. For both posterior scalings we determine the active subspace dimension d as described in Section 4.1.7. Essentially, we increase d∈ℕd , starting from d=1d=1 and solve (6) for σ02σ^2_0 every iteration, until our computed 90% CIs achieve empirical coverage of (at least) 9090% of held-out calibration data. Our goal here is to compare this minimal required subspace dimension d for both posterior scalings. Figure 3 displays the result. The generalized Bayesian posterior only required d=2d=2 active subspace directions to achieve coverage, and produces coherent CIs. Conversely, the standard Bayesian scaling of β required d=31d=31 directions and yields slightly wider, less coherent intervals. (a) β=σ−2β=σ^-2, N=50N=50, d=2d=2 (b) β=Nσ−2β=Nσ^-2, N=50N=50, d=31d=31 Figure 3: 90% confidence intervals for the, a) generalized Bayesian scaling, and b) standard Bayesian scaling. Both CIs were generated by drawing 1000 posterior samples, and propagating these through the network. The preceding results were obtained with N=50N=50 training data points. We now compare the standard and generalized scaling using N=1000N=1000. The results from the generalized Bayesian scaling (not shown) are essentially unchanged by the increase in training data, as it still yields a 2D active subspace that generates consistent CIs. We therefore focus on the 90% CI results of the standard Bayesian scaling, as displayed in Figure 4. Clearly, the quality of the mean prediction and the CIs significantly degraded. Furthermore, instead of d=31d=31, we now required 43 curvature directions to achieve coverage. Figure 4: The results for the standard Bayesian scaling β=Nσ−2β=Nσ^-2 with N=1000N=1000 and d=43d=43. The CIs were generated by drawing 1000 posterior samples, and propagating these through the network. Sampled vs linear Laplace. The underfitted results of Figure 4 are a consequence of the of propagating the Laplace posterior through the non-linear network, see also Immer et al. (2021). Instead, the posterior predictive results obtained from linearized Laplace flinf_lin, i.e. from (45), are shown in Figure 5 for N=100N=100 and N=1000N=1000 for the standard scaling. Due to the linearity of the predictive model, its distribution is Gaussian as well, with a standard deviation σf _f given by (46). We therefore perform the coverage calibration using a 95% confidence interval constructed from f(x,θ0)±2σf(x; _0)± 2 _f, which avoids Monte Carlo sampling. We further imposed a maximum subspace dimension of dmax=50d_max=50. As expected, the posterior predictive mean is improved relative to Figure 4, since the linearized predictive distribution is centered by construction at the expansion-point prediction f(x,θ0)f(x; _0). The confidence intervals are substantially narrower, confirming that the overly wide intervals observed for non-linear posterior sampling are largely caused by propagating weight samples through the non-linear network. However, the linearized predictive distribution is under resolved. For both N=100N=100 and N=1000N=1000, the calibration procedure reaches the maximum allowed dimension of 50, covering only approximately 75% and 40% of the calibration data, respectively, with nominal 95% confidence intervals. (a) N=100N=100, d=50d=50, 75% coverage (b) N=100N=100, d=50d=50, 40% coverage Figure 5: The results for the standard Bayesian scaling β=Nσ−2β=Nσ^-2 combined with the linearized Laplace predictive posterior distribution, computed using f(x,θ0)±2σf(x; _0)± 2 _f. 2.3 Regression with unknown observational noise The following test case involves the Concrete Compressive Strength dataset from the UCI machine learning repository, see Yeh (1998) for details. This is a regression test case containing N=1030N=1030 samples and 8 input features describing the concrete mixture, and a scalar output corresponding to compressive strength. We train the same architecture as the preceding test case, except that now it contains approximately D=106D=10^6 connection weights. This case therefore assesses the scalability of the proposed matrix-free low-rank generalized Laplace approximation. Since the observation noise is unknown, we estimate the noise σ via residuals: σ2≈1N∑n=1N∥yn−f(xn,θ0)∥22, σ^2≈ 1N _n=1^N y_n-f(x_n; _0) _2^2, which is approximately 0.07 for this test case, and we use this value to specify β as before. We set the number of training data N to be 90% of the 1030 available data points. The remainder (NcalN_cal) is used to calibrate d. We execute the calibration procedure until our target coverage of 90% is met, or until a preset maximum active subspace dimension dmax=50d_max=50 is exceeded. Standard vs generalized scaling. The coverage for both the generalized and standard Bayesian scaling is plotted versus d in Figure 6(a). In the case of the former, the target coverage was reached at d=3d=3, while for the latter it was never reached within the allotted range. In fact, even at d=50d=50 only approximately 60% of the calibration data was covered by the 90% CIs. (a) Sampled Laplace posterior prediction, 90% nominal coverage. (b) Linear Laplace posterior prediction, 95% nominal coverage. Figure 6: The coverage of the calibration procedure as outlined in Section 4.1.7, versus the active subspace dimension d, for both the standard Bayesian scaling β=Nσ−2β=Nσ^-2 and the generalized scaling β=σ−2β=σ^-2. Sampled vs linear Laplace. The same coverage study was conducted with the linear predictive posterior distribution, using again a 95% confidence interval constructed from f(x,θ0)±2σf(x; _0)± 2 _f. The overall results, see Figure 6(b), are similar compared to the sampled Laplace study. For the generalized scaling, with d=4d=4 active dimensions the coverage criterion is met, while the coverage of the standard scaling stagnates. Due to the multivariate input x, we display the generalized Bayesian posterior predictive results in Figure 7(a) in the form of a parity plot. Here we show the mean predicted compressive strength versus itself with 90% sampled Laplace CIs, and the mean predicted compressive strength versus the calibration data. Visually speaking, for the generalized scaling the CIs are consistent, as the mean is generally centered around the validation data, and the 90% CIs encompass most data without being overly wide. The CIs generated by the standard scaling (Figure 7(b)) are indeed overconfident. (a) Generalized mean-loss scaling, d=3d=3. (b) Standard Bayes scaling, d=50d=50. Figure 7: The parity plot of the mean compressive strength of concrete, with 90% CIs. We show the standard and generalized scaling results, both computed with sampled Laplace. Activity scores. The dominant eigenvectors of G define important directions in parameter space. However, we can also extract information about the importance of individual connections weights. Activity scores, originally derived in Constantine and Diaz (2017), are (first-order) sensitivity indices for individual parameters, which can be extracted from the active subspace eigenmodes. Analogously, using the G eigenmodes, we define the activity scores for individual connection weights as: νi=∑j=1dλjpj,i2,i=1,⋯,D. _i= _j=1^d _jp_j,i^2, i=1,·s,D. (8) Here pj,ip_j,i is the i-th component of the j-th eigenvector of G. The activity score νi _i measures how strongly parameter coordinate i contributes to the retained d-dimensional curvature-informed active subspace, weighted by the importance λj _j of each active direction. Since the connection weights have no physical interpretation, we will use (8) primarily to examine if the contribution to the curvature is broadly distributed across the network, or concentrated in a small subset of connection weights. However, rather than plotting individual νi _i, we first sort from largest to smallest and instead plot the cumulative scores κi=∑j=1iνj∑j=1Dνj,i=1,⋯,D, _i= _j=1^i _j _j=1^D _j, i=1,·s,D, (9) in a log-log plot. The results are shown in Figure 8. We see that the 1000 most-sensitive connection weights (just 0.1% of the total) rapidly saturate and are responsible for most of the curvature sensitivity. In particular, for this case the remaining 99.9% of weights account for 5% of the total sensitivity. This indicates that, in addition to being only d=3d=3 dimensional, the active subspace is also localised in the high-dimensional weight space. Figure 8: The cumulative (normalized) activity scores for our concrete compression strength feed-forward neural network with approximately D=106D=10^6 weights. 2.4 Autoregressive classification for molecular generative modelling For our final (classification) test case we consider the autoregressive molecular language model from Olivecrona et al. (2017). Here, molecules are represented as simplified molecular-input line-entry specification (SMILES) strings, a sequence of characters corresponding to atoms, generated sequentially by a recurrent neural network trained to predict a discrete distribution for the next token. In particular, we use a 3-layer LSTM network part of the REINVENT tookit (Loeffler et al. (2024)) with 256 neurons per layer and approximately 1.6M weights, trained on SMILES strings from the ChEMBL database (Gaulton et al. (2017)). The M=34M=34 tokens that make up the SMILES vocabulary and some example ChEMBL molecules can be found in the Appendix. This makes molecular generation a sequence of autoregressive token-level classification problems since, at each step, the network outputs a softmax distribution over the vocabulary, from which the next token is sampled. We study here how the low-rank posterior uncertainty in the network (for both the standard and generalized scaling) affects free-running autoregressive molecular generation. This setting separates the intrinsic stochasticity of the softmax layer from epistemic variability induced by posterior weight perturbations. Reducing computational cost via subsampling. The LSTM was trained on 1.251.25M SMILES strings, each containing a (varying) number of tokens, which results in a total number of data records of N=58.7N=58.7M. As each GvGv product used in the Lanczos algorithm is an empirical average over the dataset (see (30)), we choose to randomly subsample the data for computational efficiency, computing the eigendecomposition of G using NsubN_sub SMILES strings instead. Note a similar strategy was employed by Faller and Martin (2026). To examine the effect of this subsampling we plot the (absolute value of) the inner product between the dominant eigenvector p1p_1 and a reference vector p1refp_1^ref as a function of NsubN_sub in Figure 9. This reference vector was computed using Nsub=10kN_sub=10k, and we replicate this procedure Nrep=5N_rep=5 times at each subsampling to account for sampling noise, creating 25 possible inner products per NsubN_sub. While clear variation is present in the replica inner products for Nsub=500N_sub=500 and Nsub=1000N_sub=1000, |⟨p1,p1ref⟩| | <p_1,p_1^ref > | soon clusters near 1, indicating that reasonable alignment with the reference eigenvector can be obtained with a small fraction of the total dataset. That said, we note that there is an obvious outlier for Nsub=3000N_sub=3000. This is most-likely caused by a p1p_1 which failed to converge during the Lanczos algorithm, such that all 5 corresponding |⟨p1,p1ref⟩| | <p_1,p_1^ref > | cluster near 0.14 in this case. One simple way to flag such outliers for a given NsubN_sub is computing the Nrep×NrepN_rep× N_rep matrix with entries |⟨p1,i,p1,j⟩| | <p_1,i,p_1,j > |, i,j=1,⋯,Nrepi,j=1,·s,N_rep, p1,ip_1,i being the i-th replica eigenvector. Ill-converged vectors will cause entries away from 1. All subsequent analysis was performed using Nsub=2000N_sub=2000 SMILES, corresponding to N=93kN=93k data points. Figure 9: The inner products |⟨p1,p1ref⟩| | <p_1,p_1^ref > | as a function of Nsub∈500,1000,1500,2000,2500,3000N_sub∈\500,1000,1500,2000,2500,3000\, replicated 25 times per NsubN_sub. The reference vectors p1refp_1^ref were computed using Nsub=10.000N_sub=10.000 SMILES strings. Standard vs generalized scaling. Our first results concern the output validity, as not all generated molecules are guaranteed to be chemically correct compounds. REINVENT verifies the validity in a post processing step using the RDKit Landrum et al. (2026). Note that validity here means RDKit can parse the SMILES string into a molecular graph and the graph passes basic chemical sanitization checks. It does not imply drug-likeness, synthetic accessibility, or chemical plausibility. Table 1 shows the number of valid molecules out of a set of 10k generated SMILES strings for the baseline REINVENT model (i.e. with softmax-only uncertainty and weights fixed at θ0 _0), and the perturbed model Laplace models with standard and generalized scaling. When varying the active-subspace dimension from 1 to 3, the standard scaling does not result in a notable change in validity compared to the baseline model. For the generalized scaling, validity notably drops. Table 1: Number of valid SMILES strings (out of 10k). baseline Laplace (std. scaling) Laplace (gen. scaling) d=1d=1 8692 8709 6863 d=2d=2 8692 8753 5279 d=3d=3 8692 8732 4435 The following results were obtained using only valid SMILES strings. Figure 10 shows the kernel density estimates of the sequence length of the molecules (i.e. the number of tokens), computed from 10k generated SMILES strings. Once again only the generalized scaling shows a departure from the baseline results. Especially for d=3d=3 the generator produces longer SMILES strings, with an increased probability that the maximum allowed sequence length of 256 is reached. Finally, we examine the effect of the Laplace active subspaces on the distribution of several Quantities of Interest (QoIs), described in Table 2. The corresponding distribution are found in Figure 11. Figure 10: Kernel density estimates of the sequence length for the standard and generalized scaling with d∈1,2,3d∈\1,2,3\. The maximum allowed sequence length is 256. Table 2: Molecular quantities of interest used to compare generated SMILES distributions within REINVENT. Quantity Description Molecular weight Molecular size based on atomic masses, large values indicate large molecules. SA score Synthetic Accessibility score, higher values indicate more difficult synthesis. SlogP Estimate of how fat-soluble a molecule is. (a) Molecular weight distributions. (b) SA score distributions. (c) SlogP distributions. Figure 11: Kernel density estimates of the QoIs for the standard and generalized scaling with d∈1,2,3d∈\1,2,3\. Activity scores. The cumulative activity scores κi _i are shown in Figure 12, using d=2d=2. Note the curve saturates less quickly compared to the concrete compression strength case. Still, the first 7.6% of the activity scores already account for 95% of the total sensitivity. Figure 12: The cumulative (normalized) activity scores for the REINVENT LSTM with approximately 1.6M weights. Computed using d=2d=2. 3 Discussion We now discuss our findings and contrast our proposed method with related work. 3.1 Posterior scaling and active-subspace dimension The results show that the scaling of the generalized posterior plays a central role in the quality and efficiency of our low-rank approximation. Under the standard Bayesian scaling, the leading active directions become increasingly contracted with the number of training data, as σj2=(βλj+σ0−2)−1σ^2_j=(β _j+σ^-2_0)^-1 with β=Nσ−2β=Nσ^-2 or β=Nβ=N for regression or classification respectively. This forces coverage calibration to retain many weak-curvature directions. Note that this does not mean the standard scaling is intrinsically defective, since posterior contraction with increasing data is an expected feature of the standard Bayesian update. The standard Bayesian scaling and generalized mean-loss scaling define do different posterior objects, and in our specific context of low-rank posterior UQ for overparameterized neural networks, the generalized scaling can provide a more effective posterior approximation. To further discuss this behaviour, consider Figure 13. For the analytic regression function (7), we plot the posterior variances and the contraction ratio; ci:=σ02−σi2σ02=βλiβλi+1σ02, c_i:= σ^2_0-σ^2_iσ^2_0= β _iβ _i+ 1σ^2_0, (10) for which ci≈1c_i≈ 1 implies that the i-th direction is strongly data-informed (since then σi2→0σ^2_i→ 0), and values close to zero denote prior-dominated directions (σi2→σ02σ^2_i→σ^2_0). Under the standard Bayesian scaling (Figures 13(b)-13(d)), while the leading directions are almost fully contracted (posterior dominated), the posterior variances remain close to zero over much of the leading spectrum, due to the scaling with N: σi2≈1/(Nσ−2λi)σ^2_i≈ 1/(Nσ^-2 _i) when ci≈1c_i≈ 1. Coverage calibration therefore selects directions much deeper in the spectrum, where the posterior variances increase rapidly, as the σi2σ^2_i span multiple orders of magnitude. Consequently, nominal coverage is obtained only after adding directions that introduce substantial parameter-space variability along weak-curvature modes. This effect becomes more pronounced with more training data; in Figure 13(d) the variance along the leading active directions is reduced by an order of magnitude compared to Figure 13(b). As a result, d=43d=43 was required to meet the coverage criterion, which exacerbates the aforementioned problems of incorporating weak curvature modes with high variance. In the non-linear network, these larger perturbations may move samples outside the local region in which the Laplace approximation is accurate, resulting in wider and less stable predictive intervals, as shown in Figure 4. A better mean posterior prediction can be obtained with a lower d, at the cost of overconfident CIs, or by utilizing the linear predictive distribution. Such overconfident CIs were also observed with the UCI concrete compressive strength network with standard scaling, where 50 curvature directions were insufficient to achieve coverage. The slow increase displayed in Figure 6 also suggests a significantly higher subspace dimension is required to achieve target coverage. Moreover, the retained eigenvalue spectrum spans several orders of magnitude, indicating that the additional directions correspond to genuinely low-curvature modes. While the linear predictive distribution does remove the non-linear sampling artifacts when they occur, the excessive posterior contraction with N induced by the standard Bayesian scaling is still present for linear Laplace. This is clear from Figures 5-7. In contrast, under the generalized Bayesian scaling, sufficient uncertainty is retained in the leading active directions, regardless of N or the posterior predictive method, and desired coverage was obtained with only d=4d=4 or less in all cases. Thus, the uncertainty is concentrated in a small number of coherent active directions rather than spread across many weak-curvature modes. (a) β=σ−2β=σ^-2, N=50N=50 (b) β=Nσ−2β=Nσ^-2, N=50N=50 (c) β=σ−2β=σ^-2, N=1000N=1000 (d) β=Nσ−2β=Nσ^-2, N=1000N=1000 Figure 13: The contraction ratio (10) and posterior variances (22) for the generalized (left column) and standard (right column) Bayesian scaling, for both N=50N=50 and N=1000N=1000. These results were computed using the analytic regression function (7). The increase in the selected active subspace dimension under the standard Bayesian scaling is also computationally consequential. Although N=1000N=1000 is not large by modern machine-learning standards, the standard scaling already requires a substantially higher posterior rank to achieve the desired coverage compared to N=50N=50, which is a problem that will escalate with larger N. That said, as we showed in Section 2.4, for large datasets one could subsample the data when computing the eigenmodes of G. Still, computing these modes under the standard scaling will require a larger d, and therefore more Lanczos iterations (see Section 4.1.5). Thus, the standard Bayesian scaling not only shifts uncertainty into weaker-curvature directions, but also increases the computational cost of the low-rank posterior approximation compared to the generalized scaling. In contrast, the generalized Bayesian mean-loss scaling achieves calibrated and sharp confidence intervals with a much smaller active subspace. The proposed approach therefore provides a practical route to calibrated predictive uncertainty for regression-based neural networks. Similarly, for our autoregressive classification case, the standard scaling resulted in distributions that do not deviate from the baseline REINVENT model in any notable fashion, see Figure 11. The low-rank Laplace perturbations with d=1,2d=1,2 appear to add meaningful variation to the generated molecular distribution while remaining broadly within the regime learned by the pretrained generator. In contrast, d=3d=3 clearly crosses a stability boundary. The generator can start to produce some long (carbon-rich), formally valid yet chemically implausible molecules with extreme SlogP values. Lipinski’s Rule of Five Lipinski et al. (1997) is a rule of thumb used in drug discovery to evaluate if a chemical compound has the right properties to make it an orally active drug in humans. One violating condition of this rule is an SlogP value greater than 5, which occurs at a significantly higher probability once d=3d=3. This suggests that the addition of third active direction pushes the autoregressive dynamics outside the chemically meaningful region of the model. Overall, under standard Bayesian scaling, the low-rank Laplace generator remains close to the pretrained baseline for all tested subspace dimensions. Valid-molecule QoI distributions and sequence lengths are not meaningfully affected, indicating that the posterior samples remain in a local generative regime. The reason behind this is the same as in the regression case, and again will be become even more pronounced with increasing (subsampled) data size. In contrast, under generalized scaling, increasing the active-subspace dimension produces progressively stronger deviations from the baseline. The most pronounced effect is in the sequence-length distribution, where a higher d introduces a long right tail, suggesting that generalized posterior perturbations affect the learned stopping behaviour. The same scaling also broadens the QoI distributions, indicating increased chemical exploration at the cost of generative stability. However, this was only a sensitivity study of the effect of posterior weight perturbations on free running molecule generation. Our pretrained REINVENT model is already reasonably calibrated at the token level, and we therefore did not expect the additional uncertainty of the Laplace posterior to improve the predictive distribution. That said, low-rank Laplace methods have also been proposed for fine-tuned large language models (LLMs). For instance, the low-rank adaptation (LoRA) method of Hu et al. (2021) is a means to fine tune an LLM by introducing low-rank updates to selected weight matrices. This deterministic fine-tuning can make the model poorly calibrated, often overconfident. Laplace-LoRA, introduced by Adam. et al. (2024) places a Laplace posterior over the fine-tuned LoRA parameters to induce variability in the logits/predicted probabilities. The prior precision is tuned using gradient descent, maximizing the Laplace evidence, such that the posterior predictive distribution becomes better calibrated and achieves lower negative log likelihood. Extending our generalized-Bayes active-subspace Laplace approximation to fine-tuned language models is therefore a promising direction for future work. This would shift the role of the Laplace approximation from probing generative sensitivity, as in the present REINVENT experiments, toward correcting overconfidence in fine-tuned autoregressive models. Finally, the activity scores from Figures 8 and 12 show that the active subspace is concentrated in a very small subset of all connection weights, relative to the total weight space. However, the active subspace w=PTθw=P^Tθ itself provides an even more compact representation, as it is a linear combination of connection weights instead of a collection of individual weights. Similar observations have been made in physics-based models, see e.g. Edeling et al. (2024). 3.2 Related work In recent years several contributions to the problem of uncertainty quantification for neural networks have been made, see Psaros et al. (2023) or Wang et al. (2026) for survey articles. Here we discuss those most relevant to our work. Our work is related to the subspace inference framework of Izmailov et al. (2020), where Bayesian model averaging is performed within low-dimensional parameter subspaces constructed from random projections, SGD-trajectory PCA, or mode-connecting curves. In contrast, we construct the subspace from the dominant generalized Gauss–Newton curvature directions of the empirical mean loss. In Izmailov et al. (2020) posterior tempering is also used as a calibration device, although here the temperature is a hyper parameter selected by cross-validation, rather than being fixed by the loss function used in generalized Bayes. Our curvature-based setting exposes the different scaling issue of increasingly contracted leading active directions with N for the standard scaling. The generalized Bayesian mean-loss scaling avoids this purely data-size-driven contraction. A further distinction is that our generalized Laplace approximation yields a closed-form Gaussian posterior in the active subspace, avoiding Markov-Chain Monte Carlo sampling or variational optimization within the subspace. Most closely related is the low-rank Laplace subspace construction of Faller and Martin (2026), whose authors derive an optimal low-dimensional subspace model for approximating the predictive output covariance of a full Laplace approximation and propose a scalable compute strategy. Their construction is formulated in the linearized Laplace setting discussed previously. This linearization yields a closed-form Gaussian on both the networks weights and the output, which Faller and Martin (2026) use to find the subspace projection matrix P by minimizing the discrepancy between the output covariance of the subspace model and that of the full Laplace approximation. However, the exact construction requires access to quantities that are generally inaccessible for large neural networks, including the full (GGN) curvature matrix and an eigendecomposition of the full output covariance. In practice, scalable variants therefore rely on additional approximations such as diagonal or Kronecker-factored curvature estimates Grosse and Martens (2016). Our goal is different. Rather than computing the output covariance of the full Laplace approximation, we construct a generalized-Bayes Laplace approximation directly in a curvature informed active subspace. The posterior samples are propagated through the non-linear network itself rather than through a linearized predictive model. This distinction is important since linearized Laplace methods keep the posterior predictive mean anchored at the MAP prediction (see section 4.1.8), and can thereby avoid the underfitting behaviour observed when non-linear weight samples are propagated through the network. In our non-linear predictive setting, the size and coherence of the sampled parameter perturbations are therefore central; if the calibrated posterior requires many weak-curvature directions, the resulting samples may leave the local region in which the Laplace approximation is reliable. The generalized Bayesian scaling mitigates this effect by achieving calibrated coverage with a much smaller active subspace, thereby reducing both non-linear sampling artifacts and the number of curvature-vector products required. Also, the authors of Faller and Martin (2026) shows that a low-dimensional subspace can approximate the full predictive covariance of the linearized Laplace approximation, but the required subspace dimensions are still in the tens to hundreds in their experiments. While this is tiny compared to the size of the full weight space, our generalized-Bayes active subspace often achieves calibrated coverage with only a handful of curvature directions, preserving the intended low-rank character of the approximation. There are a number of other contributions that use the Laplace approximation of the posterior. The authors of Ritter et al. (2018) proposed a scalable Laplace approximation for neural networks using the aforementioned Kronecker-factored curvature estimates of Grosse and Martens (2016). This approach improves scalability by imposing a layer-wise Kronecker structure on the curvature, whereas our method obtains scalability by restricting the posterior to a low-rank generalized Gauss-Newton active subspace. Consequently, our approach does not require layer-wise Kronecker independence assumptions. The broader practical relevance of Laplace approximations for Bayesian deep learning has been demonstrated by Daxberger et al. (2021), whose authors review and benchmark several scalable Laplace variants and provide a dedicated laplace PyTorch library. Related active-subspace approaches based on variational inference include Jantre et al. (2024) and Nafiz et al. (2024). These methods draw on the classical active-subspace construction, in which the subspace is identified from an uncentered covariance of model gradients Constantine et al. (2014). This matrix is typically averaged over a prescribed prior weight distribution. In contrast, our active subspace can be viewed as a Monte Carlo approximation of a data-averaged local curvature matrix. Specifically, the generalized Gauss–Newton matrix is an empirical average over the training data, evaluated at the pretrained weights. Thus, unlike prior-averaged active-subspace constructions, our subspace is tied directly to the local data-fit geometry that also defines the generalized Laplace approximation. 4 Methods We now describe the methodological details of the proposed low-rank generalized Laplace framework for uncertainty quantification in neural networks. 4.1 Statistics Our method starts from a generalized Bayesian posterior defined through the empirical mean loss, rather than the summed negative log-likelihood used in standard Bayesian analysis. Around a pretrained set of weights, we then construct a Laplace approximation whose covariance is determined by the local curvature of the loss. Since the full Hessian is intractable and may be indefinite, we replace it by a positive semi-definite generalized Gauss-Newton approximation Kunstner et al. (2019) and compute its dominant eigenspace using matrix-free products computed with automatic differentiation. The resulting Gaussian posterior is restricted to this low-dimensional active subspace, with the prior variance calibrated by an empirical generalized-Bayes evidence criterion. We end with a discussion on posterior sampling strategies. 4.1.1 Generalized Bayes formulation To define a posterior update, a mechanism to connect the parameters to the data is required. In a generalized Bayesian setting this is achieved via a loss function rather than through a traditional likelihood function. This allows us to construct a generalized posterior based on the empirical mean loss while retaining the familiar structure of a prior-to-posterior belief update. We define the mean loss function as ℒN(θ,y)=1N∑n=1Nln(θ,yn). _N(θ,y)= 1N _n=1^Nl_n(θ,y_n). (11) We will consider: ln=12∑m=1M(yn(m)−fn(m))2andln=−∑m=1Mym(m)logpn(m), l_n= 12 _m=1^M (y_n^(m)-f^(m)_n )^2 l_n=- _m=1^My^(m)_m p_n^(m), (12) i.e. squared-loss and cross-entropy loss terms for regression and classification respectively. In the former yn∈ℝMy_n ^M are the (noisy) observations of the true process, whereas in the latter yny_n is the M-dimensional one-hot encoded vector with yn(m)=1y_n^(m)=1 for m=cm=c and yn(m)=0y_n^(m)=0 for m≠cm≠ c, c being the index of the correct class. Here, the neural-network outputs fn(m):=f(m)(xn,θ)f_n^(m):=f^(m)(x_n;θ) are the logits that determine the softmax probabilities pn(m)=exp(fn(m))∑jexp(fn(j)). p_n^(m)= (f_n^(m)) _j (f_n^(j)). (13) The generalized Bayesian (unnormalized) posterior is given by Bissiri et al. (2016) as: π~β(θ∣y)=exp[−βℒN(θ,y)]p(θ), π_β(θ y)= [- _N(θ,y) ]p(θ), (14) where we assume the following prior: p(θ)=(0,σ02ID). p(θ)=N (0,σ^2_0I_D ). (15) Here, IDI_D is the D-dimensional identity matrix and β>0β>0 is a so-called temperature parameter, which typically is a hyper parameter used to control how random the output of a neural networks is. Here we will use it to distinguish between generalized and standard Bayesian setting. Further note that the prior variance of the connection weights σ02∈ℝσ^2_0 will play an important role in subsequent analysis. Finally, we write the normalized posterior as πβ=π~β/Zβ _β= π_β/Z_β, with the normalizing constant ZβZ_β denoted as the marginal likelihood or the evidence. In the case of a regression problem with Gaussian noise, the following statistical model is often assumed to relate the model f to the data y; yn=f(xn,θ)+ϵn,ϵn∼(0,σ2IM). y_n=f(x_n;θ)+ _n, _n (0,σ^2I_M ). (16) Here, σ2σ^2 is the noise variance, which is assumed known, or estimated from residuals. Under this model, the standard Bayesian posterior is recovered by setting β=Nσ−2β=Nσ^-2 in (14), which is equivalent to replacing ℒNL_N with a summed loss. In contrast, the generalized posterior based on the mean squared loss is obtained using β=σ−2β=σ^-2. For classification with softmax cross-entropy, there is no analogous Gaussian noise variance. Similarly, β=Nβ=N corresponds to the summed categorical negative log-likelihood of the standard Bayesian posterior, whereas the generalized posterior based on the mean cross-entropy loss requires β=1β=1. 4.1.2 Laplace approximation of the generalized posterior We define lβ(θ)l_β(θ) as the negative log of (14), i.e.; lβ(θ):=−logπ~β=βℒN−logp(θ)∝βℒN+12σ02∥θ∥22. l_β(θ):=- π_β= _N- p(θ) _N+ 12σ^2_0 θ ^2_2. (17) The most-likely maximum a-posteriori (MAP) parameter value solves θ^=argminθlβ(θ). θ= _θl_β(θ). (18) We will use the pretrained weights θ0 _0 in place of the MAP value, which is a common approximation (also used by e.g. Ritter et al. (2018) and Immer et al. (2021)) unless the network was trained using L2 weight decay with regularization strength consistent with the values of 1/(2σ02)1/(2σ^2_0) and β. In our approach, θ0 _0 serves as the anchor point for the Laplace approximation, originally applied by MacKay (1992) to small neural networks. This approximation models the posterior distribution πβ(θ∣y) _β(θ y) as a Gaussian distribution, with a covariance structure defined from (the inverse of) the curvature of the log-likelihood function. The underpinning idea here is that the log likelihood quickly drops away (i.e. has high curvature) when the data is able to inform the parameters well. For generalized Bayes, the covariance structure is derived from the curvature of the loss function instead. To arrive at such an expression a second-order Taylor expansion of lβ(θ)l_β(θ) around θ0 _0 is made lβ(θ)≈lβ(θ0)+12(θ−θ0)T∇θ2lβ(θ0)(θ−θ0), l_β(θ)≈ l_β( _0)+ 12 (θ- _0 )^T∇^2_θl_β( _0) (θ- _0 ), (19) where under the assumption that θ0 _0 is sufficiently close to a stationary point of lβl_β, the first-order term drops out: ∇θlβ(θ0)=0 _θl_β( _0)=0. Given (19), the posterior πβ _β now becomes a Gaussian (θ0,Σ)N( _0, ): πβ(θ∣y)∝exp[−lβ(θ)]∝exp[−12(θ−θ0)TΣ−1(θ−θ0)], _β(θ y) [-l_β(θ) ] [- 12 (θ- _0 )^T ^-1 (θ- _0 ) ], (20) with Σ=(∇θ2lβ(θ0))−1=(β∇θ2ℒN(θ0)−∇θ2logp(θ))−1=(βHL+1σ02ID)−1. = (∇^2_θl_β( _0) )^-1= (β∇^2_θL_N( _0)-∇^2_θ\,\; p(θ) )^-1= (β H_L+ 1 _0^2I_D )^-1. (21) Note that here HLH_L is defined to be the Hessian of the mean loss, HL:=∇θ2ℒN(θ0)∈ℝD×DH_L:=∇^2_θL_N( _0) ^D× D. The Gaussian form is thus a consequence of truncating the Taylor expansion around θ0 _0 at the quadratic term, it is not statement that the full neural-network posterior is indeed normally distributed, which may well be multi modal and non-Gaussian. The Laplace approximation is therefore a pragmatic, local posterior estimate around one trained solution, not as a faithful global description of all plausible weights. We examine its usefulness through its posterior predictive behaviour, i.e. by the resulting posterior distribution of the neural-network output. The Laplace approximation gives a closed-form expression for the posterior on θ. Unlike UQ approaches based on (low-rank) variational inference (e.g. Jantre et al. (2024); Nafiz et al. (2024)), no additional training is therefore required. Still, computing and inverting HLH_L is computationally intractable in neural networks due to the size of θ. To overcome this limitation we will construct a computationally tractable closed-form expression for the posterior variance. 4.1.3 Low-rank posterior variance By examining (21), we observe that in the principal coordinate system spanned by the eigenvectors from HLpj=λjpjH_Lp_j= _jp_j, the posterior variances σj2 _j^2 become; σj2=1βλj+σ0−2,j=1⋯,D. _j^2= 1β _j+ _0^-2, j=1·s,D. (22) In practice, a low-rank framework is obtained by retaining only the first d≪Dd D dominant eigenvalues λj _j of an approximation to the Hessian. The quality of this approximation depends on how well the approximate Hessian captures directions of significant curvature (see Section 4.1.4), as well as on the choice of prior variance σ02 _0^2. It is not uncommon to find references which discuss Hessian approximations in the context of the sum loss, see e.g. Kunstner et al. (2019). As discussed, in the present work we use the loss from generalized Bayes rather than the summed negative log-likelihood. This corresponds to a temperature-scaled posterior where the parameter β controls how strongly the data loss is weighted relative to the prior (see (14)). For β=σ−2β=σ^-2 or β=1β=1 (mean regression / classification loss), the likelihood contribution is normalized by the number of data points. This normalization prevents the curvature eigenvalues from scaling linearly with N, such that the curvature is governed by the shape of the mean loss landscape, not by the size of the dataset. This is important in our low-rank subspace setting, since as N increases, the posterior variances (22) decrease in the standard Bayesian setting with β=Nσ−2β=Nσ^-2 or β=Nβ=N. While in this setting posterior contraction is expected as N grows, in a low-rank projected posterior this can become too aggressive, and the leading active directions become extremely narrow. Using the mean loss avoids this purely sample-size-driven contraction. It makes the eigenvalues describe the average local sensitivity of the trained model rather than the total accumulated data-fit curvature. Regardless of the choice of β, the form σj2=(βλj+1/σ02)−1 _j^2=(β _j+1/ _0^2)^-1 highlights the interplay between data and prior. Directions with large average curvature (large λj _j) lead to stronger posterior contraction, while in flat directions (λj→0 _j→ 0) the posterior variance reverts to the prior, σj2→σ02 _j^2→ _0^2. Thus, the prior variance acts as a baseline level of uncertainty in directions not informed by the data. Since σ02 _0^2 appears explicitly in the posterior variance and controls the overall scale of uncertainty, it must be chosen carefully. We estimate σ02 _0^2 using an empirical Bayes procedure in Section 4.1.6. First, we discuss common choices for the approximate Hessian, and means to efficiently extract their d dominant eigenpairs. 4.1.4 Approximate Hessians In general, an issue with using the true Hessian in the covariance expression is that it is not guaranteed to be positive semi-definite. Several approximations are therefore commonly employed, most notably the generalized Gauss–Newton matrix G, which coincides with the true Hessian only under specific assumptions Kunstner et al. (2019). Understanding the precise conditions under which these matrices agree, and how they differ otherwise, is essential for justifying their use as curvature surrogates. In the following, we briefly make these relationships explicit and discuss their implications for subspace-based uncertainty quantification. One can decompose the Hessian into a term involving second derivatives of the model output, and a (more tractable) term that depends only on first-order derivatives. A natural approach is to apply the chain rule to (11), differentiating with respect to f before taking the θ derivative. The resulting gradient will have entries (∇θℒN)i=1N∑n=1N∑m=1Mrn(m)∂fn(m)∂θi,i=1,⋯,D, ( _θL_N )_i= 1N _n=1^N _m=1^Mr_n^(m) ∂ f_n^(m)∂ _i, i=1,·s,D, (23) where rn(m)r_n^(m) are residuals, defined as: rn(m):=fn(m)−yn(m)orrn(m):=pn(m)−yn(m), r_n^(m):=f_n^(m)-y_n^(m) r_n^(m):=p_n^(m)-y^(m)_n, (24) for regression or classification. Note that for the latter, rn(m)r_n^(m) is a probability residual, which unlike its regression counterpart, is a function of all M fn(m)f_n^(m) entries (see (13)). Applying the chain rule to (23) once more gives the Hessian entries: (∇θ2ℒN)i,j=1N∑n∑m(rn(m)∂2fn(m)∂θiθj+∑k∂rn(m)∂fn(k)∂fn(k)∂θi∂fn(m)∂θj),i,j=1,⋯,D. (∇^2_θL_N )_i,j= 1N _n _m (r_n^(m) ∂^2f_n^(m)∂ _i _j+ _k ∂ r_n^(m)∂ f_n^(k) ∂ f^(k)_n∂ _i ∂ f^(m)_n∂ _j ), i,j=1,·s,D. (25) When written in matrix form we obtain: HL(θ):=∇θ2ℒN=1N∑n∑mrn(m)∇θ2fn(m)+1N∑n(∇θfn)THln(∇θfn). H_L(θ):=∇^2_θL_N= 1N _n _mr_n^(m)∇^2_θf_n^(m)+ 1N _n ( _θf_n )^TH_l_n ( _θf_n ). (26) Here, ∇θfn∈ℝM×D _θf_n ^M× D is the Jacobian of the output, and Hln:=∇f2ln∈ℝM×MH_l_n:=∇^2_fl_n ^M× M is a Hessian with (Hln)m,k=∂2ln/∂fn(m)∂fn(k)=∂rn(m)/∂fn(k) (H_l_n )_m,k=∂^2l_n/∂ f_n^(m)∂ f^(k)_n=∂ r_n^(m)/∂ f_n^(k), leading to Hln=IMandHln=diag(pn)−pnpnT H_l_n=I_M H_l_n=diag(p_n)-p_np_n^T (27) for regression and classification respectively. The generalized Gauss–Newton (GGN) matrix G is an approximation to HLH_L that captures (partial) curvature information while avoiding second-order derivatives of the model, and is obtained by simply ignoring the first term of (26); HL(θ)≈G(θ):=1N∑n(∇θfn)THln(∇θfn). H_L(θ)≈ G(θ):= 1N _n ( _θf_n )^TH_l_n ( _θf_n ). (28) This is convenient as it drops the (expensive) Hessians of fn(m)f_n^(m). Moreover, from (26) we can see that this is a well-justified approximation as argued by Kunstner et al. (2019), in the sense that the neglected term vanishes when the residuals are small. For (28) to be exact, it would be sufficient if all residuals are zero. However, this is unrealistic, and would be indicative of overfitting. The authors of Immer et al. (2021) therefore prefer an alternative justification for the use of G, namely that it is equivalent to replacing the network by its first-order Taylor expansion flinf_lin around θ0 _0, for which ∇θ2flin(θ0)=0∇^2_θf_lin( _0)=0 holds. This leads to so-called linearized Laplace predictive distributions, discussed in Section 4.1.8. Finally, we note that empirical work has shown that the loss gradient ∇θℒN _θL_N aligns strongly with the leading eigenvectors of the Hessian, see Gur-Ari et al. (2018), implying that optimization proceeds primarily within a subspace spanned by a small number of dominant curvature directions. This observation provides support for low-rank posterior approximations involving (approximate) Hessians; if training mainly explores a restricted subspace, then uncertainty quantification should focus on these same directions. Our approach leverages this structure by constructing a posterior covariance aligned with the dominant eigenspace of G(θ)G(θ). As mentioned, we use the term active subspace to denote the dominant eigenspace of the generalized Gauss-Newton curvature of the empirical mean loss, which is a slight departure from its original definition Constantine et al. (2014). That said, the overarching idea is similar, as is the linear projection to the active subspace; w0=PTθ0∈ℝd, w_0=P^T _0 ^d, (29) here applied to the connection weights. P∈ℝD×dP ^D× d contains the d≪Dd D dominant (orthonormal) eigenvectors of G. In order to form P for large neural networks, we first discuss matrix-free eigenvalue methods. 4.1.5 Matrix-free computation of dominant eigenpairs In practice, computing the dominant eigenvalues / vectors of G∈ℝD×DG ^D× D will be complicated by the dimension of G. For large enough D, even storing G in memory will not be possible. However, we can use the Lanczos algorithm Chen (2026); Lanczos (1950) to approximate the leading eigenpairs, which only requires access to matrix-vector products GvGv. These products can be computed using automatic differentiation without ever forming G. If we write Jn:=∇θfn∈ℝM×DJ_n:= _θf_n ^M× D for brevity, we obtain Gv=1N∑nJnTHlnJnv, Gv= 1N _nJ_n^TH_l_nJ_nv, (30) see (28). In general, by the chain rule, the Jacobian-vector product Jnv∈ℝM×1J_nv ^M× 1 equals Jnv=dϵf(xn,θ+ϵv)|ϵ=0, J_nv= . ddεf(x_n;θ+ε v) |_ε=0, which in turn can be computed with forward-mode automatic differentiation, see e.g. Pearlmutter (1994). Once JnvJ_nv is available at the output, pre-multiplication with the known HlnH_l_n (27) yields un:=HlnJnv∈ℝM×1u_n:=H_l_nJ_nv ^M× 1. If we define the auxiliary scalar loss function L~:=unTf L:=u_n^Tf such that its derivative equals dL~/df=und L/df=u_n, then once again by the chain rule; dL~dθ=JnTun=JnTHlnJnv∈ℝD, d Ldθ=J_n^Tu_n=J_n^TH_l_nJ_nv ^D, which equals the n-th term of (30) and can be computed using standard back propagation. Hence, each of the N contributions to GvGv can be evaluated using one forward- and one backward-pass through the network, without ever explicitly forming the Jacobian JnJ_n or the matrix G. In practice one uses mini-batches in place of a term-by-term computation. In addition, libraries such as PyTorch have dedicated subroutines for computing the required Jacobian-vector products. Access to GvGv enables scalable computation of the dominant eigenspace of G in high-dimensional settings. Briefly, we use the Lanczos algorithm, a well-known iterative method for approximating a few extremal eigenvalues and eigenvectors of a symmetric matrix G∈ℝD×DG ^D× D using only matrix–vector products. Starting from an initial unit vector v1v_1, it constructs an orthonormal basis v1,…,vl\v_1,…,v_l\ of the Krylov subspace spanv1,Gv1,…,Gl−1v1span\v_1,Gv_1,…,G^l-1v_1\ via the recursion vk+1=Gvk−bk−1vk−1−akvkbk,bk=∥Gvk−bk−1vk−1−akvk∥2,ak=vkTGvk. v_k+1= Gv_k-b_k-1v_k-1-a_kv_kb_k, b_k= Gv_k-b_k-1v_k-1-a_kv_k _2, a_k=v_k^TGv_k. (31) with b0=0b_0=0, which follows from a Gram-Schmidt procedure. This yields a low-dimensional symmetric tridiagonal matrix Tl∈ℝl×lT_l ^l× l with d<l≪Dd<l D constructed from aka_k, bkb_k coefficients (see e.g. Chen (2026)) for which Tl=VlTGVl T_l=V_l^TGV_l (32) holds, where Vl=[v1,…,vl]∈ℝD×lV_l=[v_1,…,v_l] ^D× l. The eigenvalues of TlT_l directly approximate the leading eigenvalues of G, and VlUl∈ℝD×lV_lU_l ^D× l approximate the corresponding eigenvectors, Ul∈ℝl×lU_l ^l× l being the eigenvectors of TlT_l. The first d columns of VlUlV_lU_l define our projection matrix P∈ℝD×dP ^D× d. Note that in order to bring the error down as in Figure 2, it is important to select a sufficiently high number of basis functions l>dl>d to span the Krylov subspace. We used l=d+10l=d+10. Crucially, note that the method requires only repeated evaluations of GvGv, making it well-suited for large-scale settings where G cannot be formed explicitly. 4.1.6 Empirical generalized Bayesian subspace inference for the prior variance As previously mentioned, the value of the prior variance plays an important role in (22), and must be chosen with care. However, defining a plausible prior on the connection weights θ is challenging due to their lack of physical meaning. Empirical Bayesian inference Carlin and Louis (2000) provides a solution in this context by allowing us to estimate σ02 _0^2 directly from the data y. That said, the severe overparameterization of neural networks again poses specific challenges. We found that applying empirical Bayes directly in the full, high-dimensional parameter space forces the optimization toward excessively high σ0 _0 values to compensate for the many D−dD-d flat, uninformative directions. To resolve this, we perform empirical Bayes in the projected parameter space defined by the active subspace w=PTθ∈ℝdw=P^Tθ ^d. Again, let P∈ℝD×dP ^D× d contain the dominant eigenvectors of G and define w0=PTθ0,θ(w)=θ0+P(w−w0). w_0=P^T _0, θ(w)= _0+P(w-w_0). (33) The loss restricted to the active subspace is then ℒN(w,y):=ℒN(θ(w),y) _N(w,y):=L_N(θ(w),y) Our analysis will proceed similar to Section 4.1.2. Let π~β,w π_β,w be the generalized unnormalized posterior in the w subspace, i.e.; π~β,w(w∣y,σ02)=exp[−βℒN(w,y)]p(w∣σ02), π_β,w (w y,σ^2_0 )= [- _N(w,y) ]p(w σ^2_0), (34) with the prior p(w∣σ02)=(0,σ02Id)p(w σ^2_0)=N (0,σ^2_0I_d ). We find σ02σ^2_0 by maximizing the generalized evidence: Zβ,w(y∣σ02)=∫π~β,w(w∣y,σ02)w. Z_β,w(y _0^2)= π_β,w (w y,σ^2_0 )dw. (35) The log density is given by logZβ,w(y∣σ02)=log∫exp[−lw(w∣σ02)]dw. Z_β,w(y _0^2)= [-l_w(w σ^2_0) ]\;dw. (36) with lw(w∣σ02)l_w(w σ^2_0) defined as the negative log unnormalized generalized posterior density: lw:=−logπ~β,w(w∣σ02)=βℒN(w,y)−logp(w∣σ02) l_w:=- π_β,w(w σ^2_0)= _N(w,y)- p(w σ^2_0) (37) We now insert a Taylor expansion for lwl_w around w0=PTθ0w_0=P^T _0 into (36), to give rise to a Laplace approximation; logZβ,w(y∣σ02)≈−lw(w0∣σ02)+log∫exp[−12(w−w0)T[∇w2lw(w0)](w−w0)]dw. Z_β,w(y _0^2)≈-l_w(w_0 σ^2_0)+ [- 12 (w-w_0 )^T[∇^2_wl_w(w_0)] (w-w_0 ) ]\;dw. (38) As in the full-space Laplace approximation, we assume that w0w_0 is sufficiently close to a stationary point of lwl_w, so that the first-order term in the Taylor expansion can be neglected. The second log term contains the integral of the unnormalized Laplace posterior in the active subspace, with ∇w2lw=PT∇θ2lβP=PT[β∇θ2ℒN−∇θ2logp(θ0∣σ02)]P=βPTHLP+σ0−2Id∇^2_wl_w=P^T\;∇^2_θl_β\,P=P^T [β∇^2_θL_N-∇^2_\ θ p( _0 σ^2_0) ]P=β P^TH_LP+σ^-2_0I_d, and therefore equals the reciprocal of the Gaussian normalizing constant. Defining Hw:=PTHLPH_w:=P^TH_LP, (38) then becomes; logZβ,w(y∣σ02) Z_β,w(y _0^2) ∝−lw(w0∣σ02)−12log(det[βHw+1σ02Id]) -l_w(w_0 σ^2_0)- 12 ( [β H_w+ 1 _0^2I_d ] ) ∝−βℒN(w0,y)+logp(w0∣σ02)−12log(det[βHw+1σ02Id]) - _N(w_0,y)+ p(w_0 σ^2_0)- 12 ( [β H_w+ 1 _0^2I_d ] ) ∝logp(w0∣σ02)−12log(det[βHw+1σ02Id]) p(w_0 σ^2_0)- 12 ( [β H_w+ 1 _0^2I_d ] ) ∝−12σ02∥w0∥22−dlogσ0−12∑i=1dlog(βλi+1σ02), - 12σ^2_0 w_0 ^2_2-d _0- 12 _i=1^d (β _i+ 1σ^2_0 ), (39) The loss ℒNL_N does not depend explicitly on σ02 _0^2, and is therefore considered constant in this approximation. In the last line we swap out HLH_L with an approximate Hessian discussed in Section 4.1.4 (i.e. Hw≈PTGPH_w≈ P^TGP), and the low-rank structure of HwH_w allows us to then compute the determinant using only the d dominant eigenvalues obtained from the matrix-free method of Section 4.1.5. Setting α=σ0−2α=σ^-2_0, we now find that 0 0 =dα(−α2∥w0∥22+d2logα−12∑i=1dlog(βλi+α)) = ddα (- α2 w_0 ^2_2+ d2 α- 12 _i=1^d (β _i+α ) ) =−12∥w0∥22+d2α−12∑i=1d1βλi+α. =- 12 w_0 ^2_2+ d2α- 12 _i=1^d 1β _i+α. (40) Solving for the projected weight norm; ∥w0∥22=∑i=1d1α−1βλi+α=∑i=1dσ02−σi2, w_0 ^2_2= _i=1^d 1α- 1β _i+α= _i=1^dσ^2_0-σ^2_i, (41) with σi2σ^2_i being the posterior variance (22). Therefore, the resulting σ02 _0^2 sets the prior scale such that the magnitude of the projected parameter vector is consistent with the amount of variance reduction in the active subspace. Note that performing this analysis in the full parameter space will lead to an expression involving ∥θ0∥22 _0 ^2_2 in place of ∥w0∥22 w_0 ^2_2, which will be very large for overparameterized neural networks. This then leads to an excessively large prior variance σ02σ^2_0. Although (41) defines an implicit and non-linear equation for σ02 _0^2, the corresponding function is smooth and strictly decreasing in α, and therefore admits a unique solution. This solution can be efficiently obtained using standard one-dimensional root-finding methods (we use the bisection method). 4.1.7 Calibrating d The dimension of the active subspace d is a hyperparameter that must be set. For regression problems, we will do so by evaluating the coverage of the output confidence intervals on a calibration data set (xn∗,yn∗)(x^*_n,y^*_n), n=1,⋯,Ncaln=1,·s,N_cal. Let Iγ(x)I_γ(x) be the 1−γ1-γ confidence interval for γ∈[0,1]γ∈[0,1]. Our value of d minimizes d∗=argmind1Ncal∑n=1Ncal[yn∗∈Iγ(xn∗)]≥1−γ. d^*= _d 1N_cal _n=1^N_cal 1 [y^*_n∈ I_γ(x^*_n) ]≥ 1-γ. (42) Here [A] 1[A] is an indicator function returning 1/0 if A is true/false. Note that for each d, we must perform the empirical Bayesian analysis of Section 4.1.6, such that d and σ0 _0 are estimated jointly. Calibration in the context of classification entails making sure the correct class in contained in the 1−γ1-γ set with (approximately) probability 1−γ1-γ. A similar calibration procedure can be applied to a classification network Adam. et al. (2024); Guo et al. (2017). However, if a classification network is already well calibrated, the additional stochasticity introduced via weight perturbations may diminish calibration, see the Discussion section. 4.1.8 Posterior prediction The posterior variance (22) corresponds to the Laplace approximation restricted to the dominant eigenspace of the curvature matrix. In the full parameter space, directions orthogonal to this subspace have variance σ02 _0^2, corresponding to the prior. While these directions are individually flat, they are vast in number (D−d≫dD-d d), and can collectively contribute significantly to predictive uncertainty, often resulting in overly wide confidence intervals when posterior samples are propagated through the non-linear network. Similar to e.g. Jantre et al. (2024); Faller and Martin (2026), we therefore omit these directions and obtain a low-rank approximation that concentrates uncertainty in the data-informed subspace. Sampling from this approximate posterior is then performed via θ=θ0+Pz,z∼(0,diag(σj2)),σj2=1βλj+σ0−2,j=1,⋯,d. θ= _0+Pz, z (0,diag(σ^2_j) ), σ^2_j= 1β _j+σ^-2_0, j=1,·s,d. (43) Using these samples during a forward pass we can approximate expectations such as f(x∣y)=∫f(x,θ)πβ(θ∣y)θ, f(x y)= f(x;θ) _β(θ y)dθ, (44) via Monte Carlo sampling, pushing the Laplace posterior through the network. As the Laplace posterior is an approximation, the non linear nature of the network may push the mean prediction away from the MAP prediction. This underfitting problem of sampled Laplace does not primarily mean that the Gaussian posterior family is insufficiently expressive. Rather, it refers to a degradation of the posterior predictive mean caused by propagating local Gaussian weight samples through the non-linear network. As mention in Section 4.1.4, linearized Laplace Immer et al. (2021) avoids this by making the a local linear approximation in the predictive model, thereby keeping the posterior predictive mean fixed at the MAP prediction: f(x,θ)≈flin(x,θ)=f(x,θ0)+J(x)(θ−θ0), f(x;θ)≈ f_lin(x;θ)=f(x; _0)+J(x) (θ- _0 ), (45) where J(x):=∇θf(x,θ0)∈ℝM×DJ(x):= _θf(x; _0) ^M× D. The combination of the linear output and the Gaussian distribution of θ also leads a closed-form predictive covariance Σf=JΣJT _f=J J^T, where Σ=PΛPT =P P^T, with Λ:=diag(σj2) :=diag(σ^2_j) is our low-rank Laplace posterior covariance. In the case of scalar regression the (pointwise) posterior predictive variance becomes σf2(x)=∑j=1dσj2(J(x)pj)2∈ℝ, σ^2_f(x)= _j=1^dσ^2_j (J(x)p_j )^2 , (46) with pj∈ℝDp_j ^D denoting the j-th column vector of P. Supplementary information For this article no Supplementary Materials are available. Acknowledgements We thank Dr Ketan Maheshwari at the US Department of Energy Oak Ridge National Laboratory and Dr Alessandro Tibo at AstraZeneca for helpful discussions and support with training REINVENT. We also thank Dr Laura Harbach at Brunel University London for initiating the REINVENT UQ study. Declarations Funding P.V.C. is grateful for funding from the UK Engineering and Physical Sciences Research Council under the following project “UK Consortium on Mesoscale Engineering Sciences (UKCOMES)” (Grant No. EP/R029598/1). Competing interests The authors declare no competing interests. Ethics approval and consent to participate Not applicable. Consent for publication Not applicable. Data and code availability The data used to train the neural networks, a pre-trained REINVENT model and Jupyter notebooks to reproduce the results of this article are available at Edeling (2026). Materials availability Not applicable. Author contribution W.N.E. developed the methodology, implemented the algorithms, performed the numerical experiments, analysed the results, and wrote the initial manuscript. P.V.C. contributed to the conceptual framing, interpretation of results, funding acquisition and manuscript revision. Both authors reviewed and approved the final manuscript. Declaration on the use of artificial intelligence tools The authors used OpenAI’s ChatGPT 5.5 as a tool for language editing. The tool was also used to discuss manuscript presentation, related-work positioning, and the interpretation of selected numerical results. All scientific claims, mathematical derivations, algorithms, references, numerical experiments, and conclusions were reviewed, verified, and finalized by the authors. The authors assume full responsibility for the content of the manuscript. 5 ChEMBL vocabulary Here we describe the ChEMBL vocabulary in Table 3 and display some example chemical compounds in Table 4. Table 3: Vocabulary used by the autoregressive SMILES generator. The vocabulary contains start and stop tokens, bond symbols, branch symbols, ring-closure indices, atom tokens, and charged or aromatic atom tokens. Index Token Description Index Token Description 0 $ Stop/end token 17 Br Bromine atom 1 Start token 18 C Aliphatic carbon atom 2 # Triple bond 19 Cl Chlorine atom 3 %10 Two-digit ring closure 20 F Fluorine atom 4 ( Branch opening 21 N Aliphatic nitrogen atom 5 ) Branch closing 22 O Oxygen atom 6 - Single bond 23 S Sulfur atom 7 1 Ring closure index 24 [N+] Positively charged nitrogen 8 2 Ring closure index 25 [N-] Negatively charged nitrogen 9 3 Ring closure index 26 [O-] Negatively charged oxygen 10 4 Ring closure index 27 [S+] Positively charged sulfur 11 5 Ring closure index 28 [n+] Positively charged aromatic nitrogen 12 6 Ring closure index 29 [nH] Protonated aromatic nitrogen 13 7 Ring closure index 30 c Aromatic carbon atom 14 8 Ring closure index 31 n Aromatic nitrogen atom 15 9 Ring closure index 32 o Aromatic oxygen atom 16 = Double bond 33 s Aromatic sulfur atom Table 4: 10 SMILES strings from the ChEMBL database. No. SMILES 1 C(C)C(NC(=O)C(NC(=O)C1OC2OC(C)(C)OC2C2OC(C)(C)OC12)C(C)C)C(O)=O 2 Fc1c(c1)C(=O)N1C(SCc2c(c2)N(=O)=O)=Nc2c2C1=O 3 C(C)CN1Cc2cnnn2-c2c(c2C1)-c1c1 4 CCN(C(=O)CSc1n1-c1c1)C1=C(N)N(Cc2c2)C(=O)NC1=O 5 C(C)N1CCN(C1)C(=O)N1c2c2Sc2c12 6 C1(C)Cc2c(OCc3c(c3)-c3n[nH]n3)c(Cl)c(Cl)c2C1=O 7 Cc1c(CSc2c(cn2)C(=O)Nc2c(F)c2)c1F 8 C(N1N=C(O)C2=Nc3c(Cl)c3C(=O)C2=C1O)c1ccccn1 9 C(N1C(=O)C2OC3(C)C=C(C)C=NC3C2C1=O)c1c1 10 Cc1nc(cs1)-c1c(NC(=O)c2c3nc4C(=O)NCCCn4c3c2)c1 References Adam. et al. (2024) Y. Adam., R. Maxime, W. Xi, and A. Laurence Bayesian low-rank adaptation for large language models. External Links: 2308.13111, Link Cited by: §3.1, §4.1.7. Aggarwal (2018) C.C. Aggarwal Neural networks and deep learning. Vol. 10, Springer. Cited by: §2.2. Avitabile et al. (2026) D. Avitabile, F. Cavallini, S. Dubinkina, and G.J. Lord Neural field equations with random data. SIAM/ASA Journal on Uncertainty Quantification 14 (2), p. 534–567. Cited by: §1. Bissiri et al. (2016) P.G. Bissiri, C.C. Holmes, and S.G. Walker A general framework for updating belief distributions. Journal of the Royal Statistical Society Series B: Statistical Methodology 78 (5), p. 1103–1130. Cited by: §1, §2.1, §4.1.1. Blundell et al. (2015) C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra Weight uncertainty in neural network. In International conference on machine learning, p. 1613–1622. Cited by: §1. Carlin and Louis (2000) B.P. Carlin and T.A. Louis Empirical Bayes: Past, present and future. Journal of the American Statistical Association 95 (452), p. 1286–1289. Cited by: §4.1.6. Chen (2026) T. Chen The Lanczos algorithm for matrix functions: a handbook for scientists. External Links: 2410.11090, Link Cited by: §4.1.5, §4.1.5. Cheung et al. (2011) S.H. Cheung, T.A. Oliver, E.E. Prudencio, S. Prudhomme, and R.D. Moser Bayesian uncertainty analysis with applications to turbulence modeling. Reliability Engineering & System Safety 96 (9), p. 1137–1149. Cited by: §1, §1. Constantine and Diaz (2017) P.G. Constantine and P. Diaz Global sensitivity metrics from active subspaces. Reliability Engineering & System Safety 162, p. 1–13. Cited by: §2.3. Constantine et al. (2014) P.G. Constantine, E. Dow, and Q. Wang Active subspace methods in theory and practice: applications to kriging surfaces. SIAM Journal on Scientific Computing 36 (4), p. A1500–A1524. Cited by: §1, §3.2, §4.1.4. Daxberger et al. (2021) E. Daxberger, A. Kristiadi, A. Immer, R. Eschenhagen, M. Bauer, and P. Hennig Laplace redux-effortless Bayesian deep learning. Advances in neural information processing systems 34, p. 20089–20103. Cited by: §3.2. Edeling et al. (2024) W.N. Edeling, M. Vassaux, Y. Yang, S. Wan, S. Guillas, and P.V. Coveney Global ranking of the sensitivity of interaction potential contributions within classical molecular dynamics force fields. npj Computational Materials 10 (1), p. 87. Cited by: §1, §3.1. Edeling (2026) W.N. Edeling GLAS: generalized laplace active subspaces. Note: https://github.com/wedeling/GLASGitHub repository Cited by: §2, Data and code availability. Faller and Martin (2026) J. Faller and J. Martin Low rank based subspace inference for the Laplace approximation of Bayesian neural networks. In The 29th International Conference on Artificial Intelligence and Statistics, External Links: Link Cited by: §1, §2.4, §3.2, §3.2, §4.1.8. Fu et al. (2019) H. Fu, C. Li, X. Liu, J. Gao, A. Celikyilmaz, and L. Carin Cyclical annealing schedule: A simple approach to mitigating KL vanishing. In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), p. 240–250. Cited by: §1. Gaulton et al. (2017) A. Gaulton, A. Hersey, M. Nowotka, A. Bento, J. Chambers, D. Mendez, P. Mutowo, F. Atkinson, L.J. Bellis, E. Cibrián-Uhalte, et al. The ChEMBL database in 2017. Nucleic acids research 45 (D1), p. D945–D954. Cited by: §2.4. Gelman (2013) A. Gelman Bayesian data analysis. Chapman and Hall/CRC. Cited by: §1. Gorlé and Iaccarino (2013) C. Gorlé and G. Iaccarino A framework for epistemic uncertainty quantification of turbulent scalar flux models for reynolds-averaged navier-stokes simulations. Physics of Fluids 25 (5). Cited by: §1. Grosse and Martens (2016) R. Grosse and J. Martens A kronecker-factored approximate fisher matrix for convolution layers. In International Conference on Machine Learning, p. 573–582. Cited by: §3.2, §3.2. Guo et al. (2017) C. Guo, G. Pleiss, Y. Sun, and K.Q. Weinberger On calibration of modern neural networks. In International conference on machine learning, p. 1321–1330. Cited by: §4.1.7. Gur-Ari et al. (2018) G. Gur-Ari, D.A. Roberts, and E. Dyer Gradient descent happens in a tiny subspace. arXiv preprint arXiv:1812.04754. Cited by: §4.1.4. Hu et al. (2021) E.J. Hu, Y. Shen, P. Wallis, Z. Allen-Zhu, Y. Li, S. Wang, L. Wang, and W. Chen Lora: low-rank adaptation of large language models. arXiv preprint arXiv:2106.09685. Cited by: §3.1. Immer et al. (2021) A. Immer, M. Korzepa, and M. Bauer Improving predictions of Bayesian neural nets via local linearization. In International conference on artificial intelligence and statistics, p. 703–711. Cited by: §1, §2.1, §2.2, §4.1.2, §4.1.4, §4.1.8. Izmailov et al. (2020) P. Izmailov, W.J. Maddox, P. Kirichenko, T. Garipov, D. Vetrov, and A.G. Wilson Subspace inference for Bayesian deep learning. In Uncertainty in Artificial Intelligence, p. 1169–1179. Cited by: §3.2. Jantre et al. (2024) S. Jantre, N.M. Urban, X. Qian, and B. Yoon Learning active subspaces for effective and scalable uncertainty quantification in deep neural networks. In ICASSP 2024-2024 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), p. 5330–5334. Cited by: §2.2, §3.2, §4.1.2, §4.1.8. Kingma and Ba (2014) D.P. Kingma and J. Ba Adam: a method for stochastic optimization. arXiv preprint arXiv:1412.6980. Cited by: §2.2. Kunstner et al. (2019) F. Kunstner, P. Hennig, and L. Balles Limitations of the empirical Fisher approximation for natural gradient descent. Advances in neural information processing systems 32. Cited by: §4.1.3, §4.1.4, §4.1.4, §4.1. Lanczos (1950) C. Lanczos An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. Journal of research of the National Bureau of Standards 45 (4), p. 255–282. Cited by: §4.1.5. Landrum et al. (2026) G. Landrum, P. Tosco, B. Kelley, D. Cosgrove, R. Vianello, E. Kawashima, G. Jones, A. Dalke, B. Cole, M. Swain, et al. Rdkit/rdkit: 2026_03_4. Zenodo. Cited by: §2.4. Lipinski et al. (1997) C.A. Lipinski, F. Lombardo, B.W. Dominy, and P.J. Feeney Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Advanced drug delivery reviews 23 (1-3), p. 3–25. Cited by: §3.1. Liu and Guillas (2017) X. Liu and S. Guillas Dimension reduction for Gaussian process emulation: An application to the influence of bathymetry on tsunami heights. SIAM/ASA Journal on Uncertainty Quantification 5 (1), p. 787–812. Cited by: §1. Loeffler et al. (2024) H. Loeffler, J. He, A. Tibo, J.P. Janet, A. Voronov, L.H. Mervin, and O. Engkvist Reinvent 4: Modern AI–driven generative molecule design. Journal of Cheminformatics 16 (1), p. 20. Cited by: §2.4. MacKay (1992) D.J.C. MacKay A practical Bayesian framework for backpropagation networks. Neural computation 4 (3), p. 448–472. Cited by: §1, §4.1.2. Murphy et al. (2004) J.M. Murphy, D.M.H. Sexton, D.N. Barnett, G.S. Jones, M.J. Webb, M. Collins, and D.A. Stainforth Quantification of modelling uncertainties in a large ensemble of climate change simulations. Nature 430 (7001), p. 768–772. Cited by: §1. Nafiz et al. (2024) A. Nafiz, S. Jantre, N.M. Urban, and B. Yoon Leveraging active subspaces to capture epistemic model uncertainty in deep generative models for molecular design. In 2024 IEEE 34th International Workshop on Machine Learning for Signal Processing (MLSP), p. 1–6. External Links: Link, Document Cited by: §3.2, §4.1.2. Olivecrona et al. (2017) M. Olivecrona, T. Blaschke, O. Engkvist, and H. Chen Molecular de-novo design through deep reinforcement learning. Journal of cheminformatics 9 (1), p. 48. Cited by: §2.4. Palmer (2000) T.N. Palmer Predicting uncertainty in forecasts of weather and climate. Reports on progress in Physics 63 (2), p. 71–116. Cited by: §1. Pearlmutter (1994) B.A. Pearlmutter Fast exact multiplication by the Hessian. Neural computation 6 (1), p. 147–160. Cited by: §4.1.5. Psaros et al. (2023) A.F. Psaros, X. Meng, Z. Zou, L. Guo, and G.E. Karniadakis Uncertainty quantification in scientific machine learning: Methods, metrics, and comparisons. Journal of Computational Physics 477, p. 111902. Cited by: §3.2. Ritter et al. (2018) H. Ritter, A. Botev, and D. Barber A scalable Laplace approximation for neural networks. In International Conference on Learning Representations, External Links: Link Cited by: §3.2, §4.1.2. Selva et al. (2016) J. Selva, R. Tonini, I. Molinari, M.M. Tiberti, F. Romano, A. Grezio, D. Melini, A. Piatanesi, R. Basili, and S. Lorito Quantification of source uncertainties in seismic probabilistic tsunami hazard analysis (sptha). Geophysical Journal International 205 (3), p. 1780–1803. Cited by: §1. Wang et al. (2026) K. Wang, F. Cuzzolin, K. Shariatmadar, D. Moens, and H. Hallez A review of uncertainty representation and quantification in neural networks. IEEE Transactions on Pattern Analysis and Machine Intelligence 48 (3), p. 2476–2495. Cited by: §3.2. Yeh (1998) I. Yeh Concrete Compressive Strength. Note: UCI Machine Learning RepositoryDOI: https://doi.org/10.24432/C5PK67 Cited by: §2.3.