Heteroscedastic GP Regression Overview

• Heteroscedastic GP Regression is defined by modeling input-dependent noise variance, enabling tailored uncertainty quantification across the feature space.
• Exact posterior moments are derived via a mixture representation and importance sampling, yielding precise predictions and risk-aware control applications.
• Efficient inference leverages Gaussian proposals and blockwise computations, ensuring scalability and improved calibration in real-world noisy settings.

Heteroscedastic Gaussian Process (GP) Regression generalizes standard GP regression to cases where the observation noise variance is input-dependent. Such models are essential for real-world scenarios where the sources of noise vary systematically across the feature space. Unlike homoscedastic GP regression, where the posterior remains Gaussian and all predictive moments are available in closed form, the heteroscedastic generalization introduces significant technical and computational complexity, particularly for exact Bayesian inference. Recent work (Ito, 19 Sep 2024) provides a rigorous theoretical framework that delivers exact posterior means, variances, and full distributional descriptions for heteroscedastic GP regression, and demonstrates direct applications to chance-constrained control.

1. Model Specification: Latent Structure and Priors

In heteroscedastic GP regression (HGPR), the observation process is defined as

$$y_{d,s} = f(x_d) + w_{d,s}, \qquad w_{d,s} \mid g(x_d) \sim \mathcal{N}(0,\,g(x_d)^2),$$

where $f: \mathcal{X} \subset \mathbb{R}^n \to \mathbb{R}$ is a latent function, $g(x)$ is the (positive, potentially unbounded) noise standard deviation at $x$, and $w_{d,s}$ are i.i.d. noise samples per input $x_d$, $s=1,\ldots, S$.

To enforce positivity of $g^2$, it is conventional to parameterize via the log-variance function: $h(x) = \ln g(x)^2.$ Independent Gaussian process priors are then assigned: $$f \mid x_{1:D} \sim \mathcal{N}(0, K_{D}), \qquad h \mid x_{1:D} \sim \mathcal{N}(0, L_{D}),$$ where $K_{D}$, $L_{D}$ are Gram matrices for positive-definite kernels $k$ and $\ell$. For instance, a squared-exponential kernel may be used: $$k(x, x') = \sigma_f^2 \exp\left(-\tfrac{1}{2}(x-x')^\top \Lambda^{-1} (x-x')\right),$$ with $\sigma_f>0$, $\Lambda$ diagonal.

This construction is fundamentally different from standard GP regression, as the noise variance is a (latent) function of the input rather than a global scalar hyperparameter. Other models related to this specification include latent covariate models (GPLC, (Wang et al., 2012)), composite GPs (Davis et al., 2019), models with latent precision mixtures (Lee et al., 2023), and explicit robustifications (Student-t noise, (Hartmann et al., 2017); contaminated normal likelihoods, (Iong et al., 27 Feb 2024)).

2. Exact Posterior Characterization and Non-Gaussianity

While $(f, h)$ is jointly GP-distributed, the process of observing heteroscedastic noise renders the marginal posterior $p(f(x)|D)$ non-Gaussian, as the likelihood couples $f$ and $h$ in a non-linear fashion. The exact marginal posterior at a new input $x$ is derived in two steps (Ito, 19 Sep 2024):

1. Conditional on $h$, the vector of sample means $\bar y_d = \frac{1}{S} \sum_{s=1}^S y_{d,s}$ follows

$$\bar y_d \mid f(x_d), h_d \sim \mathcal{N}(f(x_d),\, S^{-1} e^{h_d}).$$

The posterior for $f(x)$ given $h_{1:D}$ is thus multivariate normal, with mean and variance: \begin{align*} m(x\,|\,h_{1:D}) &= k_D(x)^\top \Sigma(h_{1:D})^{-1} \bar y_{1:D}, \ v(x\,|\,h_{1:D}) &= k(x,x) - k_D(x)^\top \Sigma(h_{1:D})^{-1} k_D(x), \end{align*} where $\Sigma(h_{1:D}) = K_D + S^{-1} H(h_{1:D})$; $$H(h_{1:D}) = \operatorname{diag}(e^{h_1}, \ldots, e^{h_D})$$.

2. Marginalizing over $h$ using the law of total expectation and variance, the posterior mean and variance of $f(x)$ are \begin{align*} \mu(x|D) &= \mathbb{E}h [m(x|h)], \ \sigma^2(x|D) &= \mathbb{E}_h [v(x|h)] + \operatorname{Var}_h [m(x|h)]. \end{align*} Explicitly, these are written as \begin{align*} \mu(x|D) &= k_D(x)^\top\, \mathbb{E}_h[\Sigma(h)^{-1} ]\, \bar y{1:D}, \ \sigma^2(x|D) &= k(x,x) + k_D(x)^\top \Bigl( \mathbb{E}_h [\Sigma(h)^{-1} \bar y \bar y^\top \Sigma(h)^{-1}] - \mathbb{E}_h [\Sigma(h)^{-1}]\, \bar y \bar y^\top\, \mathbb{E}_h[\Sigma(h)^{-1}] - \mathbb{E}_h [\Sigma(h)^{-1}] \Bigr) k_D(x). \end{align*} Since $p(h_{1:D}|D)$ is non-Gaussian, these expectations are generally analytically intractable—$p(f(x)|D)$ becomes a mixture of Gaussians parameterized by the latent $h_{1:D}$. Posterior CDFs take the form

$$P(f(x)\leq y \mid D) = \mathbb{E}_{h|D}\left[ \Phi\left( \frac{y - m(x|h)}{\sqrt{v(x|h)}} \right) \right ],$$

a nontrivial mixture of normal CDFs rather than a single Gaussian CDF.

This structure underpins the key insight that in HGPR the posterior is not closed under Gaussianity; all predictive quantities require expectation over a non-Gaussian $h$ posterior.

3. Inference and Approximation: Importance Sampling

Modern HGPR frameworks must efficiently approximate expectations over $h$. The importance sampling proposal of (Ito, 19 Sep 2024) is as follows:

1. Select a Gaussian proposal $q(h_{1:D}) = \mathcal{N}(\mu_q, \Sigma_q)$ with moments matched to the likelihood of the sample variances.
2. Draw $M$ i.i.d. samples $h^{(m)}$ from $q$.
3. Compute normalized importance weights

$$w_m \propto \frac{p(D\,|\,h^{(m)})\,p(h^{(m)})}{q(h^{(m)})}, \qquad \sum_{m=1}^M w_m = 1.$$

1. Use the mixture to estimate arbitrary moments as

$$\mathbb{E}[\varphi(h)] \approx \sum_{m=1}^M w_m\, \varphi(h^{(m)}).$$

This approach yields consistent Monte Carlo estimates for all moments of interest, under mild regularity conditions ($S>1$ per design site). Deterministically, the cost per prediction is $O(MD^3)$ (from repeated matrix inversions at each $h^{(m)}$), which in practice is computationally viable for moderate $D$ and moderate ensemble size $M$.

Alternative approximate inference strategies include Laplace approximations or natural gradient algorithms (critical for non-log-concave likelihoods as in Student-t noise; see (Hartmann et al., 2017)), variational methods for large-scale problems (Liu et al., 2018, Lee et al., 2023), or MCMC samplers (see (Davis et al., 2019) for full Bayesian composite GPs).

4. Theoretical Results and Downstream Inference

A key theoretical advance (Ito, 19 Sep 2024) is the closed-form derivation of all available moments of the HGPR posterior and predictive distribution. The mixture representation enables calculation of the exact mean, variance, and cumulative distribution functions, all of which are crucial for risk-aware decision-making.

These exact expressions can be directly utilized in chance-constrained control. In discrete-time systems with components governed by potentially unknown input-dependent disturbances, HGPR predictions of $f(x)$ with full uncertainty quantification permit explicit enforcement of constraints such as

$$P\left( \left| x_{t+1}^{(2)} - r_{t+1}^{(2)} \right| \leq \epsilon \,\Big|\, x_t, D \right ) \geq 1 - \delta,$$

where the set of admissible actions corresponds to quantiles of the (non-Gaussian) HGPR posterior. The method involves numerically inverting the estimated posterior CDF—computed via the importance-sampled mixture structure—for each desired quantile, permitting explicit quantile control under input-dependent stochasticity.

This property is critical in safety-critical control and reinforcement learning settings, as homoscedastic or Gaussian-only approximations can either under- or over-estimate the actual risk envelope.

5. Comparison with Related Heteroscedastic GP Paradigms

Multiple modeling strategies have been developed:

Approach	Input-dependent variance	Outlier robustness	Posterior Gaussianity
GP with GP prior on log-variance ((Ito, 19 Sep 2024), etc)	Yes	Limited	No (mixture of Gaussians)
GPLC, latent covariate (Wang et al., 2012)	Yes (via latent $w$)	Non-Gaussian noise	No (via latent integration)
Heteroscedastic Student-t (Hartmann et al., 2017)	Yes	Heavy-tailed/hard	Approximate (Laplace/MCMC)
Composite GP (Davis et al., 2019)	Yes	No	No
Mixture-precision GP (Lee et al., 2023)	Yes (mixture, matrix)	Yes (extensible)	No (variational/EM approx)

A critical insight is that any model with a nontrivial latent variance (or even arbitrary non-Gaussian, input-dependent noise model) yields intractable or non-Gaussian posteriors—one must always resort to approximations or sampling for full distributional inference.

Recent developments in scalable inference (e.g., sparse variational (Liu et al., 2018), distributed experts (Liu et al., 2018), stochastic variational EM (Iong et al., 27 Feb 2024)) enable HGPR to be used on massive datasets while retaining robust uncertainty quantification.

6. Practical Considerations and Empirical Behavior

HGPR delivers both mean predictions and predictive variances that adapt locally to the data density and observed noise structure. This behavior is contrasted with homoscedastic GP predictions, which often underestimate uncertainty in “noisy” regions and overestimate it in “quiet” or clamped regions.

For observed datasets with replications per design site, blockwise Woodbury identities yield dramatic computational savings for likelihood calculations and predictive moments (Binois et al., 2016). Real-world simulation, inventory, and epidemics management datasets have demonstrated that HGPR with efficient inference yields better predictive log-loss and uncertainty calibration than either empirical-variance smoothed approaches or stochastic kriging.

Practical success relies on:

• Consistent and stable estimation of the log-variance process (whether by GP smoothing or kernel regression).
• Reliable estimation of the importance weights, with sufficient samples $M$ to stabilize the mixture moments.
• Careful kernel and hyperparameter initialization to avoid degenerate variance surfaces.
• Fast approximate inference for large datasets, preferably with distributed architectures or mini-batch variational surrogates.

Empirical studies show that HGPR systematically improves both negative log-predictive density (calibration) and mean squared error in presence of heteroscedasticity—especially in multi-replication, large-sample, or safety-critical contexts.

7. Key Insights, Open Problems, and Future Directions

• HGPR is the standard Bayesian approach for regression with input-dependent noise, with exact theory now available for all posterior moments and the full predictive CDF as mixtures of Gaussians (Ito, 19 Sep 2024).
• Posterior distributions are non-Gaussian, mixture-valued, and require either importance sampling, Laplace- or variational-approximate inference, or (for small data) explicit MCMC.
• Modern inference leverages autonormalized importance sampling, high-efficiency blockwise matrix computation (Woodbury), and sparse/distributed approximations for scalability (Liu et al., 2018, Lee et al., 2023, Binois et al., 2016).
• The non-asymptotic mixture structure enables principled enforcement of chance constraints in control and design under uncertainty.
• HGPR methodology is broadly extensible: Student-t or contaminated normal models yield robustness to outliers (Hartmann et al., 2017, Iong et al., 27 Feb 2024), multi-output and matrix-variance extensions cover high-dimensional and structured outputs (Lee et al., 2023), and compositional or latent-covariate approaches afford non-Gaussian residual modeling (Wang et al., 2012).
• Open directions include further scalability (e.g. GPU/TPU accelerations, adaptive expert allocation (Liu et al., 2018)), principled non-stationary priors for the signal and variance processes, and rigorous quantification of the impact of approximate sampling on downstream risk-sensitive decision making.

Heteroscedastic GP regression has thus evolved into a mature, well-theorized, and flexible modeling paradigm providing comprehensive uncertainty quantification under arbitrary and unknown input-dependent noise, with tractable and theoretically principled inference methods now available for a wide range of application domains (Ito, 19 Sep 2024, Hartmann et al., 2017, Liu et al., 2018, Lee et al., 2023, Binois et al., 2016, Wang et al., 2012).