Orthogonal GP Regression

• Orthogonal GP regression is a set of methods that impose constraints between model components to prevent one part from absorbing structure meant for another.
• These techniques improve interpretability and identifiability by distinctly separating the parametric mean, stochastic residuals, and inducing subspaces.
• They offer computational benefits and robust extrapolation through optimized sparse approximations and decoupled variational frameworks.

Searching arXiv for recent and foundational papers on orthogonal Gaussian process regression and related orthogonal GP constructions. Orthogonal Gaussian Process Regression denotes a class of Gaussian-process constructions in which orthogonality constraints are imposed between distinct model components. Depending on the formulation, the constrained objects may be the stochastic residual and a parametric mean function, a GP bias term and the sensitivity directions of a mechanistic model, a principal inducing-point subspace and its residual process, or additive functional components. The common purpose is to prevent one component from absorbing structure that is meant to be attributed to another, thereby improving identifiability, interpretability, extrapolative behavior, or computational structure (Plumlee et al., 2016, Kuppa et al., 20 Feb 2026, Shi et al., 2019, Lu et al., 2022).

Construction	Orthogonality target	Representative paper
Orthogonalized kernel	GP residual vs. mean basis $g(x)$	(Plumlee et al., 2016)
Embedded model-error OGP	GP bias vs. $\nabla_\theta m(x;\theta^*)$	(Kuppa et al., 20 Feb 2026)
Sparse orthogonal VI	Principal inducing subspace vs. residual GP	(Shi et al., 2019)
Orthogonally decoupled VI	Mean residual basis vs. covariance basis	(Salimbeni et al., 2018)
Orthogonal additive kernel	Additive components $f_u$ under $p_u$	(Lu et al., 2022)
Variational orthogonal features	Spectral inducing features $\psi_m$ in $L^2$	(Burt et al., 2020)

1. Foundational formulation and the identifiability problem

The classical starting point is GP regression with a parametric mean,

$$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$

When $g(x)$ contains low-order polynomials or interactions, those directions may lie within the principal eigenspace of $k$. The stochastic part can then “soak up” variation that would otherwise be assigned to $m(\cdot)$. The resulting pointwise prediction can remain good, but $\nabla_\theta m(x;\theta^*)$0 may be poorly determined, so the estimated mean loses interpretability (Plumlee et al., 2016).

The orthogonal GP construction enforces that the stochastic component be orthogonal to the span of the mean basis over the entire domain $\nabla_\theta m(x;\theta^*)$1. In integral form, the requirement is

$\nabla_\theta m(x;\theta^*)$2

which implies

$\nabla_\theta m(x;\theta^*)$3

Defining

$\nabla_\theta m(x;\theta^*)$4

the orthogonalized kernel is

$\nabla_\theta m(x;\theta^*)$5

This is a rank-$\nabla_\theta m(x;\theta^*)$6 correction to the original kernel and can be viewed as a Schur complement. Replacing $\nabla_\theta m(x;\theta^*)$7 by $\nabla_\theta m(x;\theta^*)$8 yields generalized least squares estimation for $\nabla_\theta m(x;\theta^*)$9,

$f_u$0

with predictive mean and variance formally identical to standard GP formulas after substituting $f_u$1 and $f_u$2 for $f_u$3 and $f_u$4 (Plumlee et al., 2016).

A recurrent interpretation is that orthogonality changes the decomposition rather than the basic kriging machinery. This is particularly important in multi-fidelity simulation, where one writes $f_u$5. Standard universal kriging can produce implausible $f_u$6, whereas the orthogonal construction enforces $f_u$7 and stabilizes those physically interpretable coefficients without materially changing predictive accuracy (Plumlee et al., 2016).

2. Orthogonality in embedded model-error regression

A second line of work applies orthogonality to calibration and model discrepancy. In the standard KOH formulation,

$f_u$8

data can be fit either by changing $f_u$9 or by fitting the GP bias $p_u$0, which induces confounding between parametric learning and model-error correction. The embedded model-error framework addresses this by representing $p_u$1 in GP weight space,

$p_u$2

where $p_u$3 contains the first $p_u$4 Mercer eigenfunctions and $p_u$5. This finite-dimensional representation permits joint inference on $p_u$6 and $p_u$7 inside nonlinear models (Kuppa et al., 20 Feb 2026).

Orthogonality is defined with respect to the weighted $p_u$8 best-fit parameter

$p_u$9

At $\psi_m$0, the vanishing gradient condition yields

$\psi_m$1

After substituting the weight-space expansion and approximating the integral by quadrature or a data sum, one obtains linear constraints

$\psi_m$2

The joint posterior is

$\psi_m$3

and orthogonality may be enforced either by conditioning the GP prior on $\psi_m$4 or by adding regularization terms $\psi_m$5 to the log-prior, yielding the Regularized OGP formulation (Kuppa et al., 20 Feb 2026).

Because the combined parameter dimension $\psi_m$6 can be large, the framework uses the likelihood-informed subspace (LIS). The nonlinear procedure computes the prior-preconditioned Gauss-Newton Hessian averaged over posterior samples, extracts the leading $\psi_m$7 eigenpairs, decomposes $\psi_m$8, and samples only the $\psi_m$9-dimensional LIS while keeping complementary directions at their prior. The reported effect is a reduction of MCMC cost from $L^2$0 down to $L^2$1. In numerical examples, KOH or conventional GP embedding produces broad or diffuse $L^2$2-posteriors and strong $L^2$3–$L^2$4 confounding, whereas OGP or LOGP+LIS concentrates $L^2$5 near least-squares estimates, makes $L^2$6 and $L^2$7 nearly uncorrelated, preserves meaningful prior predictive behavior under extrapolation, and in an advection–diffusion–reaction PDE reduces maximum absolute error from $L^2$8 to $L^2$9 (Kuppa et al., 20 Feb 2026).

3. Sparse variational orthogonal decompositions

In sparse variational GP inference, orthogonality appears as a decomposition of the prior into a low-rank component and an independent residual. With inducing inputs $$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$0, one writes

$$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$1

where

$$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$2

with covariances

$$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$3

By construction $$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$4, so the prior factorizes into a principal inducing-point process and a full-rank residual GP (Shi et al., 2019).

Sparse Orthogonal Variational Inference augments the principal variables $$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$5 with “orthogonal inducing inputs” $$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$6 for the residual process, defining $$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$7 and

$$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$8

The resulting ELBO contains separate KL terms for $$y(x)=m(x)+z(x),\qquad m(x)=\beta^\top g(x),\qquad z(\cdot)\sim\mathcal{GP}(0,k(\cdot,\cdot)).$$9 and $g(x)$0. In conjugate regression, $g(x)$1 can be collapsed analytically, and when $g(x)$2, the formulation reduces to the standard SVGP bound; the paper states that the orthogonal formulation is at least as tight, and strictly tighter when $g(x)$3. The per-gradient-step cost is

$g(x)$4

and if $g(x)$5, doubling the inducing points only doubles the $g(x)$6-work, rather than giving the $g(x)$7 increase associated with a classical SVGP using $g(x)$8 points. Reported experiments show that SOLVE-GP$g(x)$9 matches SVGP$k$0 at only $k$1 cost on the Snelson toy problem, and that on several UCI datasets SOLVE-GP$k$2 achieves predictive log-likelihoods as good as or better than SVGP$k$3 (Shi et al., 2019).

A related but distinct formulation is the orthogonally decoupled variational GP. There the posterior mean and covariance use different finite subspaces, with $k$4 and

$k$5

The mean is decomposed as

$k$6

so that the $k$7-block captures the component orthogonal to the span of $k$8, while the covariance remains concentrated on the $k$9-block. The KL term remains finite-dimensional, natural-gradient updates become simple in the orthogonal coordinates, and overall time and memory scale as

$m(\cdot)$0

with $m(\cdot)$1 and $m(\cdot)$2. The reported empirical summary states that orthogonal decoupling with natural gradients converges in a few hundred iterations, versus thousands for earlier decoupled schemes, and improves RMSE, NLL, and classification performance over coupled SVGP and previous decoupled bases (Salimbeni et al., 2018).

4. Orthogonal inducing features and inter-domain bases

Orthogonality can also be imposed at the level of inducing features rather than inducing points. For a stationary kernel $m(\cdot)$3, Variational Orthogonal Features uses Bochner’s theorem and selects real-valued functions $m(\cdot)$4 that are pairwise orthogonal in $m(\cdot)$5. Defining

$m(\cdot)$6

one obtains

$m(\cdot)$7

and by rescaling $m(\cdot)$8, $m(\cdot)$9. This removes the $\nabla_\theta m(x;\theta^*)$00 inversion cost associated with dense $\nabla_\theta m(x;\theta^*)$01. With Monte Carlo estimation of the ELBO terms, the stated costs are $\nabla_\theta m(x;\theta^*)$02 for full-rank $\nabla_\theta m(x;\theta^*)$03 and $\nabla_\theta m(x;\theta^*)$04 for diagonal $\nabla_\theta m(x;\theta^*)$05, compared with $\nabla_\theta m(x;\theta^*)$06 in standard sparse GP SVI (Burt et al., 2020).

Orthogonally decoupled inter-domain methods extend the same principle to richer basis functions. In spherical inducing features, a GP is decomposed as $\nabla_\theta m(x;\theta^*)$07, where $\nabla_\theta m(x;\theta^*)$08 is represented by principal inducing variables and $\nabla_\theta m(x;\theta^*)$09 is an orthogonal residual with covariance $\nabla_\theta m(x;\theta^*)$10. For radial kernels on $\nabla_\theta m(x;\theta^*)$11, Mercer's theorem on the sphere gives a spherical-harmonic basis with diagonal prior covariance for the principal features,

$\nabla_\theta m(x;\theta^*)$12

Orthogonal features may then be chosen as additional zonal functions or spherical neural-network activations. The ELBO retains separate KL penalties for the principal and orthogonal blocks, and a minibatch costs

$\nabla_\theta m(x;\theta^*)$13

Because the method uses two smaller blocks rather than a single block of size $\nabla_\theta m(x;\theta^*)$14, it replaces $\nabla_\theta m(x;\theta^*)$15 structure by $\nabla_\theta m(x;\theta^*)$16. Reported experiments include improved ELBOs and test metrics on UCI regression benchmarks and a 1 million-point airline-delay dataset, where $\nabla_\theta m(x;\theta^*)$17 matches the performance of a standard SVGP with $\nabla_\theta m(x;\theta^*)$18 while halving the per-iterate cost (Tiao et al., 2023).

A plausible implication is that “orthogonality” in sparse GP inference serves a dual role: it is simultaneously a representational constraint, preventing redundancy between subspaces, and a computational device, enabling separate optimization of principal and residual components (Burt et al., 2020, Tiao et al., 2023).

5. Orthogonal additive kernels and functional ANOVA

In additive GP models one seeks a decomposition

$\nabla_\theta m(x;\theta^*)$19

but without further constraints this decomposition is not unique. Constants and lower-order effects can be shifted between components, and higher-dimensional interaction terms may absorb lower-order structure. The Orthogonal Additive Kernel addresses this by imposing, for each $\nabla_\theta m(x;\theta^*)$20,

$\nabla_\theta m(x;\theta^*)$21

Starting from univariate base kernels $\nabla_\theta m(x;\theta^*)$22, one constructs constrained kernels

$\nabla_\theta m(x;\theta^*)$23

and then defines

$\nabla_\theta m(x;\theta^*)$24

The OAK kernel is

$\nabla_\theta m(x;\theta^*)$25

Because the inputs factorize, different additive components are mutually orthogonal under the input measure, and the paper states that this construction exactly reproduces the functional ANOVA decomposition (Lu et al., 2022).

This orthogonality has two consequences emphasized in the paper. First, identifiability is restored: under the orthogonality constraints, the decomposition $\nabla_\theta m(x;\theta^*)$26 is uniquely determined by $\nabla_\theta m(x;\theta^*)$27. Second, sparse computation becomes more favorable because each retained component $\nabla_\theta m(x;\theta^*)$28 has effective dimension $\nabla_\theta m(x;\theta^*)$29. For a full $\nabla_\theta m(x;\theta^*)$30-dimensional squared-exponential kernel on Gaussian inputs, the number of inducing points required to reach accuracy $\nabla_\theta m(x;\theta^*)$31 scales as $\nabla_\theta m(x;\theta^*)$32, whereas OAK yields

$\nabla_\theta m(x;\theta^*)$33

The model also admits analytic Sobol sensitivity indices because orthogonality makes the posterior mean variance additive across components (Lu et al., 2022).

Empirically, the paper reports that on Pumadyn, two first-order terms and one second-order term explain $\nabla_\theta m(x;\theta^*)$34 of total variance; on SUSY, SVGP+OAK with $\nabla_\theta m(x;\theta^*)$35 and two-way interactions achieves $\nabla_\theta m(x;\theta^*)$36; on customer churn, OAK+SVGP achieves $\nabla_\theta m(x;\theta^*)$37; and on Pumadyn and Churn, OAK requires $\nabla_\theta m(x;\theta^*)$38–$\nabla_\theta m(x;\theta^*)$39 fewer inducing points than a full SE-GP to reach the same test RMSE or AUC (Lu et al., 2022).

6. Extensions, applications, and limitations

Orthogonality in GP regression is not limited to scalar-response trend correction or sparse variational approximations. In Gaussian orthogonal latent factor processes for large correlated matrices, the model

$\nabla_\theta m(x;\theta^*)$40

uses an orthonormal loading matrix to identify latent factors and to decompose the likelihood into independent projections. Because the factor processes are independent a priori and the loadings are orthonormal, the posterior distribution of the factor processes is independent as well. When the input dimension is one and Matérn kernels with half-integer smoothness are used, a continuous-time Kalman filter yields exact likelihood and posterior computations in linear time in $\nabla_\theta m(x;\theta^*)$41, and the overall complexity is reported as $\nabla_\theta m(x;\theta^*)$42 rather than $\nabla_\theta m(x;\theta^*)$43 for a full multivariate GP (Gu et al., 2020).

Several limitations recur across the literature. In embedded model-error regression, LOGP requires linearization about $\nabla_\theta m(x;\theta^*)$44; ROGP avoids this but introduces non-Gaussian priors, and the regularization parameter $\nabla_\theta m(x;\theta^*)$45 must be tuned to balance constraint strength against sampler efficiency. The same work notes that LIS must be extended to handle the ROGP regularized prior and that alternative inner products, such as Sobolev inner products, would lead to alternative orthogonality constraints (Kuppa et al., 20 Feb 2026). In GOLF, exact linear-cost Kalman computations require one-dimensional inputs; for $\nabla_\theta m(x;\theta^*)$46, the method falls back to direct $\nabla_\theta m(x;\theta^*)$47 likelihood evaluations (Gu et al., 2020).

A common misconception is that orthogonality is a single technique. The literature instead uses the term for several mathematically distinct operations: zero covariance between GP subspaces, integral orthogonality to a mean basis, orthogonality to model sensitivities, zero-mean constraints under an input measure, and orthonormal factor loadings. What unifies these constructions is not a single algorithmic template but the use of orthogonality to separate roles that standard GP models often confound. This suggests that orthogonal GP regression is best understood as a design principle for structuring GP priors, posteriors, and discrepancy terms so that identifiability and computation align with the intended scientific interpretation (Plumlee et al., 2016, Shi et al., 2019).