Robust and Conjugate Gaussian Process (RCGP)

• RCGP is a generalized-Bayesian GP regression method that replaces the Gaussian log-likelihood with a weighted score-matching loss to achieve robust conjugacy.
• It preserves exact Gaussian posterior and predictive distributions by ensuring the loss remains quadratic in the latent function values.
• RCGP robustly handles outliers and model misspecification via adaptive noise inflation and target shrinkage, all without sacrificing computational efficiency.

Searching arXiv for the primary RCGP paper and closely related extensions/comparisons. Robust and Conjugate Gaussian Process (RCGP) regression is a generalized-Bayesian reformulation of Gaussian process regression in which robustness to contamination and model misspecification is obtained by replacing the Gaussian log-likelihood in Bayes’ rule with a weighted score-matching loss, while preserving exact Gaussian posterior and predictive calculations. In the original construction, the latent prior remains $f\sim \mathcal{GP}(m,k)$, the observation model remains Gaussian for the purpose of defining the score, and the resulting posterior stays in closed form because the robust loss is quadratic in the latent function values (Altamirano et al., 2023).

1. Definition and scope

Standard GP regression assumes

$y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$

with

$f \sim \mathcal{GP}(m,k).$

At observed inputs, the latent vector $\mathbf f$ has prior

$$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$

and the Gaussian likelihood gives exact conditioning: $$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$ with

$$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$

This conjugate structure is computationally convenient, but the Gaussian noise model is fragile under outliers, heavy tails, heterogeneity, and broader misspecification (Altamirano et al., 2023).

RCGP addresses that tension by changing the updating rule rather than replacing the GP prior or introducing a non-Gaussian likelihood. The method is therefore distinct from Student-$t$-likelihood, Laplace, mixture, or heteroskedastic robust GP constructions, which typically lose exact Gaussian conditioning and require variational methods, Laplace approximations, expectation propagation, or MCMC (Altamirano et al., 2023). It is also distinct from earlier conjugate GP work whose emphasis was scalability rather than outlier robustness; for example, large-spatial-data models based on exact conjugate updating and sparse-plus-low-rank structure are conjugate but not robust in the outlier-resistant sense of RCGP (Shirota et al., 2019).

A related but separate line of work constructs auxiliary-variable representations that make non-conjugate GP models conditionally conjugate. That strategy is broad and includes robust likelihoods such as Student-$t$ and Laplace, but the conjugacy is conditional on augmentation variables rather than direct as in RCGP regression (Galy-Fajou et al., 2020).

2. Generalized-Bayes construction

The original RCGP starts from generalized Bayes: $$p_{\beta}^{L}( \mathbf{f}\mid \mathbf{x},\mathbf{y})\propto p(\mathbf{f}\mid\mathbf{x}) \exp\big(-\beta n L_n(\mathbf{f},\mathbf{y},\mathbf{x})\big).$$ Instead of the negative log-likelihood, it uses a loss derived from a weighted Fisher divergence. In conditional form, the underlying divergence is

$y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$0

with model score

$y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$1

under the Gaussian observation model. The weighted version is

$y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$2

After integration by parts, this yields the empirical loss

$y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$3

and the RCGP posterior

$y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$4

The learning-rate parameter $y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$5 is absorbed into the weight function $y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$6 (Altamirano et al., 2023).

The crucial algebraic fact is that, under Gaussian observation modeling, $y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$7 is affine in $y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$8. Hence $y_i = f(x_i) + \epsilon_i,\qquad \epsilon_i \overset{\text{i.i.d.}{\sim}} \mathcal N(0,\sigma^2),$9 is quadratic in $f \sim \mathcal{GP}(m,k).$0: $f \sim \mathcal{GP}(m,k).$1 Multiplying a Gaussian prior by the exponential of a negative quadratic leaves a Gaussian posterior. This is the exact source of conjugacy in RCGP (Altamirano et al., 2023).

3. Closed-form posterior and predictive distributions

The main closed-form result introduces a weight vector $f \sim \mathcal{GP}(m,k).$2, a weight-dependent diagonal correction $f \sim \mathcal{GP}(m,k).$3, and a shifted mean term

$f \sim \mathcal{GP}(m,k).$4

With these terms, the RCGP posterior is

$f \sim \mathcal{GP}(m,k).$5

where

$f \sim \mathcal{GP}(m,k).$6

$f \sim \mathcal{GP}(m,k).$7

At a new input $f \sim \mathcal{GP}(m,k).$8,

$f \sim \mathcal{GP}(m,k).$9

with

$\mathbf f$0

$\mathbf f$1

These expressions are formally parallel to standard GP regression, except that the identity-noise term is replaced by a weight-dependent diagonal term and the centered target vector is replaced by $\mathbf f$2 (Altamirano et al., 2023).

The preferred weight in the original paper is the inverse multiquadric (IMQ)

$\mathbf f$3

with $\mathbf f$4. The parameter $\mathbf f$5 is a soft outlier threshold, while $\mathbf f$6 is the generalized-Bayes learning-rate scale. As $\mathbf f$7, $\mathbf f$8, recovering standard GP regression (Altamirano et al., 2023).

This representation admits a direct interpretation. The weight-dependent diagonal term performs adaptive noise inflation, so observations with small weight are treated as effectively noisier. The shifted target term introduces adaptive shrinkage through $\mathbf f$9. The original paper notes that the dominant robustness effect comes from the noise term (Altamirano et al., 2023).

4. Robustness mechanism, guarantees, and misconceptions

RCGP’s robustness is formalized through the posterior influence function (PIF). If one observation $$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$0 is replaced by a contaminated value $$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$1, and $$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$2 is the contaminated dataset, the PIF is

$$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$3

A posterior is called robust if this remains bounded as $$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$4 (Altamirano et al., 2023).

For standard GP regression,

$$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$5

so the effect of arbitrarily large contamination diverges. For RCGP, assuming $$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$6,

$$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$7

Hence, if

$$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$8

the posterior influence is bounded, and RCGP is robust in this KL-sense (Altamirano et al., 2023).

Several misconceptions follow from this formulation. First, RCGP is not a heavy-tailed GP likelihood model; the observation model used to define the score is still Gaussian, and robustness comes from generalized-Bayes reweighting rather than from a Student-$$p(\mathbf{f}\mid \mathbf{x})=\mathcal{N}(\mathbf{f};\mathbf{m},K),$$9, Laplace, or contamination likelihood (Altamirano et al., 2023). Second, robustness in the formal PIF sense requires $$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$0 to depend on $$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$1; heteroskedastic GP regression can be recovered as a special case of the algebra, but it is not robust in this sense because its weight depends only on $$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$2 (Altamirano et al., 2023). Third, the method is sensitive to the centering implicit in $$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$3. The original paper explicitly notes that a poor prior mean $$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$4 can cause legitimate observations to be treated as outliers or outliers to receive too much weight (Altamirano et al., 2023).

This mean-dependence became a central criticism in later work. A comparison paper on robust Gaussian processes via relevance pursuit describes RCGP as a preweighted robustification whose success depends on the a priori weighting rule identifying likely outliers from their proximity to the prior mean or target median, and shows that this can fail when corruptions are not clearly separated in output space (Ament et al., 2024).

5. Hyperparameter learning, computation, and technical extensions

The original paper derives a closed-form pseudo marginal likelihood for RCGP, but explicitly states that maximizing it is ill-posed because the exponential loss term is not a normalized density. Instead, hyperparameters are learned by optimizing leave-one-out predictive log density,

$$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$5

with an analytic $$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$6 expression rather than a naive $$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$7 implementation. The computational complexity of exact RCGP therefore remains the same order as exact GP regression,

$$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$8

and the paper characterizes the method as operating at virtually no additional cost relative to standard GP regression (Altamirano et al., 2023).

The original construction also extends to sparse variational inference. The robust conjugate sparse variational GP (RCSVGP) preserves Gaussian algebra under inducing-point approximations, showing that the generalized-Bayes quadratic structure is compatible with standard sparse GP machinery (Altamirano et al., 2023). Later work broadened this trajectory in several directions.

For spatio-temporal problems, ST-RCGP embeds the robust update in a state-space formulation and derives Kalman-like generalized-Bayes filtering and smoothing recursions. This yields linear-in-time cost $$p(\mathbf{f}\mid \mathbf{x},\mathbf{y}) = \mathcal{N}(\mathbf{f}; \mu,\Sigma),$$9, replaces the fragile fixed center $$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$0 by the filtering predictive mean $$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$1, sets the shrinkage scale by predictive variance $$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$2, and introduces a robust hyperparameter objective $$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$3 (Laplante et al., 4 Feb 2025).

For correlated outputs, MO-RCGP extends the construction to vector-valued GPs. It replaces the scalar weight by a diagonal matrix-valued $$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$4, preserves a Gaussian posterior over the vectorized latent field, and proves bounded posterior influence not only in the contaminated channel but across all outputs. A notable modeling change is that the centering function for output $$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$5 is taken to be the conditional mean of $$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$6 given the other outputs, allowing cross-output information to support outlier detection (Rooijakkers et al., 30 Oct 2025).

For approximate large-scale computation, RCaGP combines the robust generalized-Bayes core of RCGP with computation-aware uncertainty over approximate representer weights. Its main novelty is not conjugacy itself but the joint treatment of outliers and approximation-induced uncertainty, with an added computational-uncertainty term in the predictive variance (Sinaga et al., 27 May 2025).

In Bayesian optimization, later work uses RCGP as the surrogate inside RCGP-UCB for unbounded-magnitude corruptions. The key modification is a plateau-IMQ weight that behaves exactly like standard GP regression on non-outlying observations and smoothly downweights only those outside a trusted plateau, yielding “zero-cost robustness” and regret bounds that remain sublinear under corruption-frequency constraints (Ezzerg et al., 19 Nov 2025).

6. Empirical behavior, applications, and neighboring formulations

The original empirical study evaluates RCGP on synthetic regression and on Boston, Energy, and Yacht datasets under no outliers, focused outliers, asymmetric outliers, and uniform outliers. The general pattern is that RCGP is essentially as good as standard GP when no outliers are present, while improving substantially under contamination. Reported examples include Boston training times of GP $$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$7 sec, RCGP $$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$8 sec, Student-$$\mu = \mathbf{m} +K (K+\sigma^{2} I_{n})^{-1}(\mathbf{y}-\mathbf{m}), \qquad \Sigma = K (K+ \sigma^{2} I_{n})^{-1} \sigma^{2} I_{n}.$$9 GP $t$0 sec, and mixture GP $t$1 sec, as well as Energy training times of GP $t$2, RCGP $t$3, Student-$t$4 GP $t$5, and mixture GP $t$6 (Altamirano et al., 2023).

The paper also reports a contaminated Bayesian optimization study on six-hump camel, Branin, McCormick, and Rosenbrock functions, where RCGP improves cumulative regret relative to standard GP and remains far cheaper than Student-$t$7-likelihood baselines (Altamirano et al., 2023). A qualitative finance example based on the Twitter Flash Crash further shows that ordinary GP overfits the transient crash, whereas RCGP downweights the anomalous observations and yields a smoother fit (Altamirano et al., 2023). ST-RCGP later revisits this failure mode and argues that original RCGP can still fail when the prior-mean-based center lies close to the outliers, motivating adaptive centering from the filtering predictive (Laplante et al., 4 Feb 2025).

Subsequent comparisons clarify where RCGP sits within the robust GP literature. Relevance pursuit (RRP) keeps exact Gaussian tractability conditional on per-point variance inflation, but introduces robustness by adaptively selecting a sparse subset of suspect observations and inflating only their noise variances. That paper positions RRP against RCGP and argues that RCGP is a preweighted conjugate Gaussian robustification, whereas RRP is an adaptive sparse variance-inflation GP learned by marginal-likelihood pursuit (Ament et al., 2024). This suggests a useful conceptual distinction: RCGP is a direct generalized-Bayes robustification of conjugate GP regression, while several later methods are exact-or-conditionally-exact Gaussian alternatives that target different robustness regimes.

In the broader conjugate-GP landscape, it is also important not to conflate RCGP with methods that are merely conjugate. Large-spatial-data models based on exact conjugate Bayesian updating with sparse-plus-low-rank structure are highly relevant to the conjugate side of the story, but they are not robust GP methods in the formal outlier-resistant sense (Shirota et al., 2019). Conversely, augmentation-based schemes for non-conjugate GP models show that robust likelihoods such as Student-$t$8 and Laplace can be made conditionally conjugate, but those constructions solve a different problem from the direct robust-and-conjugate regression formulation of RCGP (Galy-Fajou et al., 2020).

RCGP therefore occupies a specific position: it is a GP regression framework in which robustness is achieved through generalized-Bayes weighting, conjugacy is retained because the resulting loss remains quadratic in the latent field, and later work has specialized or extended this structure to sparse variational inference, spatio-temporal state-space models, multi-output regression, computation-aware approximation, and robust Bayesian optimization (Altamirano et al., 2023).