Authors: Vardaan Tekriwal, Mark D. Risser, Hengrui Luo, Marcus M. Noack
Organizations: Applied Math & Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, CA · Department of Electrical Engineering & Computer Sciences, University of California, Berkeley, CA · Engineering Mathematics & Statistics Division, University of California, Berkeley, CA · Department of Industrial Engineering & Operations Research, University of California, Berkeley, CA · Climate & Ecosystem Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA · Department of Statistics, Rice University, Houston, TX
Conditional density estimation is complicated by multimodality, heteroscedasticity, and strong non-Gaussianity. Gaussian processes (GPs) provide a principled nonparametric framework with calibrated uncertainty, but standard GP regression is limited by its unimodal Gaussian predictive form. We introduce the Generalized Gaussian Mixture Process (GGMP), a GP-based method for multimodal conditional density estimation in settings where each input may be associated with a complex output distribution rather than a single scalar response. GGMP combines local Gaussian mixture fitting, cross-input component alignment and per-component heteroscedastic GP training to produce a closed-form Gaussian mixture predictive density. The method is tractable, compatible with standard GP solvers and scalable methods, and avoids the exponentially large latent-assignment structure of naive multimodal GP formulations. Empirically, GGMPs improve distributional approximation on synthetic and real-world datasets with pronounced non-Gaussianity and multimodality.
Figures & tables
Figure 1: Generalized Gaussian Mixture Processes (GGMPs) capture multimodal conditional distributions while preserving closed-form inference. Many real-world processes exhibit conditioning-dependent multimodality, asymmetry, and heteroscedasticity that standard GPs cannot represent. GGMPs address this by modeling p(y∣x) as a weighted mixture of GPs. Here, predictive densities at four held-out inputs are shown against empirical conditional distributions (blue histograms). With the number of components in the mixture K=1 (orange), GGMPs reduce to a standard heteroscedastic GP. Increasing K progressively resolves multimodal structure, with the model remaining robust to overspecification as surplus components overlap with existing ones.
Figure 2: Synthetic distribution field with fixed-point slices. The background shows log10p(y∣x) , generated from a single latent mean function f(x) (black curve), while red vertical slices depict sampled local conditional densities at selected x -locations. Despite sharing the same underlying f(x) , the local distributions vary strongly in modality, spread, and asymmetry.
Model
Bhattacharyya
Symmetric KL
Wasserstein-1
L1
GP1
0.4869 ± 0.1513
4.2005 ± 1.7752
0.7589 ± 0.2126
1.3115 ± 0.1699
MDN1
0.2278 ± 0.2124
2.9952 ± 3.0088
0.3223 ± 0.3369
0.7633 ± 0.4602
MDN3
0.0169 ± 0.0049
0.1561 ± 0.0486
0.0429 ± 0.0357
0.1989 ± 0.0586
MDN5
0.0163 ± 0.0048
0.1606 ± 0.0533
0.0366 ± 0.0309
0.1936 ± 0.0522
MDN10
0.0156 ± 0.0038
0.1316 ± 0.0382
0.0470 ± 0.0415
0.1925 ± 0.0462
MDN25
0.0179 ± 0.0028
0.1357 ± 0.0446
0.0296 ± 0.0244
0.1839 ± 0.0282
Table 1: Divergence metrics for the synthetic dataset.
Model
Log Score
CRPS
PIT Mean
PIT Std
Cov 50%
Cov 90%
Cov 95%
GP1
-1.6003
1.4150
0.4774
0.2681
0.5648
0.8854
0.9330
MDN1
-0.9420
0.5585
0.4822
0.2779
0.5115
0.9143
0.9570
MDN3
-0.3572
0.4935
0.5004
0.2290
0.6710
0.9762
0.9890
MDN5
-0.3505
0.4945
0.4914
0.2286
0.6700
0.9767
0.9884
MDN10
-0.3458
0.4943
0.4988
0.2300
0.6679
0.9755
0.9886
MDN25
-0.3415
0.4934
0.5029
0.2289
0.6713
0.9766
0.9892
Table 2: Predictive scoring and calibration diagnostics on the synthetic dataset.
Model
Bhattacharyya
Symmetric KL
Wasserstein-1
L1
GP1
0.1580 ± 0.1050
2.6399 ± 1.6705
2.2896 ± 0.9970
0.6090 ± 0.1987
MDN1
0.1480 ± 0.1033
2.6022 ± 1.0942
1.7232 ± 1.0110
0.6011 ± 0.1801
MDN3
0.1428 ± 0.1011
2.5923 ± 1.2212
1.2998 ± 1.0232
0.5801 ± 0.2336
MDN5
0.1414 ± 0.1006
2.5834 ± 1.6252
1.2121 ± 0.9547
0.5244 ± 0.2482
MDN10
0.1425 ± 0.0989
2.5901 ± 0.9998
1.0829 ± 1.1029
0.5328 ± 0.2434
MDN25
0.1382 ± 0.1001
2.5602 ± 1.4353
1.0021 ± 0.9887
0.4932 ± 0.2223
Table 3: Divergence metrics for the temperature dataset.
Model
Log Score
CRPS
PIT Mean
PIT Std
Cov 50%
Cov 90%
Cov 95%
GP1
-3.7951
13.7364
0.5134
0.3050
0.4194
0.9065
0.9624
MDN1
-3.7571
7.3935
0.5201
0.4387
0.4223
0.5997
0.6600
MDN3
-3.6779
6.1454
0.5030
0.3138
0.4868
0.7547
0.8015
MDN5
-3.6770
6.0630
0.4993
0.3055
0.4889
0.8013
0.8553
MDN10
-3.6699
5.7824
0.5001
0.3043
0.4989
0.8892
0.9342
MDN25
-3.6822
5.2364
0.4817
0.3570
0.4518
0.8438
0.8911
Table 4: Predictive scoring and calibration diagnostics on the temperature extremes dataset.
Figure 3: Joint density reconstruction on held-out conditions (multivariate, K=5 ). Blue contours show the empirical test distribution in the 2D output space (Axis 1–Axis 2), and red contours show the GGMP 5 predictive distribution at the same inputs. Each panel corresponds to a different held-out condition. Overall alignment of contour shape and multimodal structure indicates that the GGMP captures the major geometry of the joint distribution.
Table 8: Predictive scoring and calibration diagnostics on the additive manufacturing dataset.
Dataset
GGMP K
Ldens(w1/K)
Ldens(w∗)
Δw1/K→w∗
Synthetic
3
− 381.08
− 380.64
0.12%
5
− 345.32
− 345.18
0.04%
10
− 307.32
− 306.12
0.39%
25
− 301.29
− 299.62
0.55%
Temperature
3
− 22,379.2
− 22,378.1
< 0.01%
5
− 22,306.1
− 22,304.9
< 0.01%
Table 9: Effect of weight optimization on the distributional log-likelihood. Δw1/K→w∗ is the relative improvement (Lopt−Lequal)/∣Lequal∣×100 .
K
L(w=1/K)
L(w∗)
L(w∗(x))
Δw∗→w∗(x)
3
− 1,937,212
− 1,881,962
− 1,877,461
0.24%
5
− 1,966,832
− 1,917,670
− 1,902,757
0.78%
10
− 1,917,485
− 1,878,141
− 1,868,909
0.49%
25
− 2,143,543
− 2,101,179
− 2,088,154
0.62%
Table 10: Distributional log-likelihood on the additive-manufacturing dataset under equal weights, optimized shared weights, and optimized input-dependent weights ( 24 ). Δw∗→w∗(x) is the relative improvement from shared to input-dependent weights, (L(w∗(x))−L(w∗))/∣L(w∗)∣×100 .
Appendix figures & tables4 assets
Supplementary material from the paper’s appendix.
Appendix
K
GMM + Align
GP Training
Weight Opt
Prediction
GGMP Total
MDN Total
1
<1
39
0
1
<41
1,318
3
9
64
3
1
77
1,241
5
21
79
9
1
110
1,212
10
50
132
9
1
192
1,178
25
144
242
12
1
399
1,198
Appendix
Table 11: Wall-clock runtimes (seconds) for the synthetic function experiment ( N=300 , N∗T=600,000 , p=1 ).
K
GMM + Align
GP Training
Weight Opt
Prediction
GGMP Total
MDN Total
1
<1
593
0
57
<651
26,613
3
177
4,842
28
58
5,105
26,452
5
329
7,150
108
63
7,650
27,130
10
717
13,317
206
70
14,310
26,828
25
1,482
20,935
1,059
82
23,558
27,164
Appendix
Table 12: Wall-clock runtimes (seconds) for the U.S. temperature extremes experiment ( N≈7,000 , N∗TGGMP≈50,000,000 , N∗TMDN=10,000,000 , p=1 ).
K
GMM + Align
GP Training
Weight Opt
Prediction
GGMP Total
MDN Total
w∗
w∗(x)
( w∗ )
1
<1
3
0
0
<1
<5
1,009
3
15
8
2
19
3
28
1,018
5
24
11
7
204
5
47
997
10
42
18
41
528
6
107
1,035
25
256
26
123
929
7
412
991
Appendix
Table 13: Wall-clock runtimes (seconds) for the additive manufacturing experiment ( N=24 , N∗T≈600,000 , p=2 ). Weight optimization is reported for shared weights ( w∗ ) and input-dependent weights ( w∗(x) ); the GGMP total uses w∗ .
Figure 4: Held-out distribution reconstruction across mixture complexity for the synthetic function. Each row shows a different number of mixture components; each column is one of four approximately equally spaced held-out inputs. In every panel, the ground-truth is shown in gray, and the model prediction is red. Increasing K improves local shape fidelity—especially multimodal structure—while low K underfits fine-scale features. When K<Ktrue , we witness an averaging phenomenon; and when K>Ktrue , we observe a distribution that looks like it could have come from a mixture with fewer components, indicating that components are likely stacked.
Gaussian Processes (GPs) are a powerful tool for modelling and quantifying uncertainty in functional relationships. However, they require practitioners to make a number of design decisions, such as the choice of the kernel and the observation model. Suboptimal choices can produce misspecified models that do not capture the underlying data generating process. We introduce Predictively Oriented Gaussian Processes (PrO-GPs), which treat predictive uncertainty as the primary inferential target and provide a robust alternative to standard GPs. Although direct computation of a PrO posterior for nonparametric models is intractable, we derive a reduced formulation and practical sampling scheme for efficient computation. Through synthetic and real data experiments, we show that PrO-GPs produce better calibrated predictive distributions under model misspecification compared to standard GP approaches.
Gaussian processes (GPs) offer a principled probabilistic model over functions, but exact inference is restricted to the linear-Gaussian regime. We establish an explicit equivalence between GPs and a class of linear diffusion models, recasting predictive sampling as an ODE with closed-form Gaussian dynamics and a likelihood-dependent guidance term that admits a simple Monte Carlo approximation. In the linear-Gaussian setting, we recover standard GP conditioning exactly; beyond conjugacy, the same machinery handles any conditioning statement admitting point-wise likelihood evaluation -- including non-linear physics, and, for the first time, natural language via large language models. Whitening isolates the irreducible non-Gaussian dynamics, minimising Wasserstein-2 transport cost and eliminating numerical stiffness. The result is a general-purpose GP inference scheme requiring no bespoke derivations. Together, these results provide a general mechanism for incorporating the full richness of real-world knowledge as conditioning information, opening a new frontier for the probabilistic modelling of real-world problems.
Henry Moss, Lachlan Astfalck, Thomas Cowperthwaite +5
Lancaster University · University of New South Wales · University of Cambridge +1
We present a theoretically grounded Gaussian process framework that leverages neural feature maps to construct expressive kernels. We show that the learned feature map can be interpreted as an optimal low-rank approximation to a Gram matrix derived from an implied RKHS, from which we establish consistency of the GP posterior. We further analyse the spectral properties of the induced kernels and introduce product feature-map kernels to address oversmoothing. This simple yet powerful approach enables fast, scalable, and accurate exact GP inference with minimal upfront work. The flexibility of kernel design supports seamless application to both regression and classification tasks across diverse data modalities, including tabular inputs and structured domains such as images. On benchmark datasets, this approach surpasses pre-existing methods in terms of accuracy and training and prediction efficiency.