Gaussian Processes with Bayesian Inference of Covariate Couplings
Authors: Mattia Rosso
Juho Ylä-Jääski
Zheyang Shen
Markus Heinonen
Maurizio Filippone
Venue: Accepted by TMLR ()

Abstract: Gaussian processes are powerful probabilistic models that are often coupled with ARD capable of uncovering the importance of individual covariates. We develop covariances characterized by affine transformations of the inputs, formalized via a precision matrix between covariates, which can uncover covariate couplings for enhanced interpretability. We study a range of couplings priors from Wishart to Horseshoe and present fully Bayesian inference of such precision matrices within sparse Gaussian processes. We empirically demonstrate the efficacy and interpretability of this approach.

Published in Transactions on Machine Learning Research (07/2025)
Gaussian Processes with Bayesian Inference of Covariate Couplings
Mattia Rosso KAUST, Saudi Arabia
Juho Ylä-Jääski Aalto University, Finland
Zheyang Shen Newcastle Universily, UK
Markus Heinonen Aalto University, Finland
Maurizio Filippone KAUST, Saudi Arabia
Reviewed on OpenReview: https: // openreview. net/ forum? id= fameaameeaal.mmmm j03
Abstract
Gaussian processes are powerful probabilistic models that are often coupled with Auto- matic Relevance Determination (ARD) capable of uncovering the importance of individual covariates. We develop covariances characterized by affine transformations of the inputs, formalized via a precision matrix between covariates, which uncover covariate couplings for enhanced interpretability. We study a range of couplings priors from Wishart to Horse- shoe and present fully Bayesian inference of such precision matrices within sparse Gaussian processes. We empirically demonstrate the efficacy and interpretability of this approach.
1 Introduction
Statistical models based Gaussian Processes (GPS) offer attractive modeling choices for various quantita- tive sciences due to their ability to impose functional priors with certain desired characteristics and to carry out principled uncertainty quantification (Rasmussen & Williams 2006). Modeling and inference of GPs has
evolved significantly in the directions of scalability for large data CCtajjar et al. 2017 Hensman et al. 2013)
Wilson & Nickisch 2015), deep learning (Damianou & Lawrence, 2013 Wilson et al.) 2016 Salimboni &
Deisenroth 2017): and generality with autodiff frameworks (Krauth et al.) 2017 Matthews al. 2017).
The choice of the covariance (kernel) function plays crucial role in specifying the function space induced by GPs. This choice is often overlooked by opting for the reputable "default" exponential ARD covariances
Neal 1996), which capture the importance of each covariate, but also assumes an azis-aligned anisotropic data structure, blind to covariate couplings (Matérn 1960p.
In contrast, affine anisotropic covariances are able to explicitly consider the linear dependencies between covariates (Matérn) 1960 Poggio Girosi 1990), which a common feature of real-world data, via the precision matrix A of the Mahalanobis distance (x- - x') A -x)x. A seminal work of Vivarelli Williams (1998) proposes parameterization of the precision based Principal Component Analysis (PCA), while Titsias & SezzaroGGredill (2013) apply mean-field variational inference over factors of such a precision matrix A. Relevant works on affinc-cvarrances GPS include non-stationary extensions (Paciorek & Schervish 2003
and applications to imaging (Kalaitzis 2009) and material sciences (Noack et al. 2020).
In this paper, our goal is to revitalise Mahalanobis distance- based covariances as more interpretable and general alternative "diagonal" ARD covariances, whereby we are able to uncover covariate couplings. This is illustrated in Fig. 1, where we refer these more general types of covariance functions as Automatic Coupling Determination (ACD) covariances. Unlike previous works considering full precision matrices A,

---

Published in Transactions on Machine Learning Research (07 (025)
(a) ARD
(b) ACD
(c) ACD A variance x10
0.2
0.29
0.10
0.0
0.00
0.01
0.1
BREA
0.0
0.1
0.00
-0.1
0.0
0.00
Figure 1: The Automatic Coupling Determination (ACD) covariance reveals rich predictive covariate couplings. Comparison between ARD diagonal precisions Σ " (a) and ACD precision matrix A mean (b) and variance (c) with graph illustrations. assume an element-wise Normal prior on A. The ACD covariance detects that the covariates (5,6,7) are close to redundant on the kinônm dataset.
study fully Bayesian scalable formulation of GPS, where carry out inference over the matrix A, thus obtaining posterior distributions over covariate couplings. An attractive feature of this approach in supervised learning problems is the possibility of obtaining information about covariate couplings that are instrumental in yiclding accurate modeling of the labels.
Our contributions are as follows: (i) GP model that determines covariate couplings through the analysis of the matrix A; (ii) an analysis of sparsity-inducing priors for the matrix from Wishart, Laplace and Horseshoe families; (iii) demonstration of the enhanced explainability of ACD covariances compared to ARD covariances; (iv) the development of fully Bayesian Markov chain Monte Carlo (MCMC) inference scheme of the couplings, while operating in scalable sparse GP framework; and (v) an empirical demonstration of the usefulness of ACD covariances.
2 Background
We consider supervised learning problems with N input-label pairs (X.y) = {(x,: bm) with x, E R and 1/ E R. We model the labels through a latent function f(x), for which we assume a GP prior (Rasmussen) & Williams 2006). Let (f(x)), ( EN denote the collection of latent variables associated with inputs = (xI....,xN)) and II P( be the likelihood function, for which we made an i.id assumption on the realization of the labels yn conditioned on the corresponding f f(xx).
2.1 Gaussian process priors
By imposing a GP prior f(x) ~ GP(O, on the latent function f(x), we are assuming that any subset of these random variables are jointly Gaussian (Rasmussen Williams 2006). The kernel function k(x,, x': 0) determines the properties of the functions that can be drawn from the prior, where O are hyper-parameeers. where k is the
The prior over the realizations of fx) at the inputs X is then p(f AN(0, Kxxjo), (xxjo) N X N covariance matrix obtained by evaluating k(x, x'; 0) at all input pairs xx,x'. For simplicity, we assume zero-mean GPs and omit the conditioning on X.
The posteriors over at inputs x., and inference or optimization over is based on the analysis of the joint
ply, f. 0) =()))))(((0))((()
With Gaussian likelihoods is possible marginalize out leading to a tractable p(y. 0) p(y f p)((p) With non-Gaussian likelihoods, further complications arise due to the lack of conjugacy (Williams & Barber
2

---

Published in Transactions on Machine Learning Research (07 (025))
1998
Opper & Winther 2000 For non-conjugate GP models, popular approaches to sample from the
posterior distribution of covariance parameters and latent variables include Murray & Adams (010)) (slice sampling) and Filippone & Girolami (0144) (pseudo-marginal MCMC).
An overarching issue with GP models is scalability, as these models generally require costly O(N3) operations involving Xxxee inverses. Tinearization techniques based on random features (Rahimi & Recht 2008 were proposed in Gredila et al. (2010), and they were later developed to operate with mini-batches within stochastic gradient optimization and to deep GPs (Cutajjr et al. 2017 In parallel line of work, sparsification techniques based inducing points (Williams & Seeger 2000 Snelson & Ghahramani 2005) were later embedded within variational formulation (Titsias, 2009); and they were extended to mini- batching Hensman et al. 2013: Krauth et al. (017)) In this paper, we consider sparse GPS, and in particular their fully Bayesian version presented in Rossi et al. (2021), where all variables are treated in a Bayesian way and inference is carried out using stochastic gradient MCMC (Chen Zhang 2004).
2.2 Fully Bayesian sparse GPs
We focus on the Bayesian sparse Gaussian process (BSGP) framework (Rossi et al. 2021 but the ACD covariance specifications generally apply to any GP implementations. In sparse GPs, we introduce a set of M inducing variables u (ui...., uM)) which denote the realization of the latent function f(x) at inducing inputs Z {z1,-.. 2v), such that um f(m)) (Candela & Rasmusssn 2005). The definition of the inducing variables implies that they multivariate Gaussian with covariance depending the inducing inputs Z. In summary, the GP prior assumption on the latent function yields the following prior specifications for the latent variables f and w:
p(f, n| X, 2, = pfflu, X, 22.0)p(ulZ,)
(2)
p(u Z. (N(,,
(3)
pf u, X, Z, ~ ~N(AU, K - AK o)).
(4)
where A K w0 This augmented model can be used for modeling tasks by introducing likelihood p(ylf). We assign priors over all remaining variables (0) and Pe(Z), notably including inducing locations Z and covariance parameters 0 (Rossi et al.) (002)). Here we assume that the priors over 0 and Z are conditionally independent, but this can be relaxed if one intends to couple Z and O in the prior. The joint becomes
p00 Z. u,u,f, = P ( u/x, Z, (y)().
(5)
Inference over latent variables and dyper-parameters is intractable, and variational inference allows one recover tractability. The introduction of inducing variables in the definition of the model enables scalable formulation leading to objective that factorizes across data, while requiring algebraic operations with the M x M matrix Keujo instead of Kxx 0- Following Rossi et al. (2021), who study the differences between fully Bayesian treatment of the variational free energy (VFE) and fully independent training conditional (FIrC) approximations, we adopl the latter to approximate the conditional p(f u, X, 2, This allows for
parameter inference over T def {u, Z, with scalable MCMC based stochastic gradients Chen & Zhang,
2004
2.3 Computational Complexity
In the sparse GP approach the matrix factorization of K,n. has complexity O(M"), where M is the number of inducing points. Additionally, multiplying Ky, with the inverse (or performing forward backward sub- stitution with the Cholesky decomposition) of K., incurs a complexity of O(NM2) where N is the number of data points. In the MCMC sampling method proposed by Rossi el al. 22021) mini-batching reduces the complexity of the O(NM²) operation to O(N' where N & N, significantly improving computational efficiency. To be more specific, the complexity for K,,. R. M is O(M2 since each of the M² entries involves matrix-vecton products with the A matrix. Similarly, for the K, E MxN matrix the complexity is O(MND2).
3

---

Published in Transactions on Machine Learning Research (07 (07/2025)
3
Related Works
3.1 Gaussian processes with Automatic Coupling Determination
The possibility of performing kernel-based modeling with the determination of the importance of inputs dates at least back to work on Automatic Relevance Determination ARD (MacKay 1995 Neal 1996). This is usually implemented by scaling input covariates within the calculation of the covariance function by some coefficients which are treated as uyper-parameters and optimized through marginal likelihood optimization.
An extension of this idea involves the use of affine transformations (rotation and stretching) of the covariates; in distancc-based covariance functions, the affine transformation implics the calculation of the so-called Ma- halanobis distance (Vivarelli & Williams, 1998 Titsias & Lazaro-Gredill 2013). In (Vivarelli & Williams, 1998), A is made positive definite by construction through the parameterization UÚ with U upper tri- angular, and it is factorized to gain insights into the dimensionality of a possible low-dimensional latent representation of the inputs. In our work, we consider a similar parameterization for A, but instead of optimizing its factors, we carry out Bayesian treatment, for which we study sparsity-inducing priors. A Also, we propose a PCA-based docompasition of A, which allows us to operate in largo-dimensional input regimes; this is done through the first d principal components Pa of the input covariance, which we use to express A P.A.PP
The work by (Titsias & Lazaro-Gredilla, 2013) considers the parametcrization A = W W without imposing any structure on W except for imposing that W maps the input covariates to a lower dimensional input space. Their work proposes a variational formulation to obtain an approximation to the posterior over A but it does not extensively study priors over W, which is the focus of our work.
3.2 Sparsity-Inducing priors for covariance/precisioo matrices
The literature on Gaussian Graphical Model (GGM) provides studies on sparsity-inducing priors for covariance
and precision matrices. Sparsity can be imposed a structured fashion by considering graph decomposability
(Banerjee & Ghosal, 2014 Lee & Lee 2021 Xiang et al. 2015 Banerjee et al., (021). or through the G-
Wishart prior, which has been introduced
conjugate prior for the precision matrix in Gaussian
framework and it is also suitable for cases where graph decomposability does not apply (Roverato) 2000
2002 Khare& & Rajaratnami
2011 Silva & Ghahramani 2009 van den Boom et al. 2022
Varions approaches have been developed to carry out inference in GGMS, including Gibbs sampling (Khare &
Rajaratnam 2011 Wang 2012) and Laplace approximations (Bancrjoe & Ghosal 2015 Other approaches,
such as those by Gan et al. (2018); Wang (2015), propose spike-and-slab sparsily-inducing priors which
typically complicate posterior sampling. Castillo et al. (2015); Li al. (2017); Sagar al. (2024) study
horseshoc priors, which perform well in practice.
Our work differs from this literature, given that we propose a model for the labels given the inputs, while attempting to uncover some couplings among covariates. The way this is done is by parameterization of the covariance function akin to the precision matrix in a GGM, and in our work we explore both matrix- variate and element-wist priors for such model parameters. In addition, we consider scalable sampling- based approaches to obtain samples from the posterior distribution over these parameters.
4 Bayesian inference of covariate couplings
In this section, after bricfly discussing covariances with Automatic Relevance Determination (ARD) (MacKay
1995 Neal, 1996), which induce some scaling of individual covariales, we present an exlension involving an affine transformation of the covariates revealing couplings among these. We discuss how this achieved by introducing a Mahalanobis distance among inputs with a prccision matrix A, and we show how to treat this in a Bayesian way by imposing matrix-variate and sparsity-inducing clement-wise priors. We term this type of covariance Automatic Coupling Determination (ACD).
4

---

Published in Transactions on Machine Learning Research (07 (075)
4.1 Automatic relevance determination
The design of covariance functions for GP models is an important part of the modeling process. Considering the space of functions f : I² + R, the choice of a covariance (x) f(x') k(x, x'; 0) determines
the prior distribution over f before observing data. A common choice is the Gaussian covariance function (Radial Basis Function (RBF)):
(x, x exp (-26(xx.xx0)
(6)
where d(x, x'; is a parametric distance function between imputs x and x'. This covariance imposes a prior over infinitely differentiable (smooth) functions. Other common covariance functions based on the distance x(x,x'; 0) include the Matérn covariance, exponential, arc-cosine; see, e.g., Shhwe--Taloor & Cristianini (2004)
for in-depth treatment.
The simplest distance form
1 d PPSOTROPIC (x, x 7z (x x') xx-x')
induces an isotropic covariance, all input features are scaled by the same length-scale parameter & and contribute equally to the distance, which assumes spherical data.
Another choice increasing model flexibility introduces covariate-specific length-scales parameters,
d (x, (x 1x-xx) (x-x)
(8)
with Σ = diag({t...., This choice gives rise to covariances suitable for ARD (MacKay
1995
Neal
1996
Intuitively, this definition is built on the assumption that, if a dimension d has a small value of the associated length-scale fa, small changes in the covariate lead to large responses in the target. The covariance induced by this choice anisotropic with an axis-aligned metric acting as a scaling of individual covariates.
4.2 Automatic coupling determination
The family of ARD covariances allows GP models to yield non- non-parametris and probabilistic mappings from inputs to labels, while simultaneously determining the importance of cach covariate if ly's are optimized or inferred. In this paper, we do not limit ourselves to assessing the relevance of each input covariate, but to automatically discover couplings among these in general way which can be readily applied to any distance-basse based covariance function.
We replace the diagonal matrix E containing the inverse length-scales with a full Positive Semi-Definit (PSD) precision matrix A = I-1 in the calculation of distances,
d (x
= (x x') TAxxTXXxXX A(x x') Σ An xf) - x/)).
(9)
yiclding the so-called Mahalanobis distance (Titsias Lazaro-Gredilla 2013 which can be interpreted as distance obtained after alfine transformation (rolation and stretching) of the inputs by the identity
Matérn 1960
Vivarelli & Williams 1998 Kalaitzis 2009)
d( A - A))
(10)
If the underlying distribution of the inpuls x is Gaussian, this operation produces an implicit whilening of the input data yielding x. While the quadratic form in [Eq. 9 has an additive form, the induced functions do not lend themselves an additive function interpretation (Vivarelli Williams] 1998 Duvenaud et al., 2011). Note that if the precision matrix has zero elements A,y O, the distance function ignores the coupling between covariates and j in the calculation of pairwise distances among inputs.
5

---

Published in Transactions on Machine Learning Research (07/2025)
Discriminative Generative modeling The parameterization of ACD covariances has apparent con- nections with Markov Random Fields (MRFs) (see, e-g., Murphy (2023)), ) where the matrix A is used to specify an adjacency structure for set of D random variables (XI,..,, Xp). MRFS offer the possibility to verify conditional independence properties of groups of random variables based on the analysis of A, while placing no other assumptions on their underlying distribution. While is lempling to think of the ACD parameterization of the covariance function as something to be used to draw conclusions on conditional independence among covariates, we are effectively not modeling their distribution. Instead, we are pushing A directly in the definition of the GP prior A() Therefore A assumes the interpretation of precision matrix inducing an affine transformation of the input, which is optimized or inferred based on the marginal likelihood (or & lower bound thereof). Thus the focus is on performing optimization or inference of A to accurately modeling the labels, with the intention of obtaining some indication of the predictive power of couplings of covariates. We leave the modeling of the input through MRFS as an interesting avenue for future work.
4.3 Precision parameterizations
The precision matrix A in the ACD covariance needs to be symmetric and PSD. The PSD constraint in the ARD covariance easy to satisfy, since working with a diagonal version of A only requires to have non-negative elements on its diagonal and consequently, oogttransformaiin of the length-scales is sufficient.
4.3.1 Lower triangular factorization
In the case of the ACD covariance, optimization or inference of A needs to be performed while preserving the PSD constraint, so that it is straightforward to operate with unconstrained optimization/Mcme //MMC sampling. Among all factorizations that can be used to express A, following (Kalaitzis 2009). a natural parameteriza- tion is via the tranglra matrix L,
(11)
This allows us to directly optimize or sample element- wise and recover A. In addition, this parameterization has some computational advantages in calculating Jacobians which are useful within MCMC, as discussed shortly.
Note that, due rotation invariance, the Cholesky factorization is one of infinitely many possible decom- positions of A. Rotation invariance implies that the product of two factors, L and L in case, gives the same A regardless of any rotation (multiplication with a unitary matrix) of the factors. In other words, given any unitary matrix U, LU is also factor of A. This argument indicates that even if we pick a given ordering of the covariates and the Cholesky factor L is a factor of A, then we would obtain the same if we applied any unitary transformation (which includes permutations) to U.
4.3.2 Computational Complexity of the ACD parameterization
Evaluating each element of K and Kuf in the sparse GP framework requires O(D2) computations for ACD covariance functions. This is the case when evaluating element-wise priors and matrix variate priors on A. To be more specific, the complexity for K.n E MMX is O(M2D²) since each of the M entries involves matrix-vector products with the A matrix. Similarly, for the K.ff E MxN matrix the complexity is OMMDD.. The quadratic scaling in D motivates the adoption of low-rank parameterizations of A, which we discuss in the next subsection. Note also that for ARD covariances these costs reduce linear in D.
4.3.3 Low-rank factorizations
The incrcased flexibility offered by the ACD formulation comes at a computational cost, which we need to deal with: going from learning & E RD length-scales in the ARD covariance to learning a full A E DDD matrix. This is why, for problems where the dimensionality D is high but we are at the same time interested in obtaining an informative precision matrix recovering the underlying structure among the D features, we tackle this problem with PCA, similarly to Vivarelli & Williams (1998) and Paciorek & Schervish (2003). In this context, PCA serves as natural way to obtain a low rank transformation of the covariates by retaining

---

Published in Transactions on Machine Learning Research (07/2025)
Table
Summary of precision priors and the range of hyperparameters studied.
Prior
Definition
Parameters
Log pdf
Wishart
p( () W(V, K)
K D, V K-
log C Ea log |Lad TT(AK
Inverse Wishart
p(A) IW(V, K)
K D, V. = Lo
log C - (2K + 1) En log Laal - I
Laplace
p(A,;) C(m,, b)
m = 0,b E (0.01,0.1.1)
log EllA,; - m
Horseshoe
p(Ac)) HLS(T)
01,
log C + 2,2 Als + log Ei A
given amount of explained variance. Focusing on the ACD distance, by applying a projection to the difference between data samples, we obtain:
(x-x) (x
(12)
where Pa the Dxd matrix obtained from the eigendecomposition of the empirical covariance matrix
x PSP
(13)
and X is the centered R NDD input matrix. To obtain Pa we select the d < D columns P corresponding to the d highest eigenvalues from S. sample x E RD can be projected down to R through P.x. As result, we learn a projected version As in this latent representation of the full precision matrix. By applying the projection in Eq. 12 we recover the precision matrix in the original space. Note that in this parameterization, computational scaling of this representation, report comparison between the ACD and ARD parameterizations for Ad-
ARD model for Aa would lead to full precision PaA.P but given the favorable prefer the added flexibility offered by an ACD. In Fig. we
5 Priors over A
As consequence of adopting the BSGP framework we need to specify prior Py (0) over covariance hyper- parameters 0. Dealing with the ACD covariance, the prior is separately placed over both the marginal variance parameter of and on the precision matrix A. While the first is simply a LogNormal distribution with a fixed mcan and variance, the prior distribution over the precision matrix A requires a deeper understanding. First of all, the Cholesky parameterization A LL" in the context of MCMC sampling introduces change of variable. We impose prior probability over non-linear transformation of L, while this is the variable that is actually sampled together with U, Z and o7.
The change of measure induced by the change of variables, requires the determinant of the Jacobian J:
p (vec L) = p(vec J A, vec
(14)
where, for simplicity, we vectorized the matricos and L so that the Jacobian of the transformation is simpler to define. By vectorizing the matrices A and L, its determinant takes particularly convenient form, which can be computed linearly in D (Magnus & Neudecker 1980
log (vec A, L) = log 2D ([(Lau)
(15)
We have identified two different families of priors P(A): (1) matrix-variate distributions over and (2) factorized scalar distributions over the single entries of the precision A defined as p(A) -Ils p(Ar)).
5.1 Matrix-variate priors
Wishart prior When dealing with matrix distributions over PSD matrices, a natural probability distribution to consider is the Wishart distribution. Beside being defined over symmetric PSD matrices.
7

---

Published in Transactions on Machine Learning Research (07/2025)
BOSTON
KIN6NM
POWERPLANT
CONCRETE
ARD
scD-ImmWsttrt
ACD-W shart
ACD Noma
ACD-Laplace (h:0
ACD-Laplace (b-0.00
Horsesose (F-1)
ACD Horsesee (1-0.1
ACD-sorrsso (v-0.01)
2.4
2.6
28
3.0
-1.35
-1.30
-1.25
2,700
2.775
2,000
3.0
3.1
EEG
WILT
BREAST
DIABETES
Rank summary
ARD
shar
ACD-WWllut
ACD-Nomal
AaDaLaplace
ACDD-aapaac (b-0.)
ACD-L MCD-Laplace (b=0.00
ACD-fforeesho (+=1)
ACD- CCDHHorsesee (r-0
(r-0.01)
033
0.4
0.5
0.036
0.0AC
0.045
0.30
0.3b
0.20
0.35
0.50
0.06
0.70
Test MNLL
Test MNLL
Test MNLL
Test MNLL
Pank
Figure 2: The ACD covariances significantly outperform ARD ones select datasets, while being competetive throughout. Test mean negative loglikelihood (MNLL) on both UCI regression benchmarks (lop) classification (bollom) benchmarks with 20% 80% error quantiles (lower is better), and rank summaries (bollom right).
the Wishart prior is a conjugate distribution precision matrices. Considering E RxD the probability density function can be expressed as:
P(A) = VAAV
C (KKD-1) exp (-21((-1A))
(16)
where C = (222) is a constant term, the scale matrix and K 2 D the degrees of freedom parameter. The Bartlett decomposition proves that imposing independent Gaussian priors on the columns of the lower-triangular matrix L (I1,-..,1) Ip) p(l,) N0, V) is equivalent to Wishart distribution over LL" as W(Ap, K). We choose K D degrees of freedom and V D such that the expected precision E[A] = In is identity.
Inverse Wishart Another prior over PSD matricos related to the Wishart distribution is the inverse Wishart. An interesting interpretation stems from the interpretation of A as a covariance matrix in the Fourier domain when Bochner's theorem applied:
KRRP-ACT (Xi, x,;o7.
(17)
TTxx =E,,sV20 cos(H²² x, b) V20, cos x, + b), \~,, Unif 2T).
(18)
Therefore, apart from viewing A as the precision matrix of the kernel, it can also be seen as covariance matrix in the frequency domain, which offers a motivation for using such prior
p(A) = IW
C (K+D++) exp
(19)
The inverse Wishart view can where C -K (K/2))- We set K= D and Ip. translate into more efficient random Fourier approximations.
5.2 Sparsity-inducing priors
Moving away from matrix-variate distributions, it is possible to encourage sparsity in A with an glement-wiss prior. Since might be interested in promoting sparsity in recovering covariance couplings to different
8

---

Published in Transactions on Machine Learning Research (07 (07/2025)
degree than in the contribution of individual covariates, separate the prior over the elements of A follows:
p(A) = ((A)): p(diag A) II p(A.) I((A.)) ()1i)
(20)
where A and diag A are the off-diagonal elements and the DE array of the diagonal elements of A, respec- tively. In this work, we assume weakly informative Gaussian prior on the diagonal of A, while we study different sparsity-promoting prior distributions for A as discussed next.
Laplace. natural way to promote sparse solutions is Li-regularization (cf. graphical lasso in Friedman ct al. (2008), which equivalent a Laplace prior. The expression in (20) becomes:
p(A) = II ij
((A |m, b) II [(Aalll,o")
(21)
where
((A,)/m, b) = = exp
(- | - mll1)
(22)
= = =CCeep
exp
20 (Ais 11)2) 1)²).
(23)
where C and C2 are the normalizing constants. We fix m μ 0, o² 1, and analyze the resulting posteriors for different sparsity coefficients b (lower increases sparsity).
Horseshoe. et al. 2009),
The Horseshoe prior has become popular probabilistic sparsity-inducing prior (Carvalho)
AnjloT N(0,0 a~ c,(11
(24)
where C+ (0, 1) is Half-Cauchy distribution for the local shrinkage o, while is the global shrinkage parameter. The Horseshoe density of a single entry Aij is
mr (Aij)
exp
E
(25)
3
2T
27
where (-) is the exponential integral function that can be approximated by elementary functions. These priors for covariance parameters have been considered, e.g,, in Oh et al. (0019) and Eriksson & Jankowiak 2021).
6 Experiments
We consider eight UCI datasets as a benchmark to evaluate the performance of GP models for regression and classification tasks. We standardize all datasets to zero mean and unit. variance, and report all results with five-fold cross-validation. Following previous works (c.g., Rossi et al. (2021 ), we report test MNLL for all data, and normalized root. mean square error (RMSE) for regression and error for classification tasks.
In all experiments, chose to approximate GPs with 500 inducing points. We ran BSGP for 10, 000 iterations with a step-size of 0.01 and mini-batch of 1, 000 data points. We evaluate performance on test data from 50 samples collected during training after 1, 500 burn-in ilerations and using a thinning of 180. We adopt gradient clipping for numerical stability and to avoid exploding gradients, which we experienced when working with the Horseshoe priors. Note that augmentation techniques might be useful in improving the numerical stability of the sampler for this prior and, more generally, avoid pathologies introduced by hierarchical models (e.g., (Eriksson & Jankowiak 2021 Betancourt & Girolami 2013)).

---

Published in Transactions on Machine Learning Research (07 (07/2025)
KINIWA POWNERPAANT
14 Arra
Figure 3: The precision matrices reveal couplings, redundancies and separabilities. The posterior mean (a) and variance (b) of precision matrices A UCI benchmark dalasets with Horseshoe prior (T = 0.1).
6.1 UCI benchmarks
With the above setup, we report results on UCT benchmarks by considering various choices of priors (See Table 11. For the kernel variance of we placed a Lognormal prior with unit variance and mean 0.05 as in Rossi et al. (2021). The proposed MCMC scheme yiclds good convergence and sampling efficiency, as
illustrated in Appendix C the supplement; see also Fig. 20 and Fig. 21 for insights on the multimodality of the posterior. Fig. 2 shows the comparative performance for the UCI benchmark datasets, including the range between the 20th and 80th percentiles over the different folds, together with rank summary. For the small data sets, we could also run full GPs and we observed a similar trend; we refer the reader to and
F Fig. 15
F Fig. 16 for a direct comparison between full GPS and BSGPs.
Interestingly, different choices of prior and prior hyper-parameters yield comparable performance. A closer inspection indicates that the element-wise Laplace prior performs worst overall, and this might be due the heavy sparsity promoted by this prior (or the lack thereof) for some hyper-parameter settings (Fgg, 5). The clement-wise Horseshoc prior, while promoting sparsity, fares slightly better than the Laplace prior. It is interesting how the inverse Wishart prior, which operates directly on A, promotes some sparsity after all, while offering relatively competitive performance.
6.2 Sparse couplings
Next, we study the precision matrices themselves. Fig. 3 shows the posterior precisions of all benchmark datasets. Notably we see strong dependencies emerging in kin8nm, eeg and wilt datasets, while powerplant. concrete, breast, and diabetes are sparsely diagonal. For the concrete dataset, the posterior precision reveals both a coupling bet ween water and superplaslicizer (a chemical component that improves the usability of concrete without increasing water content) and the importance of the covariate age. In the wilt dataset, beyond the importance of individual covariates, red values in multi-spectral remote sensing images used to classify diseased tree and other land covers. Finally, in the diabetes dataset, we observe a stronger coupling between pregnancies, glucose, and skin thickness, while unveiling the importance of individual covariates such as glucose, skin thickness, and body mass inder. We notice that the standard deviation of the elements on is larger for covariate pairs with large positive/ /neiaaivee partial covariance, while is generally small for pairs that have small partial covariance. This indicates both the relative scaling of uncortainty, and the ficxibility coupling magnitudes.
approach reveals coupling between mean green and
We provide more in-depth look into the dependencies in Fig. [page] that contrasts the precision matrix of the ARD covariance of kin8nm and breast datasets to the posterior mean and standard deviations of the precision matrix of the ACD covariances. The ACD detects that 5th and 6th covariates of kin8nm are close to redundant, and negatively coupled to 7th covariate. Less evindently, in breast we detect coupling chains over covariates such as (0,3,6,8) and (1,6,7), indicating predictive dependencies for this dataset how a different choice of the prior distribution over A can reveal different and sparser structure of the couplings. We visualize these as circular graphs along with the standard deviations of the elements of the precision matrices.
the data. Fig. 9 shows
10

---

80% 100% Published in Transactions on Machine Learning Research (07 (025))
Fig. shows ablation of comparing the posterior mean precision structures from Horseshoe prior with (0.01,0.) 0) on the concrete dataset. The Horseshoe is able to sparsify the entire structure into an ARD-like structure, while higher т = 1 reveals off- off-diagonal dependencies. To obtain more intuition into the couplings, we also visualize the covariate graphs the bottom panel of Fig. that indicate for instance the strong dependence between the 3rd and the 4th covariate. Further illustrations on all UCI data sets for the Wishart prior Fig. 9 inverse Wishart Fig. 10 and Laplace prior with b 0.1 Fig. 11 can be found in the supplement.
Horseshoe (r 0.01)
Horseshoe (r 0.1)
Horseshoe (7
3.27
0.00
-3.27
2
7
6
Figure 4: The sparsity control Horseshoe prior. The posterior mean precision matrices Horseshoe priors concrete dataset with high (left) to low (right) sparsity.
6.3 Sparsification effect
Laplace
100%
Laplace (b=0.01)
Sparsit
Horseshoe (r-0.1)
0.01
Horseshoe (r-0.01)
0.1
1.0
0%
Horseshoe (7=1)
x10
Laplace (b=0.1)
Horseshoe
100%
Normal
InwWishart
Sparsity
50%
0.01
Laplace (b-1)
1.0
Wishart
0%
0%
20%
40%
60%
Threshold of zero x10->
Sparsity [threshold of zero 10 2]
Figure 5: The sparsification of boston dataset. Left: Relationship between precision sparsity and hyperparam- eters. Right: Posterior sparsity from different priors
11

---

Published in Transactions on Machine Learning Research (07/2025)
Fig. shows the sparsity of the posterior precision matrices A in the boston dataset. Surprisingly, the Inverse Wishart prior has an intrinsic sparsifying effect. The Laplace prior sparsifies according to its hyperparameter while, for this dataset, the Horseshoe prior with = 0.1 achieves slightly more sparsity than the Horseshoe prior with == 0.01.
As a conclusion these experiments on UCI, we observe that the Horseshoe prior obtains better performance compared to the Laplace prior and it is competitive with matrix-variate priors. Also, these generally outper- form the ARD covariance. Interestingly, there seems to be some data-dependent effect connecting sparsity and performance; in data sets such boston, high sparsity seems to be associated with good performance, while for others such as eeg it is the opposite. This indicates that sparsity should perhaps be treated as parameer and learned together with the model. We leave this interesting development for future work.
6.4 Low-rank precision matrices
We also look at the effect of low-rank precision matrices, Fig. shows the posterior precision patterns learned by the Wishart prior using a PCA with rank 11, 7 or 3 in contrast to the full rank 13. The performance degrades strongly at ranks lower than 11, which is likely indicative of the intrinsic rank the dataset for this task. lso, is interesting observe how removing components, which reduces the variance explained by the PCA parameterization, affects the ability of the model to recover covariate couplings, indicating that such information is somewhat contained in the removed components. In the same figure, we also report the results obtained by an ARD modeling for the matrix A, in the PCA-based parameterization. As shown in the figure, some structure emerges even in this case, but the parameterization is less flexible, and retaining a large number of components leads to nearly diagonal precisions, preventing the possibility to reveal interesting covariate couplings.
PCA (d 13) ARD
PCA (d-11) ARD
PCA (d-7) 7) ARD
PCA (d-3) ARD
35
30
A80
(d-13) +
(d-11) +
(d-7)
PCA (d-3)
(d-13)
PCA
(d-11)
PCA
(4-7)
1011
min 8 4
10
PCA
(6-3)
2.50
2.15
3.00
3.25
Figure 6: Overly low rank degrades performance. Posterior precisions with Wishart prior of varying rank d of Eq. 12 (top) and corresponding performances (bottom) the bosten dataset. The matrices displayed are of the form P. where Aa is diagonal and full in the ARD and ACD parameterizations, respectively. Pa is obtained by pre-procsssing the training data with PCA (d components). This setup uses 200 inducing points and 3-fold cross validation.
12

---

Published in Transactions on Machine Learning Research (07 2225))
Table MoCap results on subject 09 using GP-ODE with ARD and ACD kernels. The data are projected dimensional latent space and the model is trained with dimension-wise kernels: ARD (5 lengthscales per dimension) and ACD (5 x full precision matrix per dimension). We report Test MNLL and Test Mean Squared Error (MSE) over different folds.
Metric
Method
Subject 09 (short)
GP-ODE-vanilla ARD
1.17 = 0.02
MNL
GP-ODE-vanilla ACD
1.27 = 0.17
GP-ODE-vanila ARD
10.64 + 1.58
MSE(L)
GP-ODE-vanilla ACD
14.72 + 6.95
6.5 Dependencies of motion capture data
We illustrate the capability of ACD covariances to reveal dependencies in a motion capture task, where the subjects internal connectivity is known (igg. 7). We observe a trajectory Y= (21. "x)" E NxD over N timepoints, where y, E R.D represents the noisy observation of subject state x(44) E RD at time 4. The state consists of a total of D 50 measurements across 21 body parts (Fig. 7j. We follow the GP-ODE model (Hcinonen et al. 2018 Hegde et al., 2022), where the state follows an ordinary differential equation (() f(x(t))
with vector-valued GP prior on the differential f: RD + IRD
f(x) ~ GP(O, (x(x,
(26)
where Ko E DxD is an operator-valued kernel. The most straightforward covariance function is separable one xx,x'; 0) = xxx,x' D. where we learn shared precision matrix for all outputs. As an alternative, we also consider a variant K (x,x';)) 0) = (kxx,(((x 01)....... k(z, z'; 0p)}, where each diagonal entry has its own kernel 0,, and its own precision matrix A, associated with output it.
Fig. shows the posterior shared precision mean pooled over the body parts in a human walk cycle. rich pattern of dependencies emerges. For instance, right and left wrists are strongly coupled across the body, while being negatively coupled to each other. The wrists move in large, cyclic and synchronised patterns, while the back and root have little relevance, indicating their smaller movement ranges during walking. Finally, many adjacent body parts are coupled, such as foot and tibia, and wrist and radius. Table 2 shows the performance between ARD and ACD on subject 09, where the likelihoods are similar, but ACD does perform worse in mean square error. The purpose of the experiment was to demonstrate structure learning with standard inference runs, and we did not focus on performance tuning, which ODE models are known to be finicky about (Hegde et al. 2022). As a final note, we emphasize the importance of these results in light of possible simpler alternative analyses to determine covariate couplings; for instance, we could st study the inverse of the empirical covariance of the input covariates directly. For completeness, analysis in the Appendix in Fig. 17 where the sample covariance is also obtained by retaining the first 15 principal components. Unlike this simplistic approach, which one could view as a preliminary analysis of the covariates, the results in Fig. illustrate the emergence of covariate couplings informed by the supervised learning task.
report this
7 Conclusions
In the literature of GPs, covariances equipped with ARD are popular. These materialize with the definition of set of length-scale parameters scaling the inputs, which are then optimized inferred based the marginal likelihood (or anisotropic covariances, where the distance metric among inputs is determined by PSD precision matrix. We showed that this extension provides a framework for metric learning and we discussed some interesting insights on the determination of couplings among covariates. Crucially, thanks to a fully Bayesian scalable formulation of GPs, we can operate with virtually any number of data points and obtain samples from the posterior distribution over such covariate couplings, which can be used to determine the level of confidence in their predictive power.
crrooximation/bouniee
In this work, we revisited more general definition of
13

---

Published in Transactions on Machine Learning Research (07 (07/2025)
(a] A Average Pooling
(b) Skeleton
root
lowerback
head(3)
upperneck(3)
upperback
thorax
0.6
lowerneck(3)
clavicie(2)
lowerneck
upperneck
thorax(3)
humerus(3)
head
rclavicle
upperback(3)
rhumerus
0.2
radius(1)
maaiiu
lowerback(3)
wrist(1)
rwrist
Iclavicle
root(6)
lhumerus
Iradius
-02
femur(3)
twrsst
remur
-0.4
rtibia
bibia(1)
rfoot
-0.6
feeur
Itibia
foot(2)
Ifoot
root
loverback
upperback
thorax
lowerrect
uppernech
head
relavicte
rhumerus
rradius
rwrist
Iclavicle
Ihumeent
Iradius
twrist
reemur
rtibia
rioat
Iffmur-
lubbia lfoot
Figure 7: The ACD covariance reveals highly regular coupling structure from human motion. GP-ODE model trained with shared ACD covariance latent space of 15 dimensions. Panel (a) shows the precision matrix A reporting just the average value for each group of sensors. Panel (b) shows the reference skeleton connectivity.
We also studied priors for the precision matrix A determining the input metric. We showed that element- wise Laplace and Horseshoe priors provide the highest level sparsity, while Horseshoe priors seem offer better performance. Interestingly, the inverse Wishart prior offers higher sparsity than the Wishart prior with overall comparable performance.
In order to address the quadratic scalability with respect to the number of covariates, we also revisited the work by Vivarelli & Williams (1998), which proposes low-dimensional projection of the inputs through PCA, in light of modern scalable GPS and inference.
We are currently investigating an extension our approach whereby the conclusions we can draw from the analysis of are terms of conditional independence statements. In order to do this, we plan extend our model to larget the modeling of both labels and inputs, including prior over the inputs p (Ax))A) in the form of Markov Random Field, where A now determines the conditional independence among covariates.
As future work, it would also be interesting to consider ways in which we could learn the adequate level of sparsily from data by inferring relevant prior hyper-parameters.
References
S. Banerjee and S. Ghosal. Posterior convergence rates for estimating large precision matrices using graphical models. Electronic Journal of Statistics, 8(2):2111 2137, 2014.
S. Banerjee and S. Ghosal. Bayesian structure learning in graphical models. Journal of Multinariate Analysis, 136(C):147 162, 2015.
S. Banerjee, I. Castillo, and S. Ghosal. Bayesian inference in high-dimensional models. arziv:2101.04449 2021.
D. Barry, J.-Y. Parlange, and L. Li. Approximation for the exponential integral (Theis well function). Journal of Hydrology, (1):7)287-229) 2000.
M. J. Betancourt and M. Girolami. Hamiltonian Monte Carlo for hierarchical models. ar arXiv:1312.0006, 2013.
J. Q. Candela and C. E. Rasmussen. Unifying View of Sparse Approximate Gaussian Process Regression. Journal of Machine Learning Research, 11939-1959) 2005.
14

---

Published in Transactions on Machine Learning Research (07/2025)
C. Carvalho, N. Polson, and J. Scott. Handling sparsity via the Horseshoe. In AISTATS, 2009.
I. Castillo, J. Schmidt-Hieber, and A. van der Vaart. Bayesian linear regression with sparse priors. The Annals of Statistics, 43(5):1986 - 2018, 2015.
S. Chen and D. Zhang. Robust image segmentation using FCM with spatial constraints based on new kernel-indueed distance measure. IEEE Transactions on Systems, Man and Cybernetics, Part B, 34(4): 1907 1916, 2004.
K. Cutajar, E. V. Bonilla, P. Michiardi, and M. Filippone. Random feature expansions for deep Gaussian Processes. In ICML, 2017.
A. C. Damianou and N. D. Lawrence. Deep Gaussian Processes. In AISTATS, 2013.
D. Duvenaud, H. Nickisch, and C. Rasmussen. Additive Gaussian processes. In NewrIPS, 2011.
D. Eriksson and M. Jankowiak. High-dimensionss Bayesian optimization with sparse axis- subspaces. Proceedings of the Seventh Conference on Uncertainty in Artificial Intelligence, 2021.
M. Filippone and M. Girolami. Psoudo-marginal Bayesian inference for Gaussian processes. IEFE Transac- tions on Pattern Analysis and Machine Intelligence, 36(11):2214 2226, 2014.
J. Friedman, T. Hastie, and R. Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 2008.
L. Gan, N. N. Narisetty, and F. Liang. Bayesian regularization for graphical models with unequal shrinkage. Journal of the Americon Statistical Association, 114:1218 1231, 2018.
P. Hegde, Ç. Yıldız, H. Lähdesmäki, S. Kaski, and M. Heinonen. Variational multiple shooting for Bayesian ODEs with Gaussian processes. In UAI, 2022.
M. Heinonen, C. Yildiz, H. Mannerström, J. Intosalmi, and H. Lähdesmäki. Learning unknown ODE models with Gaussian processes. In ICML, 2018.
J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian Processes for Big Data. In UAI, 2013.
Kalaitzis. Image inpainting with Gaussian Processes. Master's thesis, University of Edinburgh, 2009.
K. Khare and B. Rajaratnam. Wishart distributions for decomposable covariance graph models. The Annals of Statistics, 39(1):514- 555, 2011.
K. Krauth, F. V. Bonilla, K. Cutajar, and M. Filippone. AutoGP: Exploring the capabilities and limitations of Gaussian process models. In UAI, 2017.
M. Lázaro-Gredilla, J. Siinoncro-Candela C. E. Rasmussen, and A. R. Figuciras- Vidal. Sparse Spectrum Gaussian Process Regression. Journal of Machine Learning Research, 11:1865 1881, 2010.
K. Lee and J. Lee. Estimating large precision matrices via modified Cholesky decomposition. Stalistion Sinica, (31):173-19)) 196, 2021.
Y. Li, B. A. Craig, and A. Bhadra. The graphical horseshoe estimator for inverse covariance matrices. Journal of Computational and Graphical Statistics, 28:747 757, 2017.
D. J. C. MacKay. Probable networks and plausible predictions review of practical bayesian methods for supervised neural networks. Network: Compulation in Neural Systems, 6(3):469, 1995.
J. R. Magnus and H. Neudocker. The elimination matrix: Some lemmas and applications. SIAM Journal on Malria Analysis and Applications, 1(4):422-44 144,222-44, 1980.
B. Matérn. Spatial Variation. Springer, 1960.
15

---

Published in Transactions on Machine Learning Research (07/2025)
A. G. Matthews, M. van der Wilk, T. Nickson, K. Fujii, A. Boukouvalas, P. León- and J. Hensman. GPflow: A Gaussian Process library using TensorFlow. Journal of Machine Learning Research, 18(40):1 6, April 2017.
León-Villagré, Villlgrá, Z. Ghahramani,
K. Murphy. Probabilistic Machine learning: Advanced topic. The MIT Press, 2023.
L. Murray and R. P. Adams. Slice sampling covariance hyperparameters of latent Gaussian models. In J. D. Lafferty, C. K. I. Williams, J. Stawe-Taylor, R. S. Zemel, and A. Culotta (eds.), Advances in Newral Information Processing Systems 23: 24th Anmual Conference on Neural Information Processing Systems 2010. Proceedings of a meeting held 6-9 December 2010, Vancouver, British Columbia, Canada., pp. 1732 1740. Curran Associates, Inc., 2010.
R. M. Neal. Bayesian Learning for Neural Networks (Lecture Notes in Statistics). Springer, 1 edition, August 1996.
M. Noack, G. Doerk, R. Li, J. Streit, R. Vaia, K. Yager, and M. Fukuto. Autonomous materials discovery driven by gaussian process regression with inhomogeneous measurement noise and anisotropic kernels. Scientific reports, 2020.
C. Oh, J. M. Tomezak, E. Gavves, and M. Welling, Combinatorial Bayesian optimization using the graph cartesian product. Curran Associates Inc., Red Hook, NY, USA, 2019.
M. Opper and O. Winther. Gaussian processes for classification: Mean-field algorithms. Newral Computation, 12(11):2655-2666 2684, 2000.
C. Paciorek and M. Schervish. Nonstationary covariance functions for Gaussian Process Regression. In NIPS, 2003.
T. Poggio and F. Girosi. Networks for approximation and learning. In Proceedings of the IEEE, 1990.
A. Rahimi and B. Recht. Random Features for Large-Scale Kernel Machines. In NIPS, 2008
C. E. Rasmussen and C. Williams, Gaussian Processes for Machine Learning. MIT Press, 2006.
S. Rossi, M. Heinonen, F. Bonilla, Shen, and M. Filippone. Sparse Gaussian Processes revisited: Bayesian approaches to indueing-variable approximations. In AISTATS, 2021.
A. Roverato. Cholesky decomposition of hyper inverse Wishart matrix. Biometrika, 87(1):99-112 2000.
Roverato. Hyper inverse Wishart distribution for non-decomposable graphs and its application to bayesian inference for gaussian graphical models. Scandinavion Journal of Statistics, 3913-44114 2002.
K. Sagar, S. Banerjee, J. Datta, and A. Bhadra. Precision matrix estimation under the horseshoe-like prior siorr-penatty dual. Electronic Journal of Slatistics, 18(1):1 - 46, 2024.
H. Salimbeni and M. Deisenroth. Doubly Stochastic Variational Inference for Deep Gaussian Processes. NewrIPS, 2017.
J. Shawe-Taylor and N. Cristianini. Kernel Methods for Pattern Analysis. Cambridge University Press, New York, NY, USA, 2004.
Silva and Z. Ghahramani. The hidden life of latent variables: Bayesian learning with mixed graph models. Journal of Machine Learning Research, 10(41):1187-1238, 2009.
E. Snelson and Z. Ghahramani. Sparse Gaussian Processes using Pseudo-inputs. In NIPS, 2005.
M. Titsias and M. Lazaro-Gredille Variational Inference for Mahalanobis Distance Metrics in Gaussian Process Regression. In NIPS, 2013.
M. K. Titsias. Variational Learning of Inducing Variables in Sparse Gaussian Processes. Tn AISTATS, 2009.
16

---

Published in Transactions on Machine Learning Research (07 (07/2025)
W. van den Boom, A. Beskos, and M. D. lorio. The g- W-Wihhrt weighted proposal algorithm: Efficient posterior computation for Gaussian graphical models. Journal of Compulational and Graphical Stalistics, 31(4):1215 1224, 2022.
F. Vivarelli and C. Williams. Discovering hidden features with Gaussian processes regression. In NIPS, 1998.
H. Wang. Bayesian Graphical Lasso Models and Efficient Posterior Computation. Bayesian Analysis, 7(4): 867 886, 2012.
H. Wang. Scaling It Up: Stochastic Search Structure Learning in Graphical Models. Bayesion Amalysis, 10 (2):351 377, 2015.
C. K. I. Williams and D. Barber. Bayesian classification with Gaussian Processes. IEEE Transactions on Pattern Analysis and Machine Intelligence, 20:1342-1351, 1998.
C. K. Williams and M. Seeger. Using the Nyström method to speed up kernel machines. In NIPS, 2000.
A. Wilson and H. Nickisch. Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP). In ICML, 2015.
A. G. Wilson, Z. Hu, R. Salakhutdinov, and E. P. Xing. Deep Kernel Learning. In AISTATS, 2016.
R. Xiang, K. Khare, and M. Ghosh. High dimensional posterior convergence rates for decomposable graphical models. Electronic Journal of Statistics, 9(2):2828 2854, 2015.
17

---

Published in Transactions on Machine Learning Research (07 (025))
A Experimental details
In this section, we present details to reproduce our experimental campaign. All the experiments were conducted on Google Colab. The code to reproduce the results on the UCI and MoCap datasets can be http:///ww.a..otttttee/////aaaaa Automatic h Dwt.reenatttiinnoooooo and https: lound at grgaaaa...cc/aattrceeesuusssssee
respectively.
BSGP model We use M = 500 inducing points initialized by a k-means algorithm as commonly used in practice and we place Normal prior Pz over the inducing locations Z. For inference, we use an adaptive version of Stochastic Gradient Hamiltonian Monte Carlo (SGHMC) in which the hyperparameters are automatically tuned during a burn-in phase. We set the default hyperparameter of the number of SGHMC steps to K 10. Exclusively for regression datasets with Gaussian likelihood, we employ an Adam optimizer with a learning rate set at 0.01 for optimizing the variance of the likelihood.
ARD kernel We use the Radial Basis Function (RBF) kernel with Automatic Relevance Determination (ARD) placing a LogNormal prior with unit variance and means equal to 1 and 0.05 for the lengthscales and variance, respectively.
ACD kernel We place a LogNormal prior with unit variance and mean 0.05 over the kernel variance of while over the precision matrix A explore wide range of priors. Note that for the Horseshoe prior we use the exponential integral approximation of Barry et al. 2000).
Table 3: Parameter settings for the UCI experiments.
parameter
value
mum. of inducing points
500
mini-batch size num. iterations
1000
10500
step size
0.01
momentumm of burn-in steps
0.05
1500
num. of samples thinning interval
50
180
B Simulated dataset
In this section we carry out experiment using simulated datasets with known underlying precision matrices. In particular we assess the ability of the ESGP model using ACD kernel to recover the true precision A while fitting simple regression problems. Here it's described how the simulated regression datasets are constructed and some experiments conduced to show the behaviour of the model. We consider N input-label pairs y} = {Xn,, W with x,, E RR and yn E R defined as follows:
Xa NO,
K K Kxxli, = 07 ((x, xx,) AAx, - xx,)
(27)
~~ N0 Kxx + on1)
Once that the underlying precision A has been constructed, specifying value for the kernel variance of and another one for the Gaussian noise in observations via ox is sufficient. The regression datasct {X,y) can be used to train a BSGP model by means of the ACD kernel. Through acquiring samples of the precision matrix A, we aim to recover the original underlying precision used to generate the data. A visual insight into this experiment given Fig. 8. Note also that the novariance/proccsioon of the input covariates the identity; therefore, analyzing the covariates by themselves would not reveal any covariate couplings, while our model is capable to retrieve the couplings that were used to generate the labels.
18

---

Published in Transactions on Machine Learning Research (07 (07/202) (025)
True
Normal
Wishart
InvWishart
0.2
0.2
0.2
V mean
2
33
3
16544
5 4
16544
ACD
0.1
7 6
--2.2
4
x10
x10-?
x10-2
0
ACD A std
231
2
3
3
7 1665 5 4
0.5
3
0.0
3
0.0
3
0.0
Figure 8: Underlying sparse precision (top left) compared with and standard deviation of the A samples obtained with different priors. The dataset is made of N = 1000 samples and the labels are obtained according to Eq. 27 setting of 1 and on 0.1
C
Additional results
10
Figure Posterior precision matrix mean (a) and variance (b) with Wishart prior.
ROSTON
KINFMM
POWERPPAR ANIT
CONORETE
RFG
WET
BREAST
DIARETES
AC
Jas
ATO
Figure 10: Posterior precision matrix mean (a) and variance (b) with Inverse Wishart prior.
19

---

Published in Transactions on Machine Learning Research (07/202) 2025)
BOSTON
POWERPLANT
CONCRETE
DIABETES
50 ADD
117
ACO
Figure 11: Posterior precision matrix (a) and variance (b) with Laplace prior b b= 0.1
Laplace 0.01)
Laplace (b 0.1)
Laplace (b 1)
1.35
0..00
-1.35
0
2
2
3
Figure 12: The precision matrices of the wilt dataset using Laplace prior show a progressive level sparsity.
Table 5: Runtime comparison of BSGP training and inference on the boston dataset using ARD kernel and ACD (with Wishart prior) kernel. Values are reported seconds as mean + standard deviation across three different folds.
ARD
ACD
Training time
2082.57 + 16.61
2513.56 + 3.14
Inference time
0.12 + 0.00
0.13 + 0.00
20

---

Published in Transactions on Machine Learning Research (07/202) 2025)
Average
lnvelack [2]
lnverock 15:11
thou 12:15
15
lowomeck(31
head
upperneck())
ixpertect 11271]
13
thorax(3)
Imaa 24]
<dovice
upperback(
radius(1)
modia
lowerback(3)
rootte)
femun(33
fhess 18
mas 1200)
tibia(1)
dost 12:11
15
Punnn
rika
lbot
Figure 13: The full motion capture precision matrix (a), a pooled part-wise version (b) and the reference skeleton connectivity (c).
Predictive posterior sample
-0.6
2-
Predictive posterior sample
0.6
Table 4: UCI datasets used, including number of data- points and dimensionalities.
-.8
Dataset
N
100
200
1.25
1.00
0.75
0.50
boston
506
13
Predictive posterior sample 2
breast
683
9
diabetes
783
8
-0.3
concrete
1,030
8
wilt
4,839
5
-0.5
kin8nm
8,192
8
100
200
-0.5
0.4
0.3
-02
powerplant
9,568
4
sample
Effecsive Somple Wise (ESs)
eeg
45,730
14
1.01
526.80 80
1.00
703-17
651.88
Traces of the mean of the predictive distribu- tion for three test points on boston dataset with Inverse Wishart prior (4 chains, 200 samples represented); the ta- ble reports f and Effective Sample Size (Ess) statistics for each set of 4 chains.
Figure 14:
21

---

Published in Transactions on Machine Learning Research (07 (07/202) 2025)
BOSTON
BOSTON
0-000
2--20
r-000
0.26
0.24
Tes:
0.22
CONCRETE
CONCRETE
0.78
MRMSE
3,30
0.25
Tex:
3,00
AGD
AGO
ACB hww.ehotrrnneeeeeeeeeooe ootfornttoe
ACD
ACD
ACO
ACD
Watert Nerra
InviWetrt Wishat
Nomal
(a)
(b)
Figure 15: Comparison of full GPs (D) vs BSCPs (0) with 200 inducing points on two UCT regression data sets. metrics are MNLL in (a) and normalized RMSE in (b).
BREAST
100 QUL 64B T- 1-0 001
BREAST
Awddl Lw211 wile -nn
EnteRRtt 18%
0.32
3.0
0.10
2.0
Fext turor Hare Lura--aaaee (%)
30.0
25.0
(180
ACD
ADD
A00
A5O
ACD
AED
ACD
ARD
AED
4.3
ACD
ACD
A50
ACD
(a)
(b)
Figure 16: Comparison full GPs (O) vs BSGPS (O) with 500 inducing points on two UCI classification data sets. The metrics MNLL (a) and Error-Rate in (b).
(6) Arca Full
(b) Arce Ahorage Peoking
nor [34g
[d Stideoon
Iovertace 2-5
1890 Sack 5:13
head
upperack(3)
lovamack
doow [1:
lowemeckr3)
clavicie(2)
thorax(3)
humerus(3)
upperbbacke)
redius(3)
mudin
S0:11
lowerback(3)
wrist:1]
idoice (1:33
root(6)
hameraa 1o:30
femur!3)
hadia
30-
dasar
eible(3)
-4
foot(2)
Icst RESS
Woree
Figure 17: (a): The full estimated precision matrix APCA SPTSP22) -1 obtained applying PCA with rank 15 the training dataset; (b): The pooled version of the same matrix; (c): The reference skeleton connectivity.
22

---

Published in Transactions on Machine Learning Research (07 (07/202) 2025)
BSGP Wishart
Full GP Wishart
1.32
-1.32
Figure 18: Comparison of posterior mean of the precision matrix A on the boston dataset with Wishart prior for full GP vs BSGP with 500 inducing points.
ARD
ACD-InoWahert
MCD-WWshart
CCD-Nomaal
ACD Laplace (b-1)
USGP BGP
5
Tixt MNLL
ACD Laplace (b-0.1)
ACD Laplace (b-0.01)
ACD Horseshoe (r-1)
ACD: lorsesoe (--0.1)
ACD T.D-llorcshae (=-0.0)
5
Tex MNLL
wellyw dalaonn Aisi
Wلlll www.
Ju Man www.ehsphheeeeo was
600
1200
600
800
1003
800
1020
800
1000
650
500
1000
haratoon
Ineerione
Iesations
Iterations
Figure 19: MNLL vs iterations for BSGP with 500 inducing points and for full GP on concrete dataset. The plots show value every 10 of the 10, 000 iterations.
23

---

Published in Transactions on Machine Learning Research (07/202) 2025)
-03
-01
66
-01
13
-41
30
-41 10
-11
-13
-12
-03 -22
C
0
3
7
Figure 20: Posterior samples distribution of precision matrix entries for kin8nm dataset with Horseshoe (7 0.1) prior.
Wishart
InwWishart
Normal
Horseshoe
Laplace
0.2
0.2
0.2
PCA2 (4)
0.1
0.
0.3
0.0
0.0
0.0
0.0
-0.1
0.1
-0.1
-0.1
7.5
0.50
5.0-
0.01
0.01
0.25
0.0-
-0.25
-25
PCA-1
PCA-1
PCA-1
PCA-1
PCA-1
Figure 21: PCA representation of vectorized posterior precision matrices. Each point in the 2D space represents a posterior sample (precision matrix). (a) kin8nm dataset, (b) eeg dataset.
24


Reviews:
[{"review_id": "cUSjUAeA2S", "reviewer": "Reviewer_UvS6", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\nTo the best of my knowledge, the use of a dense matrix as the precision matrix\nis novel. Comprehensive numerical experiments are conducted to validate the\napproach.\n\n\nWeaknesses:\n\nI think the contribution of the paper is incremental. In fact, in the abstract the author(s)\nsummarized that the main contribution of the paper was\n\"We develop covariances characterized by affine transformations of the inputs,\nformalized via a precision matrix between covariates, which can uncover covariate couplings\nfor enhanced interpretability.\"\n\nFirst, I think the \"covariances characterized by affine transformations of the inputs\"\nis described in Section 3.2. Here, the author(s) first summarized in Section 3.1 that\nin the literature, the Gaussian RBF kernel usually uses\n\\begin{gather*}\nd^2_{\\mathrm{ISOTROPIC}}(\\mathbf{x},\\mathbf{x}';\\boldsymbol{\\theta})\n=\\frac{1}{\\ell^2}(\\mathbf{x}-\\mathbf{x}')^\\top(\\mathbf{x}-\\mathbf{x}').\n\\end{gather*}\nThen the paper proposed to replace this with the Mahalanobis distance\n\\begin{gather*}\n(\\mathbf{x}-\\mathbf{x}')^\\top\\boldsymbol{\\Lambda}(\\mathbf{x}-\\mathbf{x}'),\n\\end{gather*}\nwhere the precision matrix $\\boldsymbol{\\Lambda}$ could be dense and not diagonal,\nso long as it is positive semi-definite. While this increases the model complexity,\nthis kind of technique is just ordinary Gaussian distribution, and I feel\ndifficult to regard it as a big breakthrough of methodology.\n\nSecond, this paper uses empirical approaches and does not include a theoretical\ninvestigation.\n\nHowever, I acknowledge that the proposed approaches are novel, and that a paper\nof acceptable quality does not necessarily require a theoretical investigation.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper investigates Gaussian processes for nonlinear regression, employing\na dense covariance matrix in place of the conventional RBF kernels. To the best\nof my knowledge, this approach is novel. Furthermore, the sparsity of the\nprecision matrix is exploited through the adoption of tailored prior\ndistributions.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I guess that Equation (11) should be\n\\begin{gather}\n\\sum^D_{i,j}\\Lambda_{ij}(x_i-x_i')(x_j-x_j').\n\\end{gather}\nAt least, please clarify how the current form is derived.", "broader_impact_concerns": "This paper studies a families of Bayesian statistical algorithms.\nI have no concern about potential societal, ethical, environmental, economic,\nor cultural implications of the research beyond its immediate scientific and\ntechnical contributions.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "9Ti17OXk1w", "reviewer": "Reviewer_PLGc", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "To me, the main strength of the paper is the spirit of revitalising precision-based methods to capture dependencies between input features and also their relevance. These sorts of methods, which are linked with seminal works from the 90's like Neal (1996) and MacKay (1995) on the ARD approach, are also current with other approaches in the Bayesian community, (i.e. Laplace approximations with Hessian-based method to approximate the precision matrix of NN's parameter posteriors). In that sense, I'm glad to read a work of this type, but perhaps I'd love to read more details or larger connections (the one with Markov Random Fields is a great insight, but not developed at all).\n\nIn that sense, the continuation of the paper, particularly in sections 4 and 5 is correct, and I've not found particular mistakes or typos. Everything is clear and technically correct in my opinion. A possible critic could be that some subsections are somehow isolated and the flow of explanations is somehow lost sometimes. For instance, the new inclusion of Section 4.4 makes no sense to me, in the sense that yes, computational complexity is important but suddenly sparse approximation cross-covariance matrices are introduced (a bit) out of the blue when in the subsection before the formulation was focusing on a completely different aspect.\n\nExperimental results are OK, and seem sufficient to me, despite it is claimed that a great structure/coupling between features is discovered, when it's not that surprising at times.. I particularly see this in Figure 3, where I don't see really a complex coupling between features, apart from the KIN8 perhaps.  Same with Figure 4 and experiments in subsection 6.2.\n\nThe most interesting (and also classic) experiment is clearly the last one in Figure 7 and developed in subsection 6.5. I am somehow disappointed by the limitations of the method, as it is described in the paragraph before section 4.3 where it's clearly stated: *\"While it is tempting to think of the ACD parameterization of the covariance function as something to be used to draw conclusions on conditional independence among covariates, we are effectively not modeling the distribution of these\"*\n\nLast but not least, the exploration with the 4 types of priors is not extremely novel or of technical relevance to me, as it reminds me a bit of some technical reports of conjugate priors (like the ones K. Murphy sometimes has published). Additionally, I am somehow concerned about the \"full\" Bayesian approach with priors on both theta and Z. In Eq. 6 these ones are independently modelled with marginal priors, but this might not be fully representative of their conditioning. (Anyone who does optimization on both Z and theta at the same time knows how much these ones affect each other). More discussion or taken this point more carefully would be great, at least to provide insights to the reader.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper revitalises precision-based priors for Gaussian process (GP) models. The main motivation comes from the idea of capturing additional information about the coupling of input features. This is particularly inspired by the well-known *automatic relevance determination* (ARD) kernel, where *lengthscale* hyperparameters independently capture the impact of each input feature, but not their particular correlation or degree of conditioning. The work builds on top of key advancements like the use of the Mahalanobis distance metric in the GP regression case, Bayesian GP inference (with the use of priors on hyperparameters and inducing variables) and sparse approximations.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Regarding changes, I'm not keen on requiring major updates, as I think the paper is in a good shape given the current technical contribution, ideas, and empirical results provided. Moreover, I'm aware that the work has already undergone a round of improvement and I acknowledge the good changes included by the authors. To me, it would be great if some of my comments could be considered for a minor improvement of the manuscript, as well as adding more info or updating the following:\n\n- I miss details on the MCMC side, and its convergence. We are not going here for optimization-based solutions for some reason; it's important to know that from reading, but also what the consequences and computational implications of using the MCMC approach in this setting. What if we have a classification model? Does it work equally well, or would it be unfeasible?\n- Section 4.3.2 and details on PCA --- could the authors explain or reformulate a bit better why it makes sense to use it there and what the consequences are.\n- Section 5. Equation (16) and the change of variable >>> I miss the advantage here, and don't understand what is really going on.\n- Figure 3, some comments around and in the caption make me think: \"not always right?\" I would've expected more structure, maybe is worth explaining what should be expected beforehand.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "GubE5ATIQQ", "reviewer": "Reviewer_8AN7", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "As I see it the strengths of the submission include:\n- It is reasonably well written and easy to follow\n- The ACD priors introduced are intuitive and easy to understand.\n- The empirical evaluation is relatively extensive and supports the viability of the introduced priors and inference approach\n\nAs I see it the weaknesses of the submission include:\n- Improvements in predictive performance are limited and no comparison is made to other GP-based methods, e.g. deep GPs.\n- The discussion of related work could be more extensive\n- The discussion of certain empirical details is somewhat limited\n- The methodological novelty of the submission is limited", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors consider fully Bayesian Gaussian Process (GP) regression and classification and introduce priors over the precision matrix in the RBF kernel that controls the GP function prior. The authors consider a few different priors, including (Inverse) Wishart-based and sparsity-inducing Horseshoe-based priors. Inference is done using the inducing-point and SGHMC-based inference procedure from Rossi et al. (2021). In experiments the authors demonstrate that the proposed approach (Automatic Coupling Determination [ACD]) can provide small improvements in predictive performance compared to standard ARD (Automatic Relevance Determination) while simultaneously uncovering (potentially interpretable) sparse structure in the covariates if such is present.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Critical changes:**\n\n1. The lower-triangular-based factorization (Eqn 13) used in much of the submission depends on the order of the covariates. As such the induced prior depends on the order of the covariates. As far as I can tell, there is no discussion of this critical point. Does the ordering of the variables have any impact on the predictive performance and/or inferred sparse structure? This would seem to be a key theoretical and empirical question, and I think it's critical for the authors to address. Some of the priors introduced by the authors may be insensitive to the covariate ordering (e.g. the low rank prior), but in any case possible order dependence should be clearly examined and discussed in detail.\n2. Expand discussion of literature that considers hyperpriors over the hyperparameters of RBF-like kernels, e.g. [A, B, C]. For example, [B, C] both consider Horseshoe-based priors. Certainly the authors could find more relevant references with additional literature search.\n3. Rossi et al describe both VFE- and FITC-based formulations. Please specify which you use.\n4. The authors write that \"E1(·) is the exponential integral function that can be approximated by elementary functions.\" Please clarify if this is what you do in practice. If so what approximation do you use? What is the relevant citation for the approximation, if any? \n\n**Non-critical changes:**\n\n4. The authors write that \"We adopt gradient clipping for numerical stability and to avoid exploding gradients, which we experienced when working with the Horseshoe priors.\" As the authors know, this introduces another source of bias on top of FITC, SGHMC-induced bias, inducing-point-based bias, etc. It would be good to understand why this is necessary in practice. Presumably this is due in large part to to using the exponential integral formulation instead of an expanded formulation on an augmented space (i.e. without integrating out sigma in Eqn 26). Have the authors tried doing inference in the augmented space? This should avoid the potentially extreme curvature that is visible in Eqn 27. Working in the augmented space is e.g. what is done in [C], where care is also taken to use a reparameterized coordinate system (as generally described in [D]). I suspect that the authors could get similar or better empirical performance by using such an approach in favor of ad hoc gradient clipping.\n5. I find it strange that the authors only use 50 MCMC samples to make test-time predictions. Are you sure this is enough to get low-variance estimates of the various metrics you report?\n\n**References:**\n\n[A] \"Pseudo-Marginal Bayesian Inference for Gaussian Processes,\" Maurizio Filippone; Mark Girolami.\n\n[B] \"Combinatorial Bayesian Optimization using the Graph Cartesian Product,\" Changyong Oh, Jakub M. Tomczak, Efstratios Gavves, Max Welling.\n\n[C] \"High-dimensional Bayesian optimization with sparse axis-aligned subspaces,\" David Eriksson, Martin Jankowiak.\n\n[D] “Hamiltonian Monte Carlo for Hierarchical Models,” M. J. Betancourt, Mark Girolami.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)