stat.MLAug 13, 2026

On the Brittleness of Maximum Likelihood Estimation for Gaussian Process Hyperparameter Optimization

Authors: Tyler R. Johnson, Kian Ben-Jacob, Christopher P. Muller, Ramin Bostanabad

Organizations: Department of Mechanical and Aerospace Engineering, University of California, Irvine · Department of Civil and Environmental Engineering, University of California, Irvine

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

Machine learning (ML) has become an indispensable part of modern engineering design workflows. A crucial step in training an ML model is the selection of the loss function which can be systematically formulated via various techniques such as maximum likelihood estimation (MLE) and cross-validation . While MLE is one of the most popular, effective, and intuitive mechanisms for training ML models, it is brittle: if the assumptions underpinning it are not met, the trained ML model may generalize poorly. This brittleness affects even Gaussian processes (GPs) which are widely used in engineering design and are often (incorrectly) presumed to be very robust to overfitting. In this paper, we fundamentally evaluate the brittleness of MLE in the context of training GPs for probabilistic regression or classification tasks. We compare theoretically grounded metrics against MLE and propose practical solutions. Our extensive studies demonstrate the effectiveness of our solutions in downstream design tasks such as Bayesian optimization and provide a blueprint for practitioners to build accurate and robust GPs that can even outperform tabular foundation models in terms of prediction accuracy, uncertainty quantification, and inference cost. Our contributions are publicly available via GitHub at https://github.com/Bostanabad-Research-Group/GP-vs-TabPFN-vs-GPyTorch.

Figures & tables

Figure 1: 1D regression examples with n=20n{=}20 training points: Mean predictions with 95% prediction intervals are shown. Axes span x∈[−0.5,0.5]x\in[-0.5,0.5] , with tick labels omitted for compactness. The defining equations for each example are given in the project repository [ 54 ] .
Figure 2: One-dimensional classification examples following Milios et al. [ 48 ] , with two classes and gaps in the training data. Curves give the predicted probability P(y=1∣x)P(y{=}1\mid x) and shaded regions the 95% interval; filled markers denote training points of class 1 and open markers those of class 0. Class 0 probabilities are the complement of those shown and are omitted. Vertical axes show P(y=1∣x)P(y{=}1\mid x) and horizontal axes span x∈[−1,1]x\in[-1,1] , with xx tick labels omitted for compactness.
Figure 3: Evaluations based on RRMSE in Equation 14 : Noise is either ε∼N(0,(0.002×c)2)\varepsilon\sim\mathcal{N}(0,(0.002\times c)^{2}) or ε∼N(0,(0.08×c)2)\varepsilon\sim\mathcal{N}(0,(0.08\times c)^{2}) . GPs are used on CPU while TabPFN is used on GPU. All the GPs are trained with the default settings of GP+ and GPyTorch packages. GP+ (PE) and GP+ (LOO) denote cases where the default Gaussian kernel and loss function are switched to, respectively, the PE kernel and LOO-based loss in Equation 4 .
Figure 4: Evaluations based on NIS in Equation 7 : Noise is either ε∼N(0,(0.002×c)2)\varepsilon\sim\mathcal{N}(0,(0.002\times c)^{2}) or ε∼N(0,(0.08×c)2)\varepsilon\sim\mathcal{N}(0,(0.08\times c)^{2}) . GPs are used on CPU while TabPFN is used on GPU. All the GPs are trained with the default settings of GP+ and GPyTorch packages. GP+ (PE) and GP+ (LOO) denote cases where the default Gaussian kernel and loss function are switched to, respectively, the PE kernel and LOO-based loss in Equation 4 .
Figure 5: Evaluations based on NCRPS in Equation 10 : Noise is either ε∼N(0,(0.002×c)2)\varepsilon\sim\mathcal{N}(0,(0.002\times c)^{2}) or ε∼N(0,(0.08×c)2)\varepsilon\sim\mathcal{N}(0,(0.08\times c)^{2}) . GPs are used on CPU while TabPFN is used on GPU. All the GPs are trained with the default settings of GP+ and GPyTorch packages. GP+ (PE) and GP+ (LOO) denote cases where the default Gaussian kernel and loss function are switched to, respectively, the PE kernel and LOO-based loss in Equation 4 .
Figure 6: Evaluations based on RRMSE and NIS (Rosenbrock Dx=80D_{x}=80 ): Noise is either ε∼N(0,(0.002×c)2)\varepsilon\sim\mathcal{N}(0,(0.002\times c)^{2}) or ε∼N(0,(0.08×c)2)\varepsilon\sim\mathcal{N}(0,(0.08\times c)^{2}) . GPs are used on CPU while TabPFN is used on GPU. All the GPs are trained with the default settings of GP+ and GPyTorch packages. GP+ (PE) denotes the case where the default Gaussian kernel is switched to the PE kernel.
Figure 7: Evaluations based on RRMSE and NIS (Buckling Dx=4D_{x}=4 ) with and without log-scale target preprocessing: Noise is either ε∼N(0,(0.002×c)2)\varepsilon\sim\mathcal{N}(0,(0.002\times c)^{2}) or ε∼N(0,(0.08×c)2)\varepsilon\sim\mathcal{N}(0,(0.08\times c)^{2}) . GPs are used on CPU while TabPFN is used on GPU. All the GPs are trained with the default settings of GP+ and GPyTorch packages. GP+ (log) and GPyTorch (log) denote models trained after a log transform of the response, with metrics evaluated on the original scale.
Figure 8: Bayesian optimization performance for GP+ and TabPFN across benchmark problems with observation noise ε∼N(0,(0.08×c)2)\varepsilon\sim\mathcal{N}(0,(0.08\times c)^{2}) . Buckling and Borehole are maximization problems, hence best negative y values are reported.
Figure 9: Classification performance on the Electrical Grid Stability, Truss 6D, Stellar, Steel Plates Faults benchmarks for varying training set size NN . Columns correspond to the datasets in that order and rows to L\mathcal{L} , Accuracy, ECE, respectively. With the exception of using Adam as the optimizer, GP+ package default settings are used for training of all GPs. Error bars denote one standard deviation across 10 seeds.
Appendix figures & tables3 assets

Supplementary material from the paper’s appendix.

Appendix

    RRMSENIS
ProblemND xNoiseGP+GP+ (PE)GP+ (LOO)PFN 2.5PFN 2.0GPyTorchGP+GP+ (PE)GP+ (LOO)PFN 2.5PFN 2.0GPyTorch
Buckling5D x40.0020.344 ±\pm 0.140.344 ±\pm 0.140.519 ±\pm 0.180.385 ±\pm 0.260.286 ±\pm 0.140.730 ±\pm 0.156.940 ±\pm 3.206.940 ±\pm 3.2010.35 ±\pm 5.091.377 ±\pm 0.241.524 ±\pm 0.278.769 ±\pm 1.92
   0.080.537 ±\pm 0.120.537 ±\pm 0.120.851 ±\pm 0.570.459 ±\pm 0.110.423 ±\pm 0.140.759 ±\pm 0.179.781 ±\pm 3.189.781 ±\pm 3.1818.10 ±\pm 10.41.806 ±\pm 0.271.902 ±\pm 0.329.340 ±\pm 2.00
 20D x40.0020.058 ±\pm 0.020.058 ±\pm 0.020.049 ±\pm 0.020.346 ±\pm 0.270.205 ±\pm 0.250.261 ±\pm 0.070.198 ±\pm 0.060.198 ±\pm 0.060.235 ±\pm 0.160.354 ±\pm 0.100.476 ±\pm 0.151.438 ±\pm 0.58
   0.080.226 ±\pm 0.060.226 ±\pm 0.060.249 ±\pm 0.230.585 ±\pm 0.350.382 ±\pm 0.160.335 ±\pm 0.081.139 ±\pm 0.491.139 ±\pm 0.491.359 ±\pm 2.920.981 ±\pm 0.171.285 ±\pm 0.221.847 ±\pm 0.72
Borehole5D x80.0020.024 ±\pm 0.010.027 ±\pm 0.010.044 ±\pm 0.010.049 ±\pm 0.010.065 ±\pm 0.010.027 ±\pm 0.000.131 ±\pm 0.040.139 ±\pm 0.050.645 ±\pm 0.300.271 ±\pm 0.030.345 ±\pm 0.050.138 ±\pm 0.04
Appendix
Table A1 : RRMSE and NIS statistics across benchmark problems for GP+, GP+ (PE), GP+ (LOO), PFN 2.5, PFN 2.0, and GPyTorch. Entries are median ±\pm std over runs. Best models for a given problem (lowest median per metric) are in bold.
Figure A1 : Ranking of Models by NLL: These plots rank the training NLL of GP+, GP+ (PE), and GPyTorch on problems (a)-(i), where (a) is Buckling, (b) is Borehole, (c) is Wing Weight, (d) is Ackley Dx=20D_{x}=20 , (e) is Griewank Dx=20D_{x}=20 , (f) is Zakharov Dx=20D_{x}=20 , (g) is Ackley Dx=40D_{x}=40 , (h) is Dixon-Price Dx=40D_{x}=40 , and (i) is Rosenbrock Dx=80D_{x}=80 .
 GP+GP+ (PE)GP+ (LOO)PFN 2.5PFN 2.0GPyTorch
ProblemND xNoiseTrain (s)Pred (s)Train (s)Pred (s)Train (s)Pred (s)Fit (s)Inf (s)Fit (s)Inf (s)Train (s)Pred (s)
Buckling5D x40.0025.23 ±\pm 1.160.011 ±\pm 0.0024.55 ±\pm 1.010.011 ±\pm 0.0377.14 ±\pm 0.970.010 ±\pm 0.0030.17 ±\pm 0.071.05 ±\pm 0.060.17 ±\pm 0.021.23 ±\pm 0.041.54 ±\pm 0.300.006 ±\pm 0.004
   0.085.54 ±\pm 1.590.012 ±\pm 0.0024.79 ±\pm 1.350.011 ±\pm 0.0026.88 ±\pm 0.930.010 ±\pm 0.0020.16 ±\pm 0.151.05 ±\pm 0.040.15 ±\pm 0.031.21 ±\pm 0.061.63 ±\pm 0.490.008 ±\pm 0.002
Buckling20D x40.0026.81 ±\pm 1.370.016 ±\pm 0.0015.87 ±\pm 1.250.017 ±\pm 0.0036.15 ±\pm 1.360.015 ±\pm 0.0020.18 ±\pm 0.041.04 ±\pm 0.060.18 ±\pm 0.041.24 ±\pm 0.042.62 ±\pm 0.620.010 ±\pm 0.004
   0.082.92 ±\pm 1.980.017 ±\pm 0.0032.73 ±\pm 1.730.014 ±\pm 0.0026.03 ±\pm 2.150.015 ±\pm 0.0030.17 ±\pm 0.021.06 ±\pm 0.040.15 ±\pm 0.031.20 ±\pm 0.042.08 ±\pm 0.410.008 ±\pm 0.001
Borehole5D x80.0020.50 ±\pm 0.100.006 ±\pm 0.0020.68 ±\pm 0.080.013 ±\pm 0.0020.92 ±\pm 0.180.009 ±\pm 0.0020.17 ±\pm 0.021.21 ±\pm 0.060.19 ±\pm 0.041.43 ±\pm 0.052.47 ±\pm 0.210.006 ±\pm 0.001
Appendix
Table A2 : Training and inference time (median ±\pm std) for GP+, GP+ (PE), GP+ (LOO), PFN 2.5, PFN 2.0, and GPyTorch across benchmark problems.
More in this topic

Gaussian Processes

Gaussian Processes are probabilistic models of functions that use kernels to make predictions with explicit uncertainty.

See all 113 papers

Explore similar work

CardsList

  1. Sharp Asymptotic Theory of Maximum Likelihood Estimation for Gaussian Processes with an RBF Kernel

    Oct 7, 2026Ameer Qaqish, Sudipto Banerjee, Didong LiGaussian ProcessesKernel Methods

    Department of Biostatistics, University of North Carolina at Chapel Hill · Department of Biostatistics, University of California, Los Angeles

  2. Predictively Oriented Gaussian Process Posteriors

    Oct 2, 2026Callum Lau, Jeremias Knoblauch, Louis SharrockGaussian Process RegressionUncertainty Quantification

    University College London

  3. Goal-Oriented Lower-Tail Calibration of Gaussian Processes for Bayesian Optimization

    May 19, 2026Aurélien Pion, Emmanuel VazquezPost-Hoc CalibrationBayesian Optimization

    Transvalor S.A., Biot, France · Univ. Paris-Saclay, CNRS, CentraleSup´elec, L2S, Gif-sur-Yvette, France