Residual Correlation as a Diagnostic for Joint-Uncertainty Gains from GP Coregionalisation

Fangqin Zhou; Joaquin Vanschoren

cs.LGSep 24, 2026

Residual Correlation as a Diagnostic for Joint-Uncertainty Gains from GP Coregionalisation

Authors: Fangqin Zhou, Joaquin Vanschoren

Organizations: Eindhoven University of Technology

Abstract

In multi-target regression, correlated targets are often coupled through multi-output Gaussian processes with an intrinsic model of coregionalisation (GP-ICM), assuming that sharing statistical strength improves overall performance. In practice, the benefits are inconsistent. Across the settings studied, we find that the main benefit of coregionalisation is joint uncertainty quantification rather than point prediction. Raw target correlation does not predict when coupling helps; in the separable GP-ICM settings studied here, residual correlation, the cross-target dependence left unexplained by independent per-target predictors, is the strongest predictor of joint-uncertainty gains. We introduce a lightweight diagnostic, $D_{\rm logdet}=-\frac{1}{2}\log\det R_{\rm res}$ , which represents the idealised joint negative log-likelihood (NLL) gain from modelling a full rather than diagonal residual covariance and is computable from independent GPs alone. Across a controlled synthetic study, 16 multi-target benchmarks, and frozen transformer and convolutional neural network representations for keypoint regression, point prediction remains largely unchanged ( $ΔR^2\approx 0$ ). In contrast, $D_{\rm logdet}$ strongly predicts observed ICM NLL improvements ( $ρ_s=-0.83$ , $p<0.001$ ), outperforming heuristics such as the feature-to-sample ratio. We also propose Residual-ICM, which preserves independent marginal variances while adding residual-correlation structure to the joint covariance. Residual-ICM achieves the best average joint NLL among the compared methods, while the diagnostic indicates when covariance coupling is likely to be useful. The diagnostic is specific to global Gaussian residual dependence, the structure captured by separable coregionalisation.

Figures & tables

Predictor of $\Delta\mathrm{NLL}$	Spearman $\rho_{s}$	$p$
Raw target correlation $\bar{\rho}$	$-0.07$	0.80
Feature/sample ratio $p/n$ (coarse proxy)	$-0.46$	0.072
$D_{\text{logdet}}$ (principled)	$\mathbf{-0.83}$	$\mathbf{<0.001}$

Table 1: Predictors of the measured ICM benefit

\Delta\mathrm{NLL}

across 16 datasets. Since

\Delta\mathrm{NLL}=\mathrm{NLL}_{\mathrm{ICM}}-\mathrm{NLL}_{\mathrm{Indep}}

, more negative values favour ICM; hence, a negative correlation means that larger diagnostic values predict larger ICM gains.

Statistic	Value	$p$
$\rho_{s}(D_{\text{logdet}},\Delta\mathrm{NLL})$ (headline)	$-0.83$	$<0.001$
$\rho_{s}(\Delta\mathrm{NLL}/T,\;D_{\text{logdet}}/T)$ (per-target robustness)	$-0.85$	$<0.001$
partial $\rho_{s}$ , controlling $T$	$-0.74$	$0.002$
partial $\rho_{s}$ , controlling $T,p/n,R^{2}_{\mathrm{Indep}}$	$-0.61$	$0.028$
partial $\rho_{s}(T,\Delta\mathrm{NLL})$ , controlling $D_{\text{logdet}}$	$-0.05$	$0.85$

Table 2: Separating residual structure from target dimensionality across the 16 datasets. Rank correlations are against

\Delta\mathrm{NLL}

, so negative values indicate that a larger diagnostic is associated with a larger ICM improvement. Partial correlations are rank-based, with

p

computed at

\mathrm{df}=n-2-k

for

k

controls.

Dataset	$n$	$p$	$T$	$p/n$	$\bar{\rho}$	$D_{\text{logdet}}$	$\Delta\mathrm{NLL}$	$\Delta R^{2}$	$R^{2}_{\mathrm{Indep}}$
scm20d	280	61	16	0.22	0.58	6.83	$-3.590$	$-0.002$	0.546
scm1d	280	280	16	1.00	0.65	3.92	$-1.759$	$-0.006$	0.791
atp7d	280	411	6	1.47	0.63	2.02	$-1.446$	$-0.000$	0.530
atp1d	280	411	6	1.47	0.82	1.57	$-1.145$	$-0.045$	0.715
rf2	280	576	8	2.06	0.42	1.34	$-0.851$	$+0.031$	0.802
rf1	280	64	8	0.23	0.40	1.25	$-1.091$	$+0.001$	0.927

Table 3: Summary statistics for the 16 datasets.

p/n

: feature-to-sample ratio;

\bar{\rho}

: mean pairwise target correlation;

D_{\text{logdet}}

: residual-correlation diagnostic.

\Delta\mathrm{NLL}

and

\Delta R^{2}

denote ICM minus Indep, so negative

\Delta\mathrm{NLL}

favours ICM. Values are averaged over five seeds, cross-validation split, and GP optimiser restart.

Dataset	Method	NLL $\downarrow$	Gap 90 $\downarrow$	Gap 95 $\downarrow$	Sharp $\downarrow$
Tecator	GP (RBF)	$\mathbf{-0.790}$	0.027	0.012	0.165
	BLR	$-0.475$	0.010	0.006	0.174
	Ensemble	0.277	0.073	0.042	0.434
Energy	GP (RBF)	0.262	0.041	0.023	0.206
	BLR	0.323	0.017	0.020	0.316
	Ensemble	0.327	0.039	0.014	0.280

Table 4: Cross-target predictive calibration of the independent GP, Bayesian linear regression (BLR), and deep ensembles, on standardised targets and averaged over both prediction directions per dataset. Standardising targets makes NLL comparable across datasets. NLL, coverage gap, and sharpness are lower-is-better; Tecator and Energy are nonlinear tasks, PHO is near-linear.

Figure 1: Synthetic phase diagram and real-data overlay. Left: joint uncertainty $\Delta\mathrm{NLL}=\mathrm{NLL}_{\mathrm{ICM}}-\mathrm{NLL}_{\mathrm{Indep}}$ (negative favours ICM). Right: point prediction $\Delta R^{2}=R^{2}_{\mathrm{ICM}}-R^{2}_{\mathrm{Indep}}$ (positive favours ICM), which stays near zero across the plane. In both panels, colour encodes the coupling benefit , so that blue indicates a larger benefit ( $-\Delta\mathrm{NLL}$ on the left and $\Delta R^{2}$ on the right); the background shows the mean synthetic benefit over the swept $(T,p/n)$ grid, and points denote real datasets, coloured by which model wins, sized by the effect magnitude, and labelled with their signed values. Each panel’s colour scale is normalised to its own range.

Figure 2: Joint NLL improvement

\Delta\mathrm{NLL}

as a function of the number of coupled targets

T^{\prime}

(mean

\pm

one standard deviation over 20 random target subsets). Gains increase with

T^{\prime}

for atp1d, rf1, and rf2, while wq remains near zero. Point-prediction differences remain negligible ( \Cref tab:datasets_full).

Dataset	Indep	ICM r=1	ICM r=T-1	Residual-ICM	Gate
scm20d (strong)	$\phantom{+}0.00$	$-3.59$	$-5.95$	$-5.41$	$-3.59$
atp1d (strong)	$\phantom{+}0.00$	$-1.14$	$-1.56$	$-1.93$	$-1.14$
rf1 ( $R^{2}_{\mathrm{Indep}}{=}0.93$ )	$\phantom{+}0.00$	$-1.09$	$-1.21$	$-0.86$	$\phantom{+}0.00$
sarcos (ceiling)	$\phantom{+}0.00$	$+0.44$	$-0.06$	$-0.83$	$\phantom{+}0.00$
tecator (ceiling)	$\phantom{+}0.00$	$+0.78$	$+0.62$	$-0.45$	$\phantom{+}0.00$
wq (low struct.)	$\phantom{+}0.00$	$-0.14$	$-0.06$	$-0.10$	$-0.14$

Table 5: Uncertainty gains (

\Delta\mathrm{NLL}

; lower is better) for independent GPs, rank-1 ICM, rank-

(T-1)

ICM, Residual-ICM, and diagnostic-gated coupling, averaged over five seeds. Rank-

(T-1)

is the highest-rank ICM configuration evaluated; the gate uses rank-1 ICM with the leave-one-dataset-out threshold.

Goal	$T$	$p/n$	$R^{2}_{\mathrm{Indep}}$	Recommendation
Selection	any	any	any	Single-output GP-LML
Point prediction	any	any	any	Independent GPs
Joint NLL / calibration	$\geq 6$	$\geq 0.2$	$<0.9$	Residual-ICM; jointly fitted GP-ICM if a coregionalised model is required
Joint NLL / calibration	$2$ – $5$	$0.2$ – $1$	$<0.9$	Residual-ICM; expected gain is usually modest
Joint NLL / calibration	any	$<0.2$	any	Independent GPs unless $D_{\text{logdet}}$ is large
Joint NLL / calibration	any	any	$>0.9$	Gate retains independent GPs; Residual-ICM may still help when $D_{\text{logdet}}$ is large

Table 6: Practical decision map for the separable GP-ICM setting studied here. The

R^{2}_{\mathrm{Indep}}>0.9

ceiling applies only to the diagnostic gate; Residual-ICM can still benefit high-accuracy datasets when substantial residual dependence remains. Recommendations reflect average behaviour rather than guarantees.

Method	Mean $\tau$	Best $\tau$	Worst $\tau$
GP-LML (single-output)	0.706	0.956	0.357
GP-ICM (multi-output)	0.607	0.911	0.286
LogME ( You et al., 2021 )	0.527	0.733	0.143
Capacity ( $-\dim$ )	$-0.392$	$-0.036$	$-0.629$

Table 7: Representation-selection performance measured by Kendall

\tau

between evidence-based rankings and downstream

R^{2}

rankings. Higher values indicate better agreement.

Figure 3: Relationship between

D_{\text{logdet}}

and the observed ICM benefit on frozen image representations for MPII Human Pose. Each point is an encoder–subset condition; marker size indicates the number of coupled targets. The diagnostic remains strongly associated with the uncertainty gain from coregionalisation (

\rho_{s}=-0.70

Appendix figures & tables12 assets

Supplementary material from the paper’s appendix.

Appendix

Dataset	Direction	GP-LML	GP-ICM	LogME
PHO	$V_{c\max}\to J_{\max}$	0.867	0.800	0.733
PHO	$J_{\max}\to V_{c\max}$	0.956	0.911	0.600
Tecator	fat $\to$ protein	0.467	0.444	0.556
Tecator	protein $\to$ fat	0.733	0.689	0.556
Energy	heating $\to$ cooling	0.357	0.286	0.143
Energy	cooling $\to$ heating	0.857	0.514	0.571

Appendix

Table 8: Per-task Kendall

\tau

for representation selection. GP-LML wins or ties in 5 of 6 pairs. GP-ICM never exceeds GP-LML.

Variant	Default $r=1$	No cap	$r=2$	$r=T{-}1$	Optimised $\ell$
Mean $\Delta\mathrm{NLL}$	$-0.616$	$-0.627$	$-0.822$	$\mathbf{-0.928}$	$-0.814$
# datasets improving over default	—	3/16	14/16	8/16	15/16

Appendix

Table 9: ICM hardening ablation over 16 datasets and five seeds. Negative

\Delta\mathrm{NLL}

favours ICM. Relaxing the default rank-1 configuration by removing the noise cap, increasing rank, or re-optimising the shared lengthscale improves mean ICM performance.

	All 16 datasets		$T\geq 3$ (10 datasets)
Method	mean $\Delta\mathrm{NLL}$	$\sum$ positive $\Delta\mathrm{NLL}$	mean $\Delta\mathrm{NLL}$	$\sum$ positive $\Delta\mathrm{NLL}$
ICM, rank 1	$-0.62$	$1.85$	$-0.89$	$1.22$
ICM, rank 2	$-0.82$	$1.43$	$-1.20$	$0.79$
ICM, rank $T{-}1$	$-0.93$	$1.26$	$-1.38$	$0.62$
ICM, hindsight rank	$-0.94$	$1.26$	$-1.40$	$0.62$
Residual-ICM	$\mathbf{-1.13}$	$\mathbf{0.10}$	$\mathbf{-1.54}$	$\mathbf{0.01}$

Appendix

Table 10: Rank sweep against Residual-ICM over five seeds. We report the mean

\Delta\mathrm{NLL}

and the sum of positive

\Delta\mathrm{NLL}

over harmed datasets. “Hindsight rank” selects the best tested rank per dataset using the evaluation

\Delta\mathrm{NLL}

and is therefore an oracle rather than a deployable selection rule.

Dataset	Direction	$\Delta R^{2}$	95% CI	$n$
PHO	$V_{c\max}\to J_{\max}$	$+0.032$	$[+0.015,\;+0.060]$	661
PHO	$J_{\max}\to V_{c\max}$	$+0.031$	$[+0.017,\;+0.048]$	661
Energy	cooling $\to$ heating	$+0.074$	$[+0.047,\;+0.106]$	768
SH	$V_{c\max}\to J_{\max}$	$+0.076$	$[-0.292,\;+0.507]$	50
WUR	$V_{c\max}\to J_{\max}$	$-0.071$	$[-0.249,\;+0.167]$	48

Appendix

Table 11: GP cross-target recovery (

\Delta R^{2}=R^{2}_{\text{GP}}-R^{2}_{\text{linear probe}}

using one target’s encoding to predict the other; higher is better). Bootstrap 95% confidence intervals (CIs; 2000 samples). Significant at large

n

(PHO, Energy); indistinguishable from zero at

n\leq 50

(SH, WUR).

Variant	mean $\Delta\mathrm{NLL}$	improved	$\rho_{s}$ with neural $D_{\text{logdet}}$	mean $\Delta R^{2}$
post-hoc residual	$\mathbf{-0.80}$	13/16	$\mathbf{-0.88}$ ( $p<10^{-4}$ )	$0$ by construction
jointly learned full	$-0.34$	11/16	$-0.25$ ( $p=0.36$ )	$-0.066$ (12/16 reduced)

Appendix

Table 12: Neural covariance comparison over 16 datasets and three seeds.

\Delta\mathrm{NLL}

is relative to the diagonal deep ensemble; the diagnostic is computed from that ensemble’s out-of-fold predictive-standardised residuals.

Figure 4: Synthetic $(T,p/n)$ phase diagram of the joint-NLL benefit at five target-correlation levels $\rho$ . Cell values report $\Delta\mathrm{NLL}=\mathrm{NLL}_{\mathrm{ICM}}-\mathrm{NLL}_{\mathrm{Indep}}$ , so negative values favour ICM; colour encodes the equivalent benefit $-\Delta\mathrm{NLL}$ , with blue indicating larger coupling benefit. The high- $T$ , high- $p/n$ structure is absent at $\rho\leq 0.3$ , where little residual correlation remains to exploit, and emerges at larger $\rho$ without changing its overall shape.

Regime	global $D_{\text{logdet}}$	Residual-ICM (global $\mathbf{R}$ )	oracle region-wise $\mathbf{R}$
stationary	$1.76$	$-1.86$	$-1.86$
nonstationary	$0.08$	$-0.08$	$\mathbf{-1.18}$

Appendix

Table 13: Input-dependent residual correlation (

T=4

\rho=0.9

, mean over three seeds). The nonstationary regime carries the same local dependence as the stationary one but with opposite signs in the two halves of the input space. A global correlation matrix cancels it; an oracle region-wise matrix recovers it.

Residual coupling	$D_{\text{logdet}}$	Pearson $\|\rho_{p}\|$	Dist. corr.	Gaussian $\Delta\mathrm{NLL}$	$k$ -NN cond. $\Delta\mathrm{NLL}$
independent	$0.00$	$0.02$	$0.10$	$-0.00$	$+0.05$
linear	$0.22$	$0.60$	$0.56$	$\mathbf{-0.23}$	$-0.18$
nonlinear	$0.00$	$0.01$	$0.33$	$\phantom{-}0.00$	$\mathbf{-0.30}$

Appendix

Table 14: Nonlinear residual-dependence control (

T=2

, three-seed mean). The nonlinear case is Pearson-uncorrelated but detected by distance correlation; Gaussian covariance coupling gives no gain, whereas the nonparametric conditional model does.

Figure 5: Controlled collinearity sweep (

T=24

targets, equicorrelation

c

, 20 seeds). Plain

D_{\text{logdet}}

(red) saturates at its numerical floor in the rank-deficient regime

n<T

(left), becoming constant across

c

and unable to order collinear conditions. The Ledoit–Wolf estimate (blue) tracks the true value (dashed) monotonically. This explains the AFLW saturation and its mitigation by shrinkage.

Regime	#conditions	plain	Ledoit–Wolf
Tabular benchmarks	16	$\mathbf{-0.81}$	$-0.60$
MPII (moderate residual)	50	$-0.70$	$-0.71$
AFLW (collinear, saturated)	50	$-0.69$	$\mathbf{-0.83}$

Appendix

Table 15: Rank correlation

\rho_{s}(D_{\text{logdet}},\Delta\mathrm{NLL})

under the plain empirical estimator and a Ledoit–Wolf shrinkage estimator, evaluated on the 16 tabular datasets and the MPII/AFLW keypoint conditions using one seed-0 run for a like-for-like comparison. Shrinkage mitigates AFLW saturation, leaves MPII essentially unchanged, but degrades the tabular ranking.

Check	$\mathbf{R}_{z}$ (paper)	$\mathbf{R}_{\text{raw}}$ (sensitivity)
Spearman $\rho_{s}$	$-0.83$	$-0.85$
Pearson $\rho_{p}$	$-0.91$	$-0.93$
Bootstrap 95% CI (over datasets)	$[-0.99,\,-0.47]$	$[-0.98,\,-0.51]$
Leave-one-dataset-out range	$[-0.92,\,-0.80]$	$[-0.91,\,-0.81]$
Per-seed range (5 seeds)	$[-0.89,\,-0.73]$	$[-0.87,\,-0.67]$
Permutation $p$ (one-sided)	$<10^{-4}$	$<10^{-4}$

Appendix

Table 16: Robustness of

\rho_{s}(D_{\text{logdet}},\Delta\mathrm{NLL})

across the 16 datasets under predictive-standardised (

\mathbf{R}_{z}

) and raw (

\mathbf{R}_{\text{raw}}

) residual correlation.

Ceiling $c$	$0.80$	$0.85$	$0.90$	$0.926$	$0.93$	$0.935$	$1.0$
Gate mean $\Delta\mathrm{NLL}$	$-0.588$	$-0.641$	$\underline{-0.641}$	$-0.641$	$\mathbf{-0.709}$	$-0.682$	$-0.682$

Appendix

Table 17: Sensitivity of the diagnostic gate to the accuracy ceiling

c

. Mean

\Delta\mathrm{NLL}

is reported over all 16 datasets at the fixed operating point

\tau=0.5

; lower values are better. The gate is stable across

c\in[0.85,0.926]

, with the main text’s

c=0.9

inside this plateau. Residual-ICM achieves a lower mean

\Delta\mathrm{NLL}

(

-1.13

) than the gate for every ceiling value.

Explore similar work

May 13, 2026cs.LG

Safe Bayesian Optimization for Uncertain Correlation Matrices in Linear Models of Co-Regionalization

This paper extends safety guarantees for multi-task Bayesian optimization with uncertain co-regionalization matrices from intrinsic co-regionalization models to linear models of co-regionalization. The latter allows for more flexible modeling of the inter-task correlations by composing multiple features. We derive uniform error bounds for vector-valued functions sampled from a Gaussian process with a linear model of co-regionalization kernel. Furthermore, we show the potential performance gains of linear models of co-regionalization in a numerical comparison on a safe multi-task Bayesian optimization benchmark.

Jannis Lübsen, Annika Eichler

Institute of Control Systems, Hamburg University of Technology, Hamburg, Germany · Institute of Control Systems, Hamburg University of Technology · Deutsches Elektronen-Synchrotron, Hamburg, Germany

Bayesian Optimization Gaussian Process Full-Rank Correlation Matrices Robust Optimization

Ask PDF

Sep 23, 2026stat.ML

Multitask Regression with Pairwise Fusion

We study multitask regression when coefficient sharing can differ by predictor. For a given predictor, many tasks may have the same coefficient while a few differ, and the exceptional tasks need not be the same for another predictor. We describe this structure by two quantities: the number of active predictors and the total number of task coefficients that differ from the most common value for their predictor. We estimate the coefficient matrix by penalizing all pairwise coefficient differences across tasks, with an additional group penalty when predictor selection is needed. The resulting upper and lower bounds have the same dependence on these two quantities. We also consider the stronger setting in which a large set of tasks shares one entire coefficient vector. Under explicit sample-size conditions, the same pairwise estimator pools those tasks exactly, while allowing the remaining tasks to differ. Simulations and household energy data illustrate the transition between broad sharing and task-specific coefficients.

Xiaodong Li, Zhentao Li

Department of Statistics, University of California, Davis, Davis, CA 95616, USA

Coefficients Collaborative Fusion

Ask PDF

Aug 12, 2026cs.LG

A Factor Graph Approach to Scalable Multi-Output Gaussian Process Regression

Multi-output Gaussian process regression scales cubically in the number of observations times outputs, and dense kernel-matrix methods need bespoke handling whenever different outputs are observed at different inputs. We express multi-output Gaussian process regression as a Forney-style factor graph in which a nearest-neighbor chain orders a fixed candidate set of

C

inputs into a one-dimensional sequence. Along this chain, latent Matérn processes evolve through linear-Gaussian transition factors, while the linear model of coregionalization mixes

L

latent processes into

D

outputs through a deterministic mixing factor and per-output scalar observation factors. Posterior computation reduces to exact Gaussian message passing on the chain at cost

\mathcal{O}(C(DL^2 + L^3))

after chain construction, and missing observations omit their local factor without any covariance-matrix restructuring. The formulation therefore scales in the number of data samples and in the rate of missing observations, while remaining best suited to candidate sets in low input dimension.We compare the factor-graph formulation against an exact kernel-matrix baseline, a sparse-variational inducing-point baseline, and a nearest-neighbor baseline on a synthetic input-dimension sweep and on electricity time series forecasting. At low input dimension the factor-graph posterior tracks the exact kernel-matrix posterior closely, and the gap grows gradually as input dimension increases while staying competitive with both approximate baselines. On the electricity time series our factor-graph formulation matches all three baselines in forecast accuracy while scaling linearly in the number of data points, where the exact kernel-matrix method becomes infeasible and the inducing-point baseline remains substantially slower.

Wouter W. L. Nuijten, Esther G. van Pelt, Albert Podusenko +2

Eindhoven University of Technology & Lazy Dynamics B.V., the Netherlands · Eindhoven University of Technology, Eindhoven, the Netherlands · Lazy Dynamics B.V., the Netherlands

Gaussian Process Factor Graphs

Ask PDF