cs.LGAug 29, 2026
SaveInformation-Based Calibration of Uncertainty Quantification in Product-of-Experts Gaussian Process Models
Organizations: University College London, London, UK · The University of Manchester, Manchester, UK
Abstract
Gaussian process (GP) regression with a single global GP (GP-glo) incurs cubic computational cost, limiting scalability to large datasets. Product-of-experts GP models (GP-pro), which combine local GP models to capture global correlations, alleviate this computational burden. However, training local experts on disjoint data subsets can lead to overestimated posterior variances. We propose GP-pro-c, a product-of-experts GP model that calibrates these variances using an information-based method. The method exploits the monotonicity and submodularity of information gain in GPs to define a calibration ratio that reduces the posterior variance of individual local GP models. We evaluate GP-pro-c using negative log-likelihood (NLL), root mean squared error (RMSE), and expected normalised calibration error (ENCE). Experiments on four synthetic functions and six regression datasets show that GP-pro-c achieves average reductions of 2.3% in NLL and 12.0% in ENCE compared with the uncalibrated GP-pro model. The proposed method mitigates posterior variance overestimation while maintaining predictive accuracy and reducing computational complexity. GP-pro-c provides a promising approach for uncertainty estimation in scalable GP models and may serve as a useful surrogate model for Bayesian optimisation with high-dimensional and large-scale data.
Figures & tables
Figure 1. Illustration of posterior predictions produced by a single global GP model (GP-glo) and a product-of-experts GP model (GP-pro) for the one-dimensional Ackley function. The white dots represent the training data points. The plot on the left shows the posterior prediction of the GP-glo model, with the shaded region indicating the standard deviation at each input value . The plot on the right shows a GP-pro model comprising two local experts. The balltree assignment method is used to allocate training data points to the local GP models. The green and red curves, together with their corresponding shaded regions, represent the posterior predictions and standard deviations of the two local GP models. The GP-pro model employs an entropy-weighted aggregation scheme to combine the predictions of the local GP models. The blue curve and shaded region represent the aggregated posterior prediction and its associated standard deviation. Compared with the GP-glo model, the GP-pro model exhibits overestimated posterior variance, as evidenced by the wider blue uncertainty intervals. Two line graphs (left and right) show the Ackley function values, with the y-axis ranging from −4 to 2 and the x-axis representing input values from −10 to 10 in increments of 5. The left graph presents the function values and their corresponding confidence intervals predicted by a single global Gaussian process (GP) model. The right graph presents the function values and their corresponding confidence intervals predicted by a product of two local GP models. The confidence intervals produced by the product of local GP models are noticeably wider than those of the single global GP model, illustrating the issue of variance overestimation that can arise in product-of-Gaussian-process models.
| Dataset | Size | Dim | GP-pro-bt-ent | GP-pro-full | GP-full | ||||||
| RMSE | NLL | ENCE | RMSE | NLL | ENCE | RMSE | NLL | ENCE | |||
| Ackley | 500 | 1 | 0.300 | 0.267 | 0.296 | 0.290 | 0.204 | 0.252 | 0.290 | 0.204 | 0.252 |
| Ackley | 5000 | 10 | 0.487 | 0.775 | 0.245 | 0.468 | 0.663 | 0.057 | 0.468 | 0.663 | 0.057 |
| Levy | 5000 | 10 | 0.681 | 1.021 | 0.160 | 0.653 | 0.982 | 0.137 | 0.653 | 0.982 | 0.137 |
| Rastrigin | 5000 | 10 | 0.745 | 1.138 | 0.127 | 0.713 | 1.081 | 0.039 | 0.713 | 1.081 | 0.039 |
| Average | 0.553 | 0.800 | 0.207 | 0.531 | 0.732 | 0.121 | 0.531 | 0.732 | 0.121 | ||
Table 1. Average RMSE, NLL, and ENCE results for the GP-pro, GP-pro-full, and GP-full models on four synthetic functions. GP-pro-full and GP-full yield identical results, as expected, since GP-pro-full trains each expert on the entire training dataset, resulting in identical posteriors whose product recovers the global GP solution. GP-pro exhibits higher NLL and ENCE than GP-pro-full and GP-full, indicating variance overestimation in the GP-pro model.
Figure 2. Reliability diagrams for the GP-pro, GP-pro-full, and GP-full models on four synthetic functions. Well-calibrated models yield curves closer to the ideal diagonal line. The curves of GP-pro-full overlap with those of GP-full, as expected. For the Levy function, GP-full also exhibits noticeable deviation from the diagonal line, which can be attributed to kernel misspecification and the high-dimensional, nonstationary nature of the function. In contrast, GP-pro shows consistently larger deviations across all functions, indicating systematic overestimation of posterior variances in the GP-pro models. Four line graphs, labelled (a), (b), (c), and (d), show reliability diagrams for three types of Gaussian process models evaluated on the one-dimensional Ackley function, the 10-dimensional Ackley function, the 10-dimensional Levy function, and the 10-dimensional Rastrigin function, respectively. In each graph, the horizontal axis represents the Root Mean Variance (RMV), ranging from 0 to 1 in increments of 0.2, and the vertical axis represents the Root Mean Square Error (RMSE), also ranging from 0 to 1 in increments of 0.2. Each graph contains four lines: an ideal dashed diagonal line extending from the lower-left corner to the upper-right corner, and three model-specific lines corresponding to (i) a product-of-experts Gaussian process model, (ii) a product-of-experts Gaussian process model in which each local model was trained on the entire training dataset, and (iii) a single global Gaussian process model. Across all four graphs, the single global Gaussian process model produces a line that is closest to the ideal diagonal, indicating better calibration. In contrast, the product-of-experts Gaussian process model produces a line that deviates noticeably from the diagonal, reflecting poorer calibration.
| Dataset | Size | Dim | GP-pro | GP-pro-full | GP-full |
| training, test time | training, test time | training, test time | |||
| Ackley | 500 | 1 | 0.666, 0.018 | 0.677, 0.021 | 0.569, 0.013 |
| Ackley | 5000 | 10 | 3.101, 1.425 | 963.947, 13.745 | 47.905, 0.885 |
| Levy | 5000 | 10 | 3.380, 1.450 | 982.690, 14.161 | 46.986, 0.892 |
| Rastrigin | 5000 | 10 | 3.119, 1.452 | 972.641, 14.850 | 50.469, 1.053 |
| Average | 2.566, 1.086 | 729.988, 10.694 | 36.482, 0.711 |
Table 2. Average training and test times (in seconds) across four synthetic functions for the GP-pro, GP-pro-full, and GP-full models.
Figure 3. Illustration of posterior predictions before and after applying the calibration ratio to GP-pro-c models comprising two local experts for the one-dimensional Ackley function. The GP-pro-c model on the left uses the random assignment method, whereas the model on the right uses balltree assignment. Both GP-pro-c models employ an entropy-weighted aggregation scheme. The posterior predictions of the two local GP models are represented by the green and red curves, with the corresponding light green and light red shaded regions indicating the standard deviation at each input value . The blue curve and shaded region represent the aggregated posterior prediction. The darker green, red, and blue shaded regions indicate the reduced standard deviations obtained after applying the calibration ratio. Two line graphs, shown side by side, illustrate posterior predictions for the one-dimensional Ackley function. In both graphs, the y-axis shows function values ranging from -4 to 2, and the x-axis shows input values from -10 to 10 in increments of 2. Each graph contains three prediction curves and their corresponding uncertainty intervals. The first and second curves represent the posterior mean predictions of two local Gaussian process models, together with their uncertainty intervals before and after applying the calibration ratio. The third curve represents the aggregated posterior prediction obtained by combining the outputs of the two local Gaussian process models, along with its corresponding uncertainty interval. The darker shaded regions indicate the calibrated uncertainties after applying the calibration ratio, while the lighter shaded regions indicate the original uncertainties before calibration.
| Method | Training | Test-time | Test-time | Notes |
| complexity | mean | variance | ||
| Global GP | Full dataset, accurate variance | |||
| Vanilla GP-pro ( experts, points each) | Distributed experts, standard PoE | |||
| GP-pro-c (proposed) | Adds calibrated variance via with negligible cost | |||
Table 3. Computational complexity comparison.
| Model | Data assignment method | Weighting scheme | Uncertainty calibration | |
| bt-ent | balltree | entropy | no | |
| bt-var | balltree | variance | no | |
| bt-uni | balltree | uniform | no | |
| kx-ent | -means | entropy | no | |
| kx-var | -means | variance | no | |
| kx-uni | -means | uniform | no |
Table 4. Summary of the models used in the experiments for the comparison and evaluation of the GP-pro-c and GP-pro models.
Table 5. Average RMSE values across four synthetic functions and six benchmark datasets using balltree, -means, and random assignment methods, as well as entropy-, variance-, and uniform-weighted schemes for the GP-pro-c and GP-pro models. GP-pro-c models with information-based calibration achieve RMSE values comparable to those of GP-pro models, indicating that information-based calibration improves uncertainty quantification without affecting predictive accuracy.
Table 6. Average NLL values across four synthetic functions and six benchmark datasets using balltree, -means, and random assignment methods, as well as entropy-, variance-, and uniform-weighted schemes for the GP-pro-c and GP-pro models. GP-pro-c models with information-based calibration achieve lower NLL values than GP-pro models, demonstrating that information-based calibration improves uncertainty quantification without affecting predictive accuracy.
Table 7. Average ENCE values across four synthetic functions and six benchmark datasets using balltree, -means, and random assignment methods, as well as entropy-, variance-, and uniform-weighted schemes for the GP-pro-c and GP-pro models. GP-pro-c models with information-based calibration achieve lower ENCE values than GP-pro models, demonstrating improved uncertainty quantification.
Figure 4. Reliability diagrams comparing GP-pro-c and GP-pro models using balltree, -means, and random assignment methods across four synthetic functions and six regression datasets. GP-pro-c consistently yields curves closer to the ideal diagonal than GP-pro, indicating improved uncertainty calibration. The degree of deviation from the diagonal varies across datasets: smoother or lower-dimensional problems (e.g., Concrete, Airfoil, and Power) show near-diagonal curves, whereas noisier or higher-dimensional problems (e.g., Levy, Rosenbrock, Parkinsons, and Kin40K) show larger residual deviations. Ten line graphs show reliability diagrams for six variants of product-of-experts Gaussian process models evaluated on the 10-dimensional Ackley function, 10-dimensional Levy function, 10-dimensional Rastrigin function, 50-dimensional Rosenbrock function, the 8-dimensional Concrete dataset, the 5-dimensional Airfoil dataset, the 20-dimensional Parkinsons dataset, the 4-dimensional Power dataset, the 8-dimensional Kin40K dataset, and the 9-dimensional Protein dataset. \parIn each graph, the horizontal axis represents the Root Mean Variance (RMV), ranging from 0 to 1 in increments of 0.2, and the vertical axis represents the Root Mean Square Error (RMSE), also ranging from 0 to 1 in increments of 0.2. \parEach graph contains seven lines: an ideal dashed diagonal line extending from the lower-left corner to the upper-right corner, and six model-specific lines corresponding to (i) a product-of-experts Gaussian process model using the balltree assignment method, an entropy-weighted aggregation scheme, and uncertainty calibration; (ii) the same model without uncertainty calibration; (iii) a product-of-experts Gaussian-process model using the -means assignment method, an entropy-weighted aggregation scheme, and uncertainty calibration; (iv) the same model without uncertainty calibration; (v) a product-of-experts Gaussian process model using random assignment, an entropy-weighted aggregation scheme, and uncertainty calibration; and (vi) the same model without uncertainty calibration. \parAcross all graphs, the calibrated product-of-experts Gaussian process model using the balltree assignment method generally produces a curve that lies closest to the ideal diagonal, indicating better uncertainty calibration than the other model variants. Curves corresponding to uncalibrated models tend to deviate more substantially from the diagonal.
Table 8. Average training and test times (in seconds) across four synthetic functions and six regression datasets using balltree, -means, and random assignment methods, as well as entropy- and variance-weighted schemes (the uniform-weighted scheme is omitted for brevity) for the GP-pro-c and GP-pro models. Variants of GP-pro-c have training and test times similar to those of GP-pro.
Figure 5. Average NLL, RMSE, ENCE, and training time for GP-pro-c models using balltree, -means, and random assignment methods with 100-600 data points per local GP, evaluated on four synthetic functions and four benchmark datasets. Random assignment yields lower ENCE but higher RMSE and NLL due to conservative, weakly correlated local GPs. Balltree and -means preserve locality, yielding lower RMSE and NLL but higher ENCE. -means performs poorly with small cluster sizes because small, high-dimensional clusters can be heterogeneous and poorly conditioned; performance improves sharply as cluster size increases. -means can also incur slightly higher computational overhead due to uneven cluster sizes, whereas balltree and random assignments are more balanced. Overall, balltree assignment with 200-400 points per local GP provides a good trade-off between accuracy, calibration, and computational cost. Four line graphs labelled (a), (b), (c), and (d). \parIn all graphs, the horizontal x-axis shows the number of data points per local Gaussian process, ranging from 100 to 600 in increments of 100. Each graph contains three lines representing product-of-experts Gaussian process models constructed using the balltree, -means, and random assignment methods, respectively. All models use the entropy-weighted scheme with uncertainty calibration. \parIn (a), the vertical y-axis shows the average Root Mean Square Error (RMSE), ranging from 0.40 to 0.60 in increments of 0.05. All three lines decrease consistently as the number of data points increases. The balltree line remains the lowest throughout. The -means line is initially the highest at 100 data points, but from 200 data points onward it becomes the second lowest, overtaking the random assignment line. \parIn (b), the vertical y-axis shows the average Negative Log-Likelihood (NLL), ranging from 0.2 to 1.2 in increments of 0.2. All three lines decrease consistently. The balltree line remains the lowest throughout. The -means line is initially the highest at 100 data points, but from 200 data points onward it becomes the second lowest, overtaking the random assignment line. \parIn (c), the vertical y-axis shows the average Expected Normalised Calibration Error (ENCE), ranging from 0.1 to 0.5 in increments of 0.1. The balltree and random assignment lines remain the lowest, at approximately 0.2 ENCE, throughout the graph. The -means line begins at 0.45 ENCE with 100 data points, then decreases sharply to approximately 0.25 ENCE at 200 data points and remains the highest thereafter. \parIn (d), the vertical y-axis shows the average training time in seconds, ranging from 10 to 35 in increments of 5. All three lines increase exponentially as the number of data points grows. The balltree and random assignment lines remain the lowest throughout, while the -means line remains the highest.
| Dataset | Size | Dim | GP-pro-bt-ent-c | GP-pro-bt-ent | GP-pro-sm | GP-full | ||||||||
| RMSE | NLL | ENCE | RMSE | NLL | ENCE | RMSE | NLL | ENCE | RMSE | NLL | ENCE | |||
| Ackley | 5000 | 10 | 0.487 | 0.742 | 0.190 | 0.487 | 0.775 | 0.245 | 0.527 | 0.801 | 0.115 | 0.468 | 0.663 | 0.057 |
| Levy | 5000 | 10 | 0.681 | 1.012 | 0.103 | 0.681 | 1.021 | 0.160 | 0.731 | 1.097 | 0.115 | 0.653 | 0.982 | 0.137 |
| Rastrigin | 5000 | 10 | 0.745 | 1.124 | 0.060 | 0.745 | 1.138 | 0.127 | 0.801 | 1.211 | 0.116 | 0.713 | 1.081 | 0.039 |
| Rosenbrock | 50000 | 50 | 0.610 | 0.918 | 0.116 | 0.610 | 0.929 | 0.167 | 0.711 | 1.106 | 0.162 | N/A | N/A | N/A |
| Concrete | 1030 | 8 | 0.336 | 0.287 | 0.162 | 0.336 | 0.287 | 0.166 | 0.341 | 0.298 | 0.186 | 0.335 | 0.282 | 0.190 |
Table 9. Average RMSE, NLL, and ENCE values across four synthetic functions and six regression benchmark datasets, comparing GP-pro-c models with uncertainty calibration, GP-pro models without calibration, GP-pro-sm models using softmax-weight calibration, and a single global GP (GP-full). GP-pro-c models achieve lower NLL and ENCE scores, indicating that the proposed information-based calibration maintains predictive accuracy while improving uncertainty calibration.
Figure 6. Reliability diagrams for the GP-pro-c, GP-pro, GP-pro-sm, and GP-full models on four synthetic benchmark functions and six regression datasets. The GP-full models consistently exhibit curves closest to the diagonal. The GP-pro-c models demonstrate improved calibration compared to GP-pro across most datasets. On the Kin40K dataset, the GP-pro-sm model exhibits a curve closest to the diagonal. However, on the Protein dataset, GP-pro-sm exhibits noticeable deviation from the diagonal, whereas GP-pro-c maintains more stable calibration behaviour across datasets. Ten line graphs showing reliability diagrams for six types of Gaussian process models on the 10-dimensional Ackley function, 10-dimensional Levy function, 10-dimensional Rastrigin function, 50-dimensional Rosenbrock function, 8-dimensional Concrete dataset, 5-dimensional Airfoil dataset, 20-dimensional Parkinsons dataset, 4-dimensional Power dataset, 8-dimensional Kin40K dataset, and 9-dimensional Protein dataset. \parOn each line graph, the horizontal x-axis shows Root Mean Variance ranging from 0 to 1 in increments of 0.2, and the vertical y-axis shows Root Mean Square Error ranging from 0 to 1 in increments of 0.2. \parEach plot contains five lines: a dashed diagonal line from the lower-left to the upper-right corner representing perfect calibration, and four curves corresponding to (i) a product-of-experts Gaussian process model using the balltree assignment method with entropy-weighted scheme and uncertainty calibration, (ii) the same model without uncertainty calibration, (iii) the same model with softmax weight calibration, and (iv) a single global Gaussian process model. \parIn the Ackley, Levy, Rastrigin, Concrete, Airfoil, and Parkinsons datasets, the single global Gaussian process model produces curves closest to the ideal diagonal line. \parAcross all datasets except Kin40K, the balltree-based product-of-experts Gaussian process model with entropy weighting and uncertainty calibration yields curves closer to the ideal diagonal line than the corresponding model with softmax weight calibration.
| Dataset | Size | Dim | GP-pro-bt-ent-c | GP-pro-bt-ent | GP-pro-sm | GP-full |
| training, test time | training, test time | training, test time | training, test time | |||
| Ackley | 5000 | 10 | 3.101, 1.427 | 3.101, 1.425 | 3.101, 1.425 | 47.905, 0.885 |
| Levy | 5000 | 10 | 3.380, 1.452 | 3.380, 1.450 | 3.380, 1.450 | 46.986, 0.892 |
| Rastrigin | 5000 | 10 | 3.119, 1.454 | 3.119, 1.452 | 3.119, 1.453 | 50.469, 1.053 |
| Rosenbrock | 50000 | 50 | 32.522, 85.483 | 32.522, 85.375 | 32.522, 85.695 | N/A |
| Concrete | 1030 | 8 | 1.263, 0.037 | 1.263, 0.036 | 1.263, 0.037 | 1.717, 0.033 |
Table 10. Average training and test times (in seconds) across four synthetic functions and six regression benchmark datasets for GP-pro-c models with uncertainty calibration, GP-pro models without calibration, GP-pro-sm models using softmax-weight calibration, and a single global GP (GP-full).
Appendix figures & tables1 asset
Supplementary material from the paper’s appendix.
Appendix
Figure 7. Reliability diagrams for the GP-pro-c and GP-pro models using the balltree assignment method, comparing entropy-, variance-, and uniform-weighted schemes on four synthetic benchmark functions and six regression datasets. Under all weighting schemes, GP-pro-c models with uncertainty calibration produce curves closer to the ideal diagonal line than GP-pro models without uncertainty calibration. Ten line graphs show reliability diagrams for Gaussian process models on the 10-dimensional Ackley function, 10-dimensional Levy function, 10-dimensional Rastrigin function, 50-dimensional Rosenbrock function, 8-dimensional Concrete dataset, 5-dimensional Airfoil dataset, 20-dimensional Parkinsons dataset, 4-dimensional Power dataset, 8-dimensional Kin40K dataset, and 9-dimensional Protein dataset. \parEach line graph has Root Mean Variance on the horizontal axis (ranging from 0 to 1 in increments of 0.2) and Root Mean Square Error on the vertical axis (ranging from 0 to 1 in increments of 0.2). \parEach plot shows seven curves: a dashed diagonal reference line and six model variants corresponding to combinations of balltree assignment with entropy-, variance-, and uniform-weighted schemes, each evaluated with and without uncertainty calibration. \parAcross all datasets, models with uncertainty calibration consistently produce curves closer to the ideal diagonal line than their uncalibrated counterparts.