Uncertainty quantification and propagation in atomistic machine learning

2.1 Distribution functions

A probability density function (PDF), denoted as p(θ), is used to describe the probability distribution of a continuous random variable Θ. The PDF represents the relative likelihood of the random variable taking on a specific value. A cumulative distribution function (CDF), denoted as F(θ), describes the probability that a random variable Θ takes a value less than or equal to θ. The CDF can be obtained from the PDF via: F ( θ ) = P ( Θ ≤ θ ) = ∫ − ∞ θ p ( t ) d t . A quantile function (QF) maps probabilities to values of the random variable. It is defined as the inverse of the CDF, Q(p) = F⁻¹(p), that is, for a given probability p, the quantile function returns the value θ such that its probability is less than or equal to an input probability value, i.e., P(Θ ≤ θ) = p. As a concrete example, Figure 1 presents the PDF, CDF, and QF of a one-dimensional standard Gaussian distribution N ( 0,1 ) with a mean of 0 and standard deviation of 1.

Figure 1:

One-dimensional Gaussian distribution with a mean of μ = 0 and standard deviation σ = 1. The probability density function (PDF), cumulative probability function (CDF), and quantile function (QF). The QF is the inverse of the CDF, as can be seen by switching the horizontal and vertical axes.

2.2 Bayesian inference

In the Bayesian view, probability provides a quantification of uncertainty (Gelman et al. 2013). Given an observed dataset D, we are interested in obtaining the conditional probability p(y*|x*, D) of the output y* for a new input x*. From it, a statistical measure of the uncertainty in y* such as the variance can then be obtained. Given a parametric model, y = f(x; θ), this can be achieved in two steps.

First, obtaining the posterior distribution over model parameters θ using Bayes’ theorem (Gelman et al. 2013):

The prior distribution p(θ) represents our prior information as to which parameters θ are likely to have generated the outputs before observing any data, based on previous knowledge, experience, or physical limitations. The effects of the observed data D come from the likelihood function, p(D|θ), which quantifies the plausibility of D for different realizations of θ. The denominator,

is called the marginal likelihood, also known as the evidence in the context of Bayesian statistics, which ensures that the posterior is a proper probability and thus integrates into one. It represents the likelihood of the observed data D, considering all possible values of the parameter θ weighted by their prior distribution. Marginalization means evaluating this equation to obtain the marginal likelihood. Bayes’ theorem converts the prior probability over model parameters into the posterior probability by incorporating the evidence provided by the observed data.

Second, with the posterior over θ, we can obtain the predictive distribution for a new data point (x*, y*) as

where p(y*|x*, θ) is the likelihood given the new data point. It is related to the likelihood function for a dataset p(D|θ) in Eq. (1), and further discussion on this will be provided in Section 2.3. From the predictive distribution, we can readily obtain, e.g., the mean as the final prediction and the variance as a point estimate of the uncertainty.

Although theoretically sound, a major practical limiting factor of the full Bayesian approach lies in the complexity of evaluating the posterior. In particular, the marginal likelihood in Eq. (2) can only be analytically evaluated for simple models like linear regression. Numerical techniques such as sampling methods and variational inference have to be undertaken to evaluate the predictive distribution for more complicated models (Bishop 2006).

One can sample the posterior distribution using Monte Carlo (MC) methods, such as Markov chain Monte Carlo (MCMC), and then obtain the predictive mean and uncertainty (Neal 1993, 2003). Sampling methods are accurate, flexible, and can be applied to a wide range of models. However, they are still computationally intensive because a large number of samples might be needed for convergence; thus, they are mainly used for small-scale problems (Bishop 2006).

Alternatively, the variational inference approach tackles the challenge by employing another distribution q(θ) to approximate the true posterior p(θ|D), and then using q(θ) to evaluate the predictive distribution in Eq. (3). The approximate distribution q(θ) is typically much simpler, and it is optimized to resemble the true posterior, e.g., by minimizing the Kullback–Leibler divergence (Kullback and Leibler 1951; MacKay 1992b) between q(θ) and p(θ|D). Although not exact, variational inference offers a computationally efficient approach to evaluate the predictive distribution. The MC dropout method to obtain uncertainty in NN models proposed by Gal (2016) adopts this approach, and it will be further discussed in Section 3.1.1.

2.3 Maximum likelihood

Maximum likelihood estimation (MLE) is a frequentist approach to estimate the optimal parameters θ of a model y = f(x; θ). It is equivalent to the least-squares parameter optimization technique widely used in science and engineering. MLE provides a point estimate of the parameters and thus predictive uncertainty cannot be directly quantified from a model trained using MLE alone. However, MLE serves as a fundamental concept in parameter estimation and forms the basis of many other UQ methods.

In MLE, we focus on the likelihood function p(D|θ), and do not care about the prior and posterior (see Eq. (1)). We assume that the observed output y is given by the model prediction f(x; θ) with an additive error ϵ, i.e., y = f(x; θ) + ϵ. A Gaussian distribution with zero mean is a reasonable choice for the error, p ( ϵ ) = N ( ϵ | 0 , σ 2 ) , where σ² is the variance. This is equivalent to the Gaussian distribution in which y is regarded as the random variable with the model prediction y ˆ = f ( x ; θ ) as its mean and σ² as the variance (Bishop 2006):

In other words, this is the likelihood of θ for a single data point (x, y). Now consider the observed dataset D where each data point is drawn independently from the distribution in Eq. (4). We can obtain the likelihood function for the dataset as the product of the likelihood for each data point:

The optimal model parameters can thus be obtained by maximizing Eq. (5) with respect to (w.r.t.) θ, which is equivalent to minimizing the negative log-likelihood (NLL):

because the logarithm function is monotonically increasing. This transformation converts a product of probabilities into a sum of log probabilities, which is often more convenient to optimize.

We note that in MLE, the variance σ² is modeled as a single constant, albeit unknown. Consequently, minimizing the NLL is equivalent to the familiar least-squares minimization w.r.t. θ using the loss:

Recall that the model parameters θ are implicitly indicated in the model prediction y ˆ = f ( x ; θ ) . With the optimal point estimate of the parameters, θ_opt, we can then get model prediction as y* = f(x*; θ_opt) for any new input x*. Again, we note that no information on uncertainty can be directly obtained from this approach only.

2.4 Evidence approximation

While the full Bayesian approach can produce predictive uncertainty, it can be computationally demanding to obtain. On the other hand, it is straightforward to get a point estimate of the optimal model parameters using MLE, but the uncertainty cannot be quantified. The evidence approximation approach (Gull 1989; MacKay 1992a), also known as the empirical Bayes (Bernardo and Smith 2009), generalized maximum likelihood (Wahba 1985), or type 2 maximum likelihood (Berger 1985), lies between the two extremes.

Evidence approximation is a method that looks at how well a model fits the data overall. It focuses on the marginal likelihood in Eq. (2), which provides the model evidence of observing the data marginalized over the parameters, meaning that it calculates the probability of seeing the observed data under all possible parameter values of the model. In evidence approximation, the prior p(θ) in Eq. (2) is further parameterized using a set of hyperparameters ξ, becoming p(θ|ξ). Consequently, the marginal likelihood is (Bishop 2006):

The introduction of the hyperparameters ξ increases the model’s capacity, meaning that the model has increased flexibility to capture the underlying structure present in the data. In practice, the prior distribution p(θ|ξ) is often chosen to be conjugate to the likelihood p(D|θ, ξ). In other words, p(θ|ξ) is specifically selected to match the form of the likelihood p(D|θ, ξ) such that the integration in Eq. (8) has a closed-form solution, thereby simplifying the process to obtain p(D|ξ). For example, if the likelihood is the Binomial distribution, a conjugate prior for it is the Beta distribution, and then Eq. (8) can be analytically evaluated to obtain the model evidence, which is a Beta distribution as well (DeGroot and Schervish 2012).

In a full Bayesian setting, after obtaining p(D|ξ), one would also marginalize over the hyperparameters ξ to obtain p(D) and then perform Bayesian inference. However, this marginalization typically does not have an analytical solution. Instead, in evidence approximation, we get a point estimate of ξ by maximizing the model evidence p(D|ξ) w.r.t. ξ. Then, the model evidence, prior, and likelihood are all evaluated using the optimal hyperparameters ξ_opt to perform Bayesian inference using Eqs. (1) and (3).

Alternatively, one can obtain the predictive uncertainty directly from p(θ|ξ_opt), provided that the prior is chosen to be of a specific form, e.g., as a high-order distribution on top of the likelihood. This is the evidential regression approach recently proposed by Amini et al. (2020). In Section 3.1.3, we will further discuss this approach and explain how to get the uncertainty.

3.1 Probabilistic approach

3.1.1 Monte Carlo dropout

The dropout technique was originally proposed by Srivastava et al. (2014) as a regularization technique to alleviate overfitting in deep NN models. It was adapted by Gal and Ghahramani (2016) to approximate the Bayesian approach for UQ mentioned in Section 2.2. Their method, known as MC dropout, can be theoretically viewed as sampling from a Bernoulli prior distribution over the weights of the NN, and then taking advantage of the variational inference technique to approximate the posterior distribution. Practically, dropout is used at both training and inference time, allowing the model to estimate uncertainty by considering multiple predictions from different subsets of the network.

The model can be trained by minimizing a loss between its predictions and the corresponding reference values to obtain the optimal NN parameters. At each training step, dropout randomly sets the outputs of a fraction of the nodes in an NN to zero (e.g., the dashed nodes in Figure 3a), effectively creating an ensemble of thinned sub-networks. Once trained, we use the model in a similar way. Multiple forward passes are performed through the NN, each with different nodes being dropped, to obtain multiple predictions. The predictions are then averaged to obtain the final prediction and their variance is computed as the predictive uncertainty. In this sense, MC dropout can also be thought of as a frequentist ensemble approach to be discussed in Section 3.2.

Figure 3:

Schematic illustration of the probabilistic UQ approaches: (a) MC dropout; (b) MVE; and (c) evidential regression.

It is important to note that MC dropout is an approximation and may not capture the full posterior distribution over the NN’s weights. Nevertheless, it has gained popularity as a practical technique in atomistic ML. For example, Wen and Tadmor (2020) have developed a dropout NN model for carbon allotropes and shown that the obtained uncertainty can reliably distinguish diamonds from graphene and graphite.

3.2 Ensemble approach

The ensemble approach, characterized by its simplicity and diverse construction methods, combines multiple models to create a more robust predictive model, surpassing individual model limitations (Zhou 2012). The frequentist ensemble approach is easy to implement and can be applied to a large number of regression algorithms. It might be computationally expensive when compared to other UQ methods, but they can be naively paralleled. A couple of methods exist to construct an ensemble, such as parameter initialization, bootstrapping, and subsampling (Figure 4).

Figure 4:

Schematic illustration of the ensemble UQ approaches: (a) parameter initialization, (b) bootstrapping, and (c) subsampling.

The first type of method involves fitting models with different parameter initializations to the same dataset. Lakshminarayanan et al. (2017) proposed the use of ensembles for estimating the uncertainty of NNs. NNs of the same structure are created, but their parameters are initialized to be different. Each member of the ensemble is trained on the entire training set, meaning that all members have access to the same data.

The second type of method focuses on fitting the same model to different datasets. Bootstrapping is a widely used method to generate multiple derived datasets from a given dataset. Each subset, called a bootstrap sample, is created by randomly sampling the same number of data points from the original dataset with replacement (Hastie et al. 2009). This means that each data point has an equal probability of being selected, and the same data point can be included multiple times in the bootstrap sample. This process is repeated M times, resulting in M bootstrap samples. Then, M models are trained separately, each using one of the bootstrap samples.

Subsampling (Politis and Romano 1994; Politis et al. 1999) is an alternative to bootstrapping to create multiple derived datasets from a given dataset. It is similar to bootstrapping but with a key difference: subsampling is performed without replacement and thus each data point can only appear at most once in each subset. As a result, each sample consists of fewer data points than the original dataset.

The final ensemble prediction and uncertainty are obtained by combining the outputs y ˆ 1 , y ˆ 2 , … , y ˆ M of all members. The mean,

gives the final prediction, and the variance,

gives the predictive uncertainty.

3.3 Feature space distance approach

Another category of UQ approach is based on some distance measure in the model’s feature/latent space. It assumes that data points resembling each other are positioned closer to one another in the feature space. Therefore, for a given test data point, if it is close to the training data in the feature space, the predictive uncertainty is low; otherwise, the predictive uncertainty is high.

In this approach, the uncertainty is typically obtained in a two-step process. First, given a primary model like an NN, a point estimate of the optimal model parameters is obtained using a training technique such as MLE. Second, construct another model to measure the distance between a test data point and the training data in the primary model’s feature space. When dealing with NNs, the feature space can be chosen as the last but one layer or other internal layers that are considered suitable.

3.3.1 LTAU

Loss trajectory analysis for uncertainty (LTAU) measures the distance in the Euclidean space. It begins by training an NN model to predict an atomic property y, chosen to be the forces on atoms in Vita et al. (2024). During training, besides optimizing the model parameters, the error ϵ i = ‖ y i − y ˆ i ‖ 2 between the model prediction y ˆ i and its corresponding reference y _i for each atom i is recorded as

(18) T i = ϵ i 1 , ϵ i 2 , … , ϵ i E ,

where the super index denotes the training epoch at which the error is logged, and E is the total number of epochs to train the model.

After training, we get a set of errors T _i along the loss trajectory for each data point i and then convert the errors to the model’s confidence score for that data point via

(19) p i = P ϵ i e ∈ T i ≤ a t o l ,

which means the ratio of data points in T _i whose value is smaller than or equal to a tolerance atol. A reasonable choice for atol would be the mean absolute error (MAE) between the predictions from the trained primary model and their references. The value of p _i is in the range of [0, 1]. The uncertainty of each data point is calculated as

(20) δ i = 1 − p i ,

which can be interpreted as the probability that the model’s prediction on data point i will have an error larger than the MAE.

The uncertainty for a new data point, j, is obtained by averaging the uncertainties of its nearest K neighbors in the training data (Figure 5a):

(21) δ j = 1 K ∑ i ∈ N j δ i ,

where N _j denotes the set of neighbors. The nearest neighbors can be determined by performing a similarity search based on Euclidean distance in the feature space.

Figure 5:

Schematic illustration of the UQ approaches based on feature space distance. (a) LTAU assigns the average uncertainty of neighboring training data (circles) to a test point (square). (b) GMM models the density of the training data in the feature space, and a test point has low uncertainty if it is located in a dense region.

LTAU has been successfully applied to tune the training–validation gap in NN potentials for carbon materials and predict the errors in relaxation trajectories of catalysts (Vita et al. 2024).

4.1 Calibration

Calibration measures the statistical consistency between the predictions and observations, a property that depends on both the predictions and observations (Gneiting et al. 2007). Uncertainty calibration has been extensively studied in the context of classification problems. Perfect calibration in this setting means that the confidence assigned to a class equals the probability of the prediction belonging to that class (Guo et al. 2017; Scalia et al. 2020). For instance, if we have 10 predictions, each with a confidence of 0.8, we expect 8 of them to be correctly classified.

Uncertainty calibration for regression is less intuitive because the model predicts continuous values, rather than discrete labels as in classification. Nevertheless, following the groundbreaking work by Gneiting et al. (2007), methods that extend the uncertainty calibration approach for classification have been proposed for regression problems (Kuleshov et al. 2018; Levi et al. 2022) and adopted in atomistic ML. Here, we discuss two such methods: interval based and error based regression uncertainty calibration. As a technical note, we will use δ to denote uncertainty in general. However, for some models, uncertainty is represented by the variance σ² (see Table 2). Therefore, these two notations will be used interchangeably when appropriate.

4.1.1 Interval based approach

Loosely speaking, in a regression setting, calibration means that a model prediction should fall in a given confidence interval γ% approximately γ% of the time (Kuleshov et al. 2018). For example, the model in the top panel of Figure 7 is not calibrated because only 20 % (2 out of 10) of the time the predictions are within the 90 % confidence interval, while the one in the bottom is calibrated. Formally, according to Kuleshov et al. (2018), for a given calibration dataset D cal = { ( x i , y i ) } i = 1 N , a regression model is calibrated if

(25) ∑ i = 1 N I y i ≤ F i − 1 ( p ) N → p   for all p ∈ [ 0,1 ] ,

as N → ∞. Here, F _i  = P(Y < y _i ) denotes the CDF of the random variable Y, and F i − 1 is the corresponding quantile function (see Section 2.1). I [ c ] is the indicator function, which evaluates to 1 if the condition c is true and 0 otherwise. In other words, Eq. (25) means that, for a calibrated model, the empirically observed CDF from the data and the expected CDF by the model should match as the dataset size goes to infinity.

Figure 7:

Illustration of a non-calibrated regression model (top) and a calibrated one (bottom). The plot is inspired by Kuleshov et al. (2018) but generated using arbitrary data.

In practice, Eq. (25) is evaluated for a selected number of p values, and the calibration curve is used to check the calibration level, which plots the observed CDF from the data versus the expected CDF by the model (Figure 8a). With this, the calibration curve can be obtained as follows (Kuleshov et al. 2018):

1. Discretize the expected CDF to a set of M values 0 < p₁ < p₂… < p _M  < 1. The data is assumed to be generated from a Gaussian, y ∼ N ( y ˆ , σ 2 ) , where y ˆ and σ² are the predictions and the associated uncertainty, respectively (Kuleshov et al. 2018; Tran et al. 2020). For example, if an ensemble method is used, then y ˆ and σ² are the ensemble mean and variance, respectively. For a Gaussian with mean y ˆ and variance σ², the expected CDF can be readily obtained (see Section 2.1). We note that assuming a Gaussian distribution may not accurately reflect the true nature of the data in all cases.

2. For each expected CDF p _j , compute the corresponding expected model output y _j . As mentioned above, the model output is assumed to follow a Gaussian; therefore, y _j can be readily computed using the quantile function, y _j  = Q(p _j ) = F⁻¹(p _j ), discussed in Section 2.1.

3. For each expected CDF p _j , compute the corresponding observed CDF p ˜ j . With y _j obtained in the previous step, p ˜ j is obtained as the empirical frequency p ˜ j = | { y i | y i ≤ y j , i = 1,2 , … N } | / N (left of Eq. (25)), where |⋅| denotes the size of a set, and N is the total number of data points in the calibration dataset D_cal. In other words, p ˜ j is computed as the fraction of data points whose prediction y _i is smaller than or equal to y _j .

4. Create the calibration curve by plotting ( p j , p ˜ j ) pairs for j = 1, 2, …M.

Figure 8:

Calibration curves and probability density for the interval approach. (a) Calibration curves for perfectly calibrated (diagonal grey), over-confident (green) and under-confident (purple) models. (b) The Gaussian probability densities correspond to the calibration curves in (a).

The calibration curve provides rich information. First, for a perfectly calibrated model as defined in Eq. (25), the calibration curve should be a diagonal line, meaning that the observed CDF from the data and the expected CDF by the model match with each other. Therefore, a model’s calibration could be qualified by the closeness of its calibration curve to the diagonal line. In addition, the shape of the calibration curve could yield other insights into the predictive uncertainty of a model. A calibration curve that is above the diagonal line at low expected CDF but below at high expected CDF suggests that the model is over-confident. To understand this, let’s focus on the low expected CDF region, e.g., at 0.2 in Figure 8a. Here, the expected CDF is smaller than the observed CDF, meaning that the variance in the Gaussian distribution used to construct the model is smaller than the variance in the observed data (Figure 8b). With a smaller expected variance (uncertainty), the model is over-confident. On the other hand, an under-confident model has a calibration curve that is below the diagonal line at a small expected CDF but above a large expected CDF.

In Kuleshov et al. (2018) and Levi et al. (2022), the observed CDFs are interpreted as observed confidence intervals. Because of the use of CDF, the interval here means − ∞ , q j . This is different from the commonly used notion of confidence interval, which is typically specified as an interval around the mean. For example, the 68 % confidence interval of a Gaussian distribution is μ ± σ. Thus, to avoid confusion, we do not use “confidence interval” but directly use CDF as is also done in Tran et al. (2020). Nevertheless, it is possible to interpret the calibration curve as the commonly used confidence interval. Instead of CDF, one can employ the probability density function (PDF), considering symmetric intervals μ ± γ _j of varying confidence level 0 < γ₁ < γ₂… < γ _M  < 1 around the mean and examining the empirical frequency of the observed data belonging to each interval (Scalia et al. 2020).

4.1.2 Error based approach

The error based approach directly compares the predicted uncertainty σ² and the expected square error between the model prediction y ˆ and the observed data y, stating that, for a calibrated model, the predicted uncertainty and the expected error should match (Levi et al. 2022). Formally, a regression model is calibrated if

(26) E x , y ( y ˆ − y ) 2 | σ 2 = u → u

for any chosen positive u, where the expectation E is taken over the joint distribution of x and y. From the definition, no average over points with different values of u is needed; thus, in principle, for each data point, one can correctly predict the expected error. In practice, however, binning is performed to empirically evaluate Eq. (26).

Similar to the calibration curve in the interval based approach, here, a reliability diagram can be created to diagnose the calibration level of a model, as follows (Levi et al. 2022):

1. Sort the data points according to their predicted uncertainty σ², and then divide them into M bins, B₁, B₂, …B _M . For simplicity, the bin boundaries can be equally located from the minimum to the maximum of the uncertainty (Figure 9a).

2. For each bin B _j , calculate the root-mean variance (RMV), RMV ( j ) = 1 | B j | ∑ i ∈ B j σ i 2 , and the root-mean-square error (RMSE), RMSE ( j ) = 1 | B j | ∑ i ∈ B j ( y ˆ i − y i ) 2 , where |B _j | is the number of data points in bin j.

3. Plot the RMSE(j) against the RMV(j) for each bin j. This plot is the reliability diagram (Figure 9b).

Figure 9:

Reliability diagram for error based calibration approach. (a) Binning an example dataset of 400 data points into 10 equally separated intervals. (b) Calibration curve, where each blue dot represents the RMSE and RMV for each bin. A perfectly calibrated model should follow the diagonal line.

According to Eq. (26), if a model is perfectly calibrated, the RMV and RMSE should be equal for each bin; therefore, the calibration curve should be a straight diagonal line in Figure 9b. A larger deviation from the diagonal line suggests a more poorly calibrated model. It is important to note that, unlike the interval based approach, here, the reliability diagram is not constrained to be within [0, 1], but instead ranges between 0 and the maximum RMV or the maximum RMSE value. Consequently, directly comparing the reliability diagrams of different models is not appropriate unless some form of normalization is performed. Furthermore, the choice of the number of bins can have a significant impact on the results. For instance, in Figure 9a, we chose 10 bins, and the last bin B₁₀ consists of only three data points, which is insufficient to obtain reliable statistics for RMV and RMSE. One could consider using a smaller number of bins; however, if the number of bins is too small, the details might be averaged out.

4.2 Metrics

Calibration plots offer a qualitative way to assess UQ methods. For quantitative comparison, scoring metrics that can assign numerical scores to each UQ method become necessary. Various metrics have been proposed, and we focus on the important ones that evaluate a UQ method from four different perspectives: calibration, precision, accuracy, and efficiency.

4.3 Benchmark studies

Using metrics as those discussed in Section 4.2, several benchmark studies have attempted to evaluate the performance of various UQ methods on diverse atomistic ML datasets (Hirschfeld et al. 2020; Hu et al. 2022; Scalia et al. 2020; Tan et al. 2023; Tran et al. 2020; Varivoda et al. 2023). A major goal is to identify UQ methods that can perform well across metrics and datasets, thus providing practical guidance for selecting appropriate ones for chemical and materials applications. Several general observations can be made from these benchmark studies.

First, the UQ methods perform differently on different metrics; a method that works well on one metric can fall short on another. For example, Tran et al. (2020) have trained various ML models to predict the adsorption energies of small molecules on metal surfaces calculated from DFT. Although all trained ML models reported in Tran et al. (2020) have similar accuracy (MAE of ∼0.20 eV) and miscalibration area (∼0.13, Figure 10), their performance on precision varies a lot. For example, MC dropout is much sharper than NN ensemble (SHA: 0.09 versus 0.14), but has a lower dispersion (C _v : 0.82 versus 1.06) (see Section 3, Tran et al. (2020), for more details). Even for UQ methods that have, for example, a similar miscalibration area as illustrated in Figure 10, the shape of the calibration curves can be drastically different, indicating different modes of miscalibration. Tan et al. (2023) observed similar behavior using the rMD17 dataset (Christensen and Von Lilienfeld 2020) of energies of small molecules.

Figure 10:

Comparison of calibration curves displaying miscalibration area of various UQ methods. NN ensemble, here, refers to NNs trained with random parameter initializations. GP: Gaussian process, a probabilistic modeling approach that inherently models predictive uncertainty (Rasmussen 2003). Images adapted from Tran et al. (2020) under a CC BY 4.0 DEED license.

Second, performance varies by dataset; for a given metric, different UQ methods can have varying error levels across datasets. For example, Tan et al. (2023) observed that, for the rMD17 dataset of energies of small molecules, NN ensemble exhibits smaller miscalibration error than the evidential regression method (Figure 11b). Varivoda et al. (2023) have found similar behaviors for the dataset of formation energies of crystals (Jain et al. 2013) and the dataset of surface adsorption energies of metal alloys (Mamun et al. 2019). However, they found that, for the dataset of band gaps of MOFs (Rosen et al. 2021), the miscalibration errors are the same for the ensemble and evidential regression methods.

Figure 11:

Comparison of various UQ methods on the rMD17 dataset. (a) Predicted uncertainties as a function of squared errors of atomic forces and (b) computational cost, MAE of atomic forces, Spearman correlation, and miscalibration area of various UQ methods. Ensemble, here, refers to NNs trained with random parameter initialization. Images adapted from Tan et al. (2023) under a CC BY 4.0 license.

Third, ensemble methods appear to be a reliable UQ approach in general. For example, Tan et al. (2023) compared the ensemble method versus the MVE, evidential regression, and GMM deterministic methods for both in-domain (rMD17 and ammonia) and out-of-domain (silica glass) tasks. Using metrics such as MAE, Spearman correlation, and miscalibration area (Figure 11), they have found that single-deterministic methods struggle to consistently perform better across each in-domain and out-of-domain task and that the ensemble method still remains the most reliable choice. For ensemble methods, different ways to generate the ensemble can result in different performances. For example, Scalia et al. (2020) have compared three ensemble methods – NN ensemble with different parameter initialization, bootstrapping, and MC dropout – using the miscalibration area, ENCE, sharpness, and dispersion metrics (Recall from 3.1.1 that MC dropout can be regarded as an ensemble method in practice). Evaluated on various MoleculeNet benchmarking datasets (Wu et al. 2018), their finds indicate that NN ensemble with different parameter initialization and bootstrapping consistently outperform MC dropout (Figure 12).

Figure 12:

Comparison of error based calibration curves of three UQ methods on the QM9 dataset. NN ensemble, here, refers to NNs trained with random parameter initializations. Image adapted from Scalia et al. (2020). Copyright 2020, American Chemical Society.

Despite the robustness of the ensemble approach, there are contradictory studies showing that ensemble methods may not always be the most reliable choice. For example, Hirschfeld et al. (2020) have demonstrated that the stacking methods that sequentially combine multiple weak models to produce the uncertainty estimate can be more consistent than the ensemble approach. Heid et al. (2024) have shown that the uncertainty of a single prediction obtained from an ensemble approach cannot be directly correlated with the absolute error per atom. This is because the absolute error is distributed along a normal distribution, with its width determined by the uncertainty arising from model variance. To address this, they developed an approach that uses locally aggregated uncertainties to identify high-error local substructures, enabling the resolution of absolute errors on an atomic scale.

Fourth, off-the-shelf metrics may not be directly applicable. The metrics, particularly those related to calibration discussed in Section 4.2, have been primarily developed within the ML community using non-chemical datasets. Their suitability for chemical and materials problems is not guaranteed. For instance, the error based calibration approach conventionally performs binning directly on the predicted uncertainty, represented by the variance σ² in Figure 9a. Chemical and materials data can have numerous small variance values. Therefore, it might be more appropriate to conduct binning after transforming the variances using a logarithmic function to avoid cluttering the bins at small variance values, as demonstrated in Figure 11a.

Benchmark studies so far have been valuable in highlighting how different UQ methods can have varied performances depending on the choice of metrics and dataset. However, these studies consistently indicate that the question of which UQ method to select in practice remains unresolved. This appears to be discouraging. Nevertheless, there are guiding principles based on ease of use, efficiency, and uncertainty propagation. Ensemble methods provide a strong baseline and are straightforward to implement, making them a go-to choice for many applications. However, single-pass models (e.g., MVE, evidential, and GMM) are generally more efficient than ensemble methods (Figure 11b), so they can be good choices for resource-bounded applications. Another consideration is UP; the chosen propagation method may make certain UQ methods more suitable. This will be further discussed in Section 5.

4.4 Recalibration

So far, we have discussed ways for evaluating UQ methods in terms of their calibration, precision, accuracy, and efficiency, and we have provided some examples. But if a model’s performance is unsatisfactory, are there any approaches to improve it? The answer is yes, and, here, we will concentrate on calibration.

Informally, the calibration problem can be described as follows: given a trained model U that can generate predictive uncertainty δ = U(x), we train another recalibration model R such that the output of the composed functions δ ˆ = R ( U ( x ) ) is calibrated. The model R should be trained on a separate recalibration dataset D_cal distinct from the datasets used for model parameter optimization or hyperparameter tuning. Below, we discuss two recalibration methods, both of which are applicable to the interval based calibration approach introduced in Section 4.1.1.

Variance scaling. For the interval based calibration, the model prediction, y, is set to take the form of a Gaussian  y ∼ N ( μ , σ 2 ) . The model can be under-confident or over-confident depending on the scale of the variance σ². To recalibrate it, that is, making the calibration curve move toward the diagonal line in Figure 8a, we can train a linear model R ( σ 2 ) : σ ˆ 2 = a σ 2 + b to scale the variance. The parameters a and b can be determined by, e.g., minimizing a calibration NLL loss (Tan et al. 2023), 1 2 ∑ i = 1 N log [ 2 π ( a σ 2 + b ) ] + ( y i − y ˆ i ) 2 / ( a σ 2 + b ) , using a recalibration dataset D_cal consisting of N data points. Once the optimal a and b are obtained, the new variance σ ˆ 2 for the Gaussian is known for every data point, which can then be used to regenerate the calibration plots.

Isotonic regression. The variance scaling approach still assumes that the model prediction follows a Gaussian distribution. The true observed data distribution, however, may not be Gaussian. Kuleshov et al. (2018) proposed a recalibration approach based on isotonic regression, which is effective even for non-Gaussian cases. Given a recalibration dataset D_cal, this approach begins by transforming it into a processed recalibration dataset D ˆ cal = { ( q j , q ˜ j ) } j = 1 M , where q _j and q ˜ j are obtained using the same procedures described in Section 4.1.1. Using D ˆ cal , an isotonic regression model q ˜ j = R ( q j ) is then trained. Isotonic regression is a technique for fitting a free-form line to map a sequence of inputs (q _j in this case) to a sequence of observations ( q ˜ j in this case) such that the fitted line is non-decreasing everywhere and lies as close to the observations as possible (Fielding et al. 1974). Isotonic regression is chosen because it accounts for the fact that the true calibration curve is monotonically increasing. Once the isotonic regression model is trained, new uncertainty δ ˆ can be generated, and new calibration plot can be produced.

5.1 Bayesian propagation

Given the posterior over model parameters p(θ|D) (see Section 2.2), the distribution of a QoI z can be written as

where p(z|x, θ) is the likelihood of θ to observe z given the composed model z = h(x; θ). The integration, generally, cannot be analytically evaluated for most QoI in MD and microkinetic simulations. To address this, again, sampling techniques can be employed. Eq. (33) can be approximated by (Angelikopoulos et al. 2012)

Assuming the likelihood p(z|x, θ) is a Gaussian with h(x; θ) as its mean, then the predictive mean and variance of Eq. (34) can be respectively expressed as

and

The sampling-based Bayesian UP means selecting multiple sets of model parameters θ, computing the QoI z _i using each parameter set, and then calculating their mean as the final prediction and their variance as the uncertainty. The sampling of the parameters can be done via MC methods, such as MCMC (Berg 2004) and transitional MCMC (Ching and Chen 2007).

Practically, it can also be viewed as an ensemble approach, where each realization of the parameters is a member of the ensemble. Therefore, for the UQ methods discussed in Section 3, this UP method can be directly applied to models trained with MC dropout and the ensemble approach, but not for the others.

Using an NN interatomic potential trained with MC dropout in MD simulations, Wen and Tadmor (2020) adopted this sampling-based approach to propagate the uncertainty in atomic forces to the mechanical stress on monolayer graphene (Figure 13). As expected, the uncertainty in stress increases as the graphene layer is compressed or stretched from its equilibrium lattice parameter of 2.466 Å.

Figure 13:

Uncertainty propagation in molecular dynamics computation of mechanical stresses in a monolayer graphene. Error bars indicate uncertainty level. Adapted from Wen and Tadmor (2020) with a CC BY 4.0 license.

Recently, various studies have investigated the effects of uncertainty from ML models on microkinetic modeling. Li et al. (2023) was able to propagate the uncertainty in stoichiometric coefficients and parameters in the Arrhenius law (Arrhenius 1889a,b) to the concentration of chemical species in biodiesel production reaction systems, using the so-called Bayesian chemical reaction NNs. In a microkinetic study of ethanol steam reforming reactions, Xu and Yang (2023) investigated the propagation of errors in binding energy predicted by ML models to kinetic properties such as reaction rates. Their results demonstrated that the preferred reaction pathway varies depending on the used ML model.

5.2 Linearized propagation

Linearized UP is a general approach that can be used together with all UQ methods discussed in Section 3. Let δ _y be the uncertainty associated with the output y of an ML model y = f(x; θ). Then for a QoI z that is a function of the model output, z = g(y), we do a first-order Taylor expansion at y₀ to obtain z = g ( y 0 ) + ∂ g ∂ y ( y − y 0 ) . With this linearization, the uncertainty in z can be expressed as (Arras 1998):

meaning that the uncertainty in y can be propagated to z by multiplying the gradient. In general, if g is a function of multiple independent inputs, i.e., z = g(y₁, y₂, …, y _M ), the uncertainty in z can be written as (Arras 1998):

Once the uncertainty is obtained from a UQ method, it can be readily propagated to the QoI z, using Eq. (37) or Eq. (38). However, there are two challenges in applying this approach in practice. First, it may not be immediately obvious how to compute the gradients of g, which typically represent physics-based simulation techniques such as MD and microkinetic modeling. Second, while the approach is exact for a linear function g, it can introduce significant errors when applied to a nonlinear function (Cho et al. 2015).

The linearized UP approach becomes very attractive if the challenges are overcome. For example, Wen and Tadmor (2020) reformulated the integration algorithm in an MD simulation, and managed to propagate the uncertainty in atomic forces to stresses. As seen from Figure 13, the predicted mean stress and uncertainty agree very well with those from the sampling-based Bayesian approach. The linearized approach is computationally more efficient than the Bayesian approach, since it only needs one model evaluation, while the Bayesian approach requires multiple model evaluations.

5.3 Sensitivity analysis

Sensitivity analysis examines how the uncertainty in a model’s output can be apportioned to various sources of uncertainty in its inputs. Ideally, uncertainty and sensitivity analyses should be conducted in tandem; by working together, they can pinpoint the most critical input parameters, offering crucial insights into the model’s reliability and providing guidance for further model refinement. Although these techniques have not been widely employed together with ML models yet, we expect this to happen soon in the near future. To illustrate their potential, we introduce some of their usage with classical non-ML models.

Given a QoI z = g(y), we can perturb the input by some amount Δy and observe the change in the output Δz. This change, Δz, can be interpreted as the propagated uncertainty if we set the uncertainty in y as Δy. Typically, the input y lies in an M-dimensional parameter space, i.e., y = [y₁, y₂, …, y _M ]; then, depending on how Δy is chosen, we get local sensitivity and global sensitivity. Local sensitivity refers to perturbing each individual parameter y _i and observing its effects on z separately. Global sensitivity refers to exploring the entire parameter space simultaneously, considering interactions between the parameters.

Sensitivity analysis is an integral part of microkinetic modeling (Motagamwala and Dumesic 2020). Local sensitivity analysis like the derivative-based technique (Döpking et al. 2018) and global sensitivity analysis such as the Sobol’ method (Sobol 2001) and the Morris method (Morris 1991) have been widely employed. For example, Bensberg and Reiher (2024) have recently leveraged these techniques to automatically refine structures, reaction paths, and energies in chemical reaction networks, and successfully identified a small number of elementary reactions and compounds that are essential for reliably describing the kinetics of the Eschenmoser–Claisen rearrangement reactions (Wick et al. 1964) of allyl alcohol and of furfuryl alcohol. In addition, their Morris sensitivity analysis also provides the uncertainty in the predicted concentrations. Similarly, Kreitz et al. (2021) applied global uncertainty assessment and sensitivity analysis to explore parametric uncertainties in microkinetic models for CO₂ hydrogenation on the (111) surface of Ni. By systematically generating numerous mechanisms, they demonstrated how UQ can identify feasible models and optimize predictions within the uncertainty space.

Sensitivity analysis has also been applied to quantify the uncertainty in QoI obtained from MD simulations. Information-theoretic approaches provide a powerful framework for this purpose (Kurniawan et al. 2022). For a QoI z that can be obtained from an MD simulation, an upper bound for the uncertainty in z can be obtained as (Dupuis et al. 2016; Pantazis and Katsoulakis 2013; Tsourtis et al. 2015):

upon parameter perturbation Δθ, where Var denotes the variance, and I ( θ ) is the Fisher information, which measures the amount of information that the trained model carries about its parameter θ (Cover and Thomas 2012; Wen 2019). This provides an efficient way to investigate the reliability of MD simulations. Using this approach, Wen et al. (2017) studied the thickness of an MoS₂ sheet and found that Eq. (39) provides a tight bound, demonstrating high reliability of the MD predictions. Although Fisher information can provide useful insights into the uncertainty in MD simulations, the overall analysis is restricted to the perturbations only in the vicinity of the equilibrium model parameters.

5.4 Data uncertainty and propagation

The discussions in the above sections have primarily focused on uncertainties and their propagation in ML model. However, uncertainty is also inherent in chemical and materials data itself. For example, for data generated using DFT, the commonly employed semi-local density functional approximations have well-known intrinsic errors (Cohen et al. 2008; Perdew and Zunger 1981) that affect the accuracy and reliability of the results. When data is derived from a single semi-local density functional approximation, there is no effective way to address data uncertainty, leaving its impact on downstream properties often unaddressed. Below, we discuss two approaches that apply systematic statistical analysis to assess data uncertainty and its propagation to downstream properties, particularly in the context of microkinetic modeling.

The first method, developed by Wellendorff et al. (2012), is an ensemble-based approach that uses the same density functional but introduces tunable parameters. This method builds on the BEEF (Bayesian error estimation functional) family (Mortensen et al. 2005), incorporating an exchange-correlation functional (BEEF-vdW) with non-local correlation terms and tunable parameters to accurately capture interactions such as van der Waals forces, which are critical for surface science and catalysis (Wellendorff et al. 2012). Uncertainty estimates are then derived from an ensemble of functionals generated by parameter perturbation, with the ensemble’s standard deviation defining the uncertainty.

The second method, proposed by Walker et al. (2016), is also an ensemble-based approach. However, rather than perturbing a single functional, it employs multiple density functionals from different rungs of Perdew’s “Jacob’s Ladder” (Perdew and Schmidt 2001). In this approach, the energies of intermediate and transition states are calculated using various density functionals chosen for their applicability to the target system. A latent model using factor analysis (Rencher and Christensen 2012), is then developed to capture shared predictions across functionals and their unique variations. Finally, the predicted energies from the latent model are refined using a secondary probabilistic model to ensure consistency with reference data.

For both methods, once the corrected energies are obtained, they can be used in calculating downstream QoIs, and the associated uncertainty in the energies can be propagated using MC sampling as in Eq. (36). Using this approach, the first method has been widely applied in computational catalysis, from assessing reaction rate reliability (Lu et al. 2022) to identifying reaction mechanisms (Kreitz et al. 2023) and analyzing surface coverage (Wang et al. 2019). The second method has been successfully applied to determine reaction pathways in catalysis, such as identifying dominant mechanisms in the water-gas shift reaction through uncertainty propagation to turnover frequency calculations (Fricke et al. 2022; Walker et al. 2016, 2018).

While both methods can provide uncertainty in the data, the limitations should be noted. The methodology overall relies on fitting to experimental reference data or computed data with higher accuracy, such as the G3/99 dataset of experimental molecular formation energies (Curtiss et al. 2000) and the S22 dataset of intermolecular interaction energies (Jurečka et al. 2006). The reference data themselves, however, carry uncertainty. Such uncertainty is often assumed but not guaranteed to be negligible compared to DFT calculations (Wang et al. 2021). Additionally, the density functionals used or developed here are fitted to specific material properties, and thus their accuracy may not be transferable to other materials or properties (Mardirossian and Head-Gordon 2017; Medvedev et al. 2017). Thus, if the functionals in the ensemble perform poorly for a particular material or property, the resulting uncertainty, and any conclusions based on them, can be unreliable. This suggests that the pursuit of widely applicable density functionals, supported by thorough benchmark studies (Kim et al. 2024; Sheldon et al. 2021, 2024; Szaro et al. 2023), is equally vital for advancing methods to improve data uncertainty management.