Variational Outlier-Robust Gaussian Process Regression with Generative Modeling

IEEE SIGNAL PROCESSING LETTERS, VOL. XX, NO. XX, 2026

1

Variational Outlier-Robust Gaussian Process Regression with Generative Modeling

arXiv:2608.16606v1 [cs.LG] 17 Aug 2026

Arslan Majal and Aamir Hussain Chughtai

Abstract—Outliers can substantially distort Gaussian process regression (GPR) due to its conventional Gaussian observation likelihood, leading to inaccurate model learning and prediction. To address this limitation, this article introduces a generative GPR model that captures observation-specific contamination and adaptively mitigates the influence of outliers. Subsequently, a variational generalized expectation-maximization procedure is used to learn the latent variables and GPR model parameters. Experiments on synthetic and real datasets under different contamination settings demonstrate that the proposed method remains competitive with—and in several cases outperforms— robust GPR baselines in prediction accuracy. Moreover, the proposed method shares the cubic computational scaling of the compared GPR methods. Index Terms—Gaussian process regression, outlier-robust regression, variational Bayesian inference, robust signal processing, uncertainty quantification, expectation-maximization

I. I NTRODUCTION AUSSIAN process regression (GPR) [1] is a Bayesian nonparametric regression method that provides flexible nonlinear modeling and principled uncertainty quantification, while remaining effective in data-limited settings. These properties make GPR useful in numerous signal-processing and learning applications. Standard GPR assumes Gaussian observation noise, making it sensitive to outliers caused by sensor faults, calibration errors, transmission noise, or unexpected disturbances. Under this model, large residuals are penalized quadratically, allowing even a few corrupted observations to substantially distort the estimated latent function [2]. This motivates robust GPR methods that reduce the influence of unreliable observations. Outliers in GPR are commonly handled by modifying the observation likelihood or the corresponding loss function. Existing robust approaches often rely on fixed likelihood modifications or robust losses with user-defined tuning parameters, including Huber- and Hampel-type penalties, heavytailed Student-t likelihoods, and generalized-Bayes likelihoodinduced losses [3]–[5]. Although tractable, these approaches can be sensitive to prescribed thresholds, loss parameters, or distributional shape parameters that control the influence of atypical observations. To overcome these limitations, more expressive adaptive models have been proposed. Instead of fixing the robust cost or likelihood parameters a priori, these methods adapt them from the observed data, thereby providing data-driven robustness.

G

A. Majal is with Zunitech, LLC, Houston, USA. A. H. Chughtai is with the Institute of Data Science, University of Engineering and Technology, Lahore, Pakistan.

The main challenge, however, is to make the outlier model sufficiently expressive while retaining tractable inference. To this end, existing approaches introduce latent variables, hierarchical likelihoods, mixture-based noise models, or optimizationbased robust formulations, including graduated non-convexity (GNC), to account for corrupted measurements [6]–[13]. A relevant class of adaptive methods employs generative models with variational inference to identify and mitigate outliers [7], [11]. Building on this idea, adaptive selective observation rejecting (ASOR) [10] learns outlier statistics rather than completely discarding suspected observations, thereby adapting their influence according to the learned contamination characteristics. ASOR has demonstrated favorable estimation accuracy and computational efficiency relative to state-of-theart GNC methods in robotic perception problems [13]. Its conjugate hierarchical structure enables closed-form variational updates for the model parameters [10], with related ideas extended to robust filtering in [12]. Motivated by the adaptive outlier treatment in ASOR, we propose a formulation tailored to outlier-robust GPR. The key development is to integrate observation-specific precision modeling with the latent GP function and its associated noise, mean, and kernel parameters within a unified model. Since exact Bayesian inference is intractable, we employ a variational generalized expectation-maximization (EM) procedure to learn the latent variables and model parameters. Closed-form updates are obtained for all quantities except the kernel hyperparameters, which are optimized using the gradient-based method [9]. To evaluate performance, the proposed method is extensively compared with representative GPR baselines on synthetic and real-data regression problems under diverse contamination settings.

II. P ROPOSED M ETHODOLOGY Notation: Nonbold symbols denote scalars or scalar-valued functions, bold lowercase symbols denote vectors, and bold uppercase symbols denote matrices. Moreover, (·)⊤ denotes transpose, [A]ij denotes the (i, j)th entry of A, and diag(·) and tr(·) denote the diagonal and trace operators, respectively. Consider n input–output samples arranged as   X =   i ⊤ (x )    

x11 ..

···

x1dx ..

xi1 ..

···

xidx ..

xn1

· · · xndx





   Y =      ,         

y1j ..

···

yij ..

···

y11 ..

···

yi1 ..

···

yn1

· · · ynj yj

y1dy .. yidy ..

· · · yndy

      i ⊤ (y )    

IEEE SIGNAL PROCESSING LETTERS, VOL. XX, NO. XX, 2026

2

where X ∈ Rn×dx contains the input samples and Y ∈ Rn×dy where A > 1 and B > 0 are fixed hyperparameters. Similarly, the nominal observation-noise variance is assigned contains the corresponding outputs. The ith rows are (xi )⊤ and the inverse-Gamma prior (y i )⊤ , respectively, while yj denotes the jth output column.   Assuming independent output components and contamination 2 2 2 ν0 s0j p(σ ) = IG σ ∝ (σj2 )−(ν0 +2)/2 e−s0j /2σj . (7) , j j processes, the problem is decomposed into dy scalar-output 2 2 GPR models, where the jth model is trained using X and yj . where ν > 0 and s > 0 are fixed hyperparameters. 0 0j The fixed prior hyperparameters associated with output A. Scalar Observation Model dimension j are collected as For output dimension j, let fij ≜ fj (xi ) denote the latent Pj ≜ {a, {θij }ni=1 , A, B, ν0 , s0j } . (8) function value at the ith input. The observation model is yij = fij + εij ,

i = 1, . . . , n,

C. Variational Generalized-EM Inference

where εij denotes the observation error and σj2 denotes its nominal variance. Following the ASOR construction, a positive latent precision multiplier Iij > 0 is introduced such that
p(εij | Iij , σj2 ) = N εij | 0, σj2 /Iij ,

For output dimension j, define I j ≜ [I1j , . . . , Inj ]⊤ , Zj ≜ {fj , I j }, Θj ≜ {mj , κj , σj2 , bj }. (9) For j = 1, . . . , dy , the joint posterior satisfies 2

p(Zj , Θj | yj , X, Pj ) ∝ p(fj | X, mj , κj )p(σj )p(bj ) where N (x | µ, v) denotes a Gaussian density in x with mean n Y µ and variance v. Accordingly, the observation likelihood is " # × p(Iij | bj )p(yij | fij , Iij , σj2 ). I ij 1/2 i=1 p(yij | fij , Iij , σj2 ) = (2πσj2 )−1/2 Iij exp − 2 (yij − fij )2 . (10) 2σj We use flat non-informative priors for m and κ , whose j j (1) The value Iij = 1 recovers the nominal observation model, constant terms are omitted from (10). Since the exact joint whereas Iij ̸= 1 adaptively modifies the observation variance posterior does not admit closed form, we employ variational generalized EM [14]. to account for outliers. The E-step uses the mean-field factorization Let fj ≜ [f1j , . . . , fnj ]⊤ . The GP prior is n Y p(fj | X, mj , κj ) = N (fj | mj 1n , Kj ) , q (Z ) = q (f ) qj (Iij ). (11) j j j j (2) i=1 [Kj ]iℓ = kj (xi , xℓ ; κj ). For ξ ∈ {fj , I1j , ..., Inj }, the coordinate-optimal update is where 1n ∈ Rn denotes the vector of ones and, for the multivariate case, N (x | µ, Σ) denotes a Gaussian density log qj⋆ (ξ) = Eqj (Zj \ξ) [log p(yj , Zj | X, Θj , Pj )] + c, (12) with mean µ and covariance matrix Σ. The covariance matrix is constructed using the automatic relevance determination where c is independent of ξ. Applying (12) to fj and each Iij gives the closed-form updates below. At each generalized(ARD) squared-exponential kernel Pdx EM iteration, these factors are updated sequentially, and their 2 2 2 kj (xi , xℓ ; κj ) = σf,j e−0.5 r=1 (xir −xℓr ) /ℓjr , (3) product gives qj (Zj ). 2 The M-step aims to solve κj ≜ [σf,j , ℓj1 , . . . , ℓjdx ]⊤ . 2 Here, σf,j is the signal variance and ℓjr is the length scale associated with the rth input dimension.

B. Hierarchical Priors Iij is assigned the unit-spike-and-Gamma prior p(Iij | bj ) = (1 − θij )G(Iij | a, bj ) + θij δ(Iij − 1),

(4)

where δ(·) denotes the Dirac delta function and G(I | a, b) denotes the Gamma density with shape a > 0 and rate b > 0: G(I | a, b) ≜

ba a−1 −bI I e . Γ(a)

(13)

Θj

where Qj (Θj ) ≜ Eqj (Zj ) [log p(yj , Zj , Θj | X, Pj )]. Closedform updates are obtained for mj , σj2 , and bj , whereas κj is updated using gradient descent with backtracking. The priors on σj2 and bj yield maximum a posteriori (MAP) updates, whereas the uninformative priors on mj and κj yield maximumlikelihood updates. D. Variational E-Step: Latent-Function Update

(5)

where Γ(·) denotes the Gamma function. The parameter θij ∈ (0, 1) is the prior probability that yij belongs to the nominal branch. To retain Bayesian conjugacy, the output-specific Gamma rate bj is assigned the Gamma hyperprior p(bj ) = G(bj | A, B),

Θ⋆j ∈ arg max Qj (Θj ).

(6)

For fixed output j, define wij ≜ Eqj (Iij ) [Iij ], Λj ≜ diag(w1j , . . . , wnj )/σj2 .

(14)

Applying (12) to fj gives qj (fj ) = N (fj | µf,j , Cf,j ) ,
−1 Cf,j = Kj−1 + Λj , µf,j = Cf,j

Kj−1 mj 1n + Λj yj

(15)

.

IEEE SIGNAL PROCESSING LETTERS, VOL. XX, NO. XX, 2026

3

E. Variational E-Step: Precision-Multiplier Update

Algorithm 1: ASOR-GPR d

y Input: Training data (X, Y ); prior hyperparameters {Pj }j=1 ;

Using the moments of qj (fj ), define Sij ≜ (yij − µf,ij )2 + [Cf,j ]ii , Rij ≜ Sij /σj2 ,

(16)

where µf,ij is the ith entry of µf,j . Let e a ≜ a + 1/2, ebij ≜ bj + Rij /2.

(17)

Applying (12) to Iij gives   qj (Iij ) = (1 − Ωij )G Iij | e a, ebij + Ωij δ(Iij − 1),

(18)

where  #−1 1 − θij Γ(e a) baj Rij Ωij = 1 + . exp a θij Γ(a) ebeij 2 "

(19)

wij = Ωij + (1 − Ωij )e a/ebij .

(20)

F. Generalized M-Step: Parameter Updates With the variational quantities fixed, Qj (Θj ) separates into terms involving σj2 , bj , and (mj , κj ). Collecting the terms that depend on σj2 and setting their derivative to zero gives the MAP update Sj + s0j . n + ν0 + 2

(21)

Pn where Sj ≜ i=1 wij Sij . Similarly, the MAP update of the Gamma rate parameter is Pn A − 1 + a i=1 (1 − Ωij ) bj ← . (22) Pn B + i=1 (1 − Ωij )e a/ebij For fixed κj , maximizing the GP-dependent terms of Qj (Θj ) with respect to mj gives −1 −1 ⊤ mj ← (1⊤ n Kj µf,j )/(1n Kj 1n ).

(23)

For the kernel update, define  ⊤ ϕj ≜ log σf,j log ℓj1 · · · log ℓjdx ,  ⊤ κj (ϕj ) ≜ exp(2ϕj1 ) exp(ϕj2 ) · · · exp(ϕj,dx +1 ) . (24) This logarithmic parameterization enforces positivity of the signal variance and length scales while allowing optimization over ϕj . For fixed mj , maximizing the terms of Qj (Θj ) that depend on κj through Kj is equivalent to minimizing Jj (ϕj ) ≜

Output: Converged independent-output ASOR-GPR models. for j ← 1 to dy do Extract yj , initialize Θj , set ϕj from κj using (24), set wij ← 1, i = 1, . . . , n, and construct Kj repeat Variational E-step: Update qj (fj ) using (14) and (15) Update {qj (Iij ), Ωij , wij }n i=1 using (16)–(20) Generalized M-step: Update σj2 , bj , and mj using (21)–(23) Update ϕj by gradient descent with backtracking using (25)– (30) Update κj using (24) and reconstruct Kj using (3) until convergence Perform a final variational E-step and store the converged quantities Return The converged ASOR-GPR models

The corresponding robustness weight is

σj2 ←

dy initial parameters {Θj }j=1 .

1 1 tr Kj−1 Dj + log |Kj | . 2 2

(25)

where ⊤

Dj ≜ Cf,j + (µf,j − mj 1n ) (µf,j − mj 1n ) .

The required kernel derivatives are   ∂Kj ∂Kj (xir − xℓr )2 = 2Kj , = [Kj ]iℓ . (28) ∂ϕj1 ∂ϕj,r+1 iℓ ℓ2jr For a trial step size αj > 0, the normalized-gradient candidate is gj ϕcand = ϕold , (29) j j − αj max{1, ∥gj ∥2 } where gj ≜ ∇ϕj Jj . A backtracking line search halves αj until the Armijo condition Jj (ϕcand ) ≤ Jj (ϕold j j ) − c1 αj

To obtain a descent direction for minimizing Jj (ϕj ), its gradient components are  
∂Kj ∂Jj 1 = tr Kj−1 − Kj−1 Dj Kj−1 . (27) ∂ϕjr 2 ∂ϕjr

(30)

is satisfied, where c1 is the Armijo constant. The accepted value is ϕj ← ϕcand , followed by κj ← κj (ϕj ) and Kj ← j Kj (κj ). Iterations terminate when the relative change in the weighted residual diagnostic Sj falls below a prescribed tolerance or the maximum iteration count is reached. G. Posterior Prediction b j , and σ Let m b j, κ bj2 denote the learned parameters, and let   bf,j . b f,j , C qj (fj ) = N fj | µ (31) For a test input x∗ , let f∗j ≜ fj (x∗ ) and define 
n b∗j ≜ kj xi , x∗ ; κ cj ≜ Kj (b b j i=1 , K κj ) , k b bj ) . k∗∗j ≜ kj (x∗ , x∗ ; κ

(32)

Marginalizing the conditional GP distribution over (31) gives qj (f∗j ) = N (f∗j | µ∗j , v∗j ) , b⊤ K c−1 (µ b f,j − m µ∗j = m bj + k b j 1n ) , ∗j

(26)

∥gj ∥22 max{1, ∥gj ∥2 }

j

b∗j + k b⊤ K b∗j . b⊤ K bf,j K c−1 k c−1 k c−1 C v∗j = b k∗∗j − k ∗j ∗j j j j (33) For a nominal future observation y∗j , the predictive distribution is obtained as
qj (y∗j ) = N y∗j | µ∗j , v∗j + σ bj2 . (34)

IEEE SIGNAL PROCESSING LETTERS, VOL. XX, NO. XX, 2026

4

Oracle ASOR-GPR GMM-GPR Student-t GPR RCGPR Standard GPR

5

4

RMSE

RMSE

3

3.5

Oracle ASOR-GPR GMM-GPR Student-t GPR RCGPR Standard GPR

3 2.5

RMSE

5

4

3

2

Oracle ASOR-GPR GMM-GPR Student-t GPR RCGPR Standard GPR

2 1.5

2 1 1

1

0.5

0

0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0 0.0

0.1

0.2

0.3

(a) Synthetic data: uniform outliers

0.4

0.5

0.6

0.7

pout

pout

(b) Air Quality: uniform outliers

0.8

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

pout

(c) Energy Efficiency: Gaussian outliers

Fig. 1. Prediction performance under different outlier settings: (a) synthetic data with uniform outliers over [0, 10]; (b) Air Quality data with uniform outliers over [−5, 15]; and (c) Energy Efficiency data with zero-mean Gaussian outliers of variance 25.

Fig. 1a considers outliers drawn uniformly from [0, 10]. ASOR-GPR remains competitive across the contamination All experiments were conducted in MATLAB R2026a on an range and achieves the lowest RMSE among the compared AMD Ryzen 9 9950X desktop with 64 GB RAM. We compare non-oracle methods at most contamination levels. For real-data ASOR-GPR with standard GPR [1], robust conjugate GPR evaluation, we use the Air Quality [15] and Energy Efficiency (RCGPR) [5], Student-t GPR [4], Gaussian-mixture-likelihood [16] datasets, standardized using training-set statistics. Using GPR (GMM-GPR) [9], and an Oracle GPR supplied with all the same probability sweep, Fig. 1b considers uniform outliers available ground-truth contamination information. over [−5, 15], while Fig. 1c considers zero-mean Gaussian All methods use independent scalar-output models with the outliers with variance 25. ASOR-GPR provides clear gains (0) ARD squared-exponential kernel, initialized with σf,j = 1 under asymmetric contamination in Fig. 1b and remains com(0) and ℓjr = 0.2dmed , where dmed is the median pairwise petitive with robust likelihood-based methods under symmetric training-input distance. All non-oracle methods are initialized Gaussian outliers in Fig. 1c. Similar trends are observed over with nominal-noise variance 0.5. For synthetic data, this is a wider range of contamination intensity. This reflects the deliberately misspecified relative to the true variance 0.25. For effectiveness of the proposed adaptive inference of observationthe real datasets, the true nominal-noise variance is unknown. specific precisions and nominal/outlier associations. Predictive experiments use ntrain = ntest = 100. I Outliers are injected entry-wise only into the training targets M EDIAN SERIAL MODEL - FITTINGTABLE TIME ( S ) OVER 30 M ONTE C ARLO TRIALS . over 30 Monte Carlo trials. Root mean square error (RMSE) Method ntrain for synthetic data is evaluated against the known latent function 200 400 600 800 1000 at the test inputs. For real data, the original observations are Oracle GPR 0.10 0.70 1.99 4.37 6.76 Std. GPR 0.19 1.46 4.24 9.43 15.05 treated as nominal, with outliers injected only into the training RCGPR 0.50 4.74 15.58 36.45 62.73 Student-t GPR 3.25 16.23 41.39 80.75 123.53 targets and RMSE evaluated against the unmodified test targets, ASOR-GPR 24.60 31.13 85.52 212.41 335.96 which retain their inherent observation noise. GMM-GPR 31.96 151.23 376.66 644.76 963.38 For runtime evaluation, Table I reports median model-fitting Baseline-specific parameters follow their respective reftimes with increasing training-set size at pout = 0.4 under zeroerences unless otherwise stated. ASOR-GPR uses a = 1, (0) mean Gaussian contamination with variance 10. Oracle GPR θij = 0.5, A = 10, B = 0.45, ν0 = 3, bj = 20, and standard GPR have the lowest runtimes, with the Oracle and s0j = 1.25 in all real-data and synthetic experiments. benefiting from known contamination information. RCGPR and Up to 1000 variational iterations are run, terminating when −5 Student-t GPR incur moderate computational cost. ASOR-GPR relative change in Sj falls below 10 . Kernel hyperparameters and GMM-GPR exhibit the highest runtimes, with ASORare optimized by normalized gradient descent with Armijo GPR consistently faster than GMM-GPR over the reported backtracking using initial step αj = 0.05, Armijo constant −4 −7 training sizes. This relative computational behavior remains c1 = 10 , minimum step αmin = 10 , and at most 25 inner broadly consistent across different contamination probabilities. iterations. All methods involve dense covariance-matrix operations and For the synthetic experiment, dx = 2 and dy = 4. Inputs 2 retain cubic complexity, while runtime differences arise from are Latin-hypercube sampled uniformly over [−3, 3] . Defining zir = (xir + 3)/6, r = 1, 2, the clean outputs are generated as their method-specific inference procedures. III. N UMERICAL E XPERIMENTS

gj (xi ) = 10 sin(πzi1 zi2 ) + 5 + 0.75 sin((0.8 + 0.1j)xi1 ) + 0.50 cos((0.6 + 0.05j)xi2 ) + 0.10j xi1 ,  ⊤ i ge(x ) = Cy⊤ g1 (xi ) · · · gdy (xi ) ,   Cy = Idy + 0.1 1dy 1⊤ dy − I dy . (35) Here, ge(xi ) denotes the mixed clean output vector. The outputs are standardized to zero mean and unit variance before adding nominal Gaussian noise and entry-wise contamination, with pout ∈ {0, 0.1, . . . , 0.8}.

IV. C ONCLUSION We proposed an outlier-robust GPR approach that embeds adaptive observation-specific precision modeling within a hierarchical generative framework. Variational generalized-EM inference enables observation-specific contamination effects and the underlying GP model to be learned jointly. Results on synthetic and real datasets show that ASOR-GPR achieves competitive predictive performance with enhanced robustness across a range of contamination settings, while preserving the cubic computational scaling of the standard GPR.

IEEE SIGNAL PROCESSING LETTERS, VOL. XX, NO. XX, 2026

R EFERENCES [1] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning. MIT Press, 2006. [2] P. J. Huber, “Robust estimation of a location parameter,” The Annals of Mathematical Statistics, vol. 35, no. 1, pp. 73–101, 1964. [3] L. Chang, K. Li, and B. Hu, “Huber’s m-estimation-based process uncertainty robust filter for integrated INS/GPS,” IEEE Sensors Journal, vol. 15, no. 6, pp. 3367–3374, 2015. [4] P. Jylänki, J. Vanhatalo, and A. Vehtari, “Robust Gaussian process regression with a Student-t likelihood,” JMLR, 2011. [5] M. Altamirano, F.-X. Briol, and J. Knoblauch, “Robust and conjugate Gaussian process regression,” ICML, 2024. [6] H. Wang, H. Li, J. Fang, and H. Wang, “Robust Gaussian Kalman filter with outlier detection,” IEEE Signal Processing Letters, vol. 25, no. 8, pp. 1236–1240, 2018. [7] A. H. Chughtai, M. Tahir, and M. Uppal, “Outlier-robust filtering for nonlinear systems with selective observations rejection,” IEEE Sensors Journal, vol. 22, no. 7, pp. 6887–6897, 2022. [8] A. Nakabayashi and G. Ueno, “Nonlinear filtering method using a switching error model for outlier-contaminated observations,” IEEE Transactions on Automatic Control, vol. 65, no. 7, pp. 3150–3156, 2019. [9] A. Daemi, H. Kodamana, and B. Huang, “Gaussian process modelling with Gaussian mixture likelihood,” Journal of Process Control, vol. 81, pp. 209–220, 2019. [10] A. H. Chughtai, M. Tahir, and M. Uppal, “Bayesian heuristics for robust spatial perception,” IEEE Transactions on Instrumentation and Measurement, vol. 73, pp. 1–11, 2024. [11] ——, “EMORF/S: EM-based outlier-robust filtering and smoothing with correlated measurement noise,” IEEE Transactions on Signal Processing, vol. 72, pp. 4318–4331, 2024. [12] A. Majal, A. H. Chughtai, and M. Tahir, “EMORF-II: Adaptive EM-based outlier-robust filtering with correlated measurement noise,” in 2025 IEEE 35th International Workshop on Machine Learning for Signal Processing (MLSP), 2025, pp. 1–6. [13] H. Yang, P. Antonante, V. Tzoumas, and L. Carlone, “Graduated nonconvexity for robust spatial perception: From non-minimal solvers to global outlier rejection,” IEEE Robotics and Automation Letters (RA-L), vol. 5, no. 2, pp. 1127–1134, 2020. [14] V. Šmı́dl and A. Quinn, The Variational Bayes Method in Signal Processing. Springer, 2006. [15] S. Vito, “Air Quality,” UCI Machine Learning Repository, 2008, DOI: https://doi.org/10.24432/C59K5F. [16] A. Tsanas and A. Xifara, “Energy Efficiency,” UCI Machine Learning Repository, 2012, DOI: https://doi.org/10.24432/C51307.

5

Related documents

• Variational Outlier-Robust Gaussian Process Regression with Generative Modeling#625181
• Variational Outlier-Robust Gaussian Process Regression with Generative Modeling#622908
• An Analysis of Posterior Collapse, Parameterization and Initialization in Variational Deep Gaussian Processes#307024
• Manufacturing Outliers#682950
• A Generalized Tangent Approximation based Variational Inference Framework for Strongly Super-Gaussian Likelihoods#1000498