Variational Elliptical Processes
Authors: Maria Margareta Bånkestad
Jens Sjölund
Jalil Taghia
Thomas B. Schön
Venue: Accepted by TMLR ()

Abstract: We present elliptical processes—a family of non-parametric probabilistic models that subsumes Gaussian processes and Student's t processes. This generalization includes a range of new heavy-tailed behaviors while retaining computational tractability. Elliptical processes are based on a representation of elliptical distributions as a continuous mixture of Gaussian distributions. We parameterize this mixture distribution as a spline normalizing flow, which we train using variational inference. The proposed form of the variational posterior enables a sparse variational elliptical process applicable to large-scale problems. We highlight advantages compared to Gaussian processes through regression and classification experiments. Elliptical processes can supersede Gaussian processes in several settings, including cases where the likelihood is non-Gaussian or when accurate tail modeling is essential.

Jens Sjölund Published in Transactions on Machine Learning Research (08/2023)
Variational Elliptical Processes
Maria Bånkestad
maria. bww.kettadddl.....m
RISE Research Institutes of Sweden and Uppsola University, Sweden
cens.sjolund@ii.c
Department of Information Technology, Uppsala University, Sweden
Jalil Taghia
jubil. www.llaahh@@rrcsoonnccommm
Ericsson Research, Sweden
Thomas B. Schön
shomas.schon@ii wu.se
Department of Information Technology, Uppsala University, Sweden
Reviewed on OpenReview: https: openreview. / net/ forum? id= dd=gy19Toqbd48
Abstract
We present elliptical processes-4 family of non-parametrie probabilistic models that sub- sumes Gaussian processes and Student's processes. This generalization includes range of new heavy-tailed behaviors while retaining computational tractability. Elliptical processes based representation of elliptical distributions as continuous mixture of Gaussian distributions. We parameterize this mixture distribution we train using variational inference. The proposed form of the variational posterior enables a sparse variational elliptical process applicable to large-scale problems. We highlight ad- vantages compared to Gaussian processes through regression and classification experiments. Elliptical processes can supersede Gaussian processes in several settings, including cases where the likelihood is non-Gaussian or when accurate tail modcling is essential.
spline normalizing flow, which
1 Introduction
EP-posterior
Systems for autonnomou decision-making are increasingly dependent on predictive models. To ensure safcty and reliability, is essential that these models capture uncertainties and risks accurately. Gaussian processes (GPs) offer framework for probabilistic modeling that is widely used partly because it provides uncertainty estimates. However, these estimates are only reliable to the extent that the model is correctly specified, i.e., that the assumptions of Gaussianity hold truc. On the contrary, heavy-tailed data arise in many real-world applications, including finance (Mandelbrot 1963) signal processing (Zoubir et al. 2012). and geostatistics DDggle et al. 1998). We use combination of normalizing flows and modern variational inference techniques to extend the modeling capabilities of GPs to the more general class of elliptical processes (EPs).
G7-posterior
0.0
Figure 1: Posterior distributions of elliptical and Gaussian processes with equal kernel hyperpa- rameters and covariance. The shaded areas are credible intervals of the posterior processes. The elliptical credible intervals are wider due to the process' heavier tail.
Elliptical processcs. Elliptical procosses subsume Gaussian processes and Student's t processes (Rasmussen & Williams 2006 Shah et. al.] 2014). They based

---

Published in Transactions on Machine Learning Research (08/2023)
the elliptical distribution-a scale-mixture of Gaussian distributions that is attractive mainly because it can describe heavy-tailed distributions while retaining most. of the Gaussian distribution's computational tractability (Fang et al. 1990). We use a normalizing flow (Rezende & Mohamed] 2015 Papamakarios et al, 2021 to model the continuous scalc-mixture, which provides an added ficxibility that can benefit a range of applications. We explore using elliptical processes as both a prior (over functions) and a likelihood, as well as the combination thereof. We also explore using EPs as a variational posterior that can adapt its shape to match complex posterior distributions.
Variational inference. Variational inference is tool for approximate inference that uses optimization to find a member a predefined family of distributions that is close to the target distribution (Wainwright ct al. 2008 Blei et al.) 2017). Significant advances in the last decade have made variational inference the method of choice for scalable approximate inference in complex parametric models (Ranganath et al. 2014
Hoffman al. 2013 Kingma & Welling 2014 Rezende et al. 2014).
It is thus not surprising that the quest for more expressive and scalable variations of Gaussian processes has gone hand-in-hand with the developments in variational inference. For instance, sparse GPs use variational inference to select inducing points to approximate the prior (Titsias 2009). Inducing points are a common building block in deep probabilistic models such as deep Gaussian processes (Damianou & Lawrence 2013
Salimbeni et al.] 2019) and can also be applied in Bayesian neural networks (Maroñas et al., 2021 Ober
& Aitchison 2021). Similarly, the combination of inducing points and variational inference enables scal-
able approximate inference in models with non-Gaussian likelihoods (Hensman al.) 2013). such as when
performing GP classification (Hensman et al. 2015 Wilson et al., 2016b).
However, the closeness of the variational dis- tribution to the larget distribution is bounded by the flexibility of the variational distribution. Consequently, the success of deep (neural net- work) models has inspired various suggestions on flexible yet tractable variat ional distributions, of- ten based on parameterized transformations of simple base distribution (Tan et al. 2016 In particular, models using a composition of in- vertible transformations, known as normalizing flows, have been especially popular (Rezende Mohamed 2015 Papamakarios et al. 2021).
plx,, X2)
R
0.3
0.2
0.1
0.0
-2
0
-2
X1
Figure 2: Left: A contour plot of elliptical two- dimensional, correlated distribution with zero means. The name derives from its elliptical level sets. Right: Three examples of one-dimensional elliptical distributions with zero means and varying tail-heaviness. Elliptical dis- tributions are symmetric around the mean [[[] μ..
Our contributions. We propose adapta- tion of elliptical distributions and processes in the same spirit as modern Gaussian processes. Our ED construction similar to that Maune- Deschamps et al. (2017) but aims to turn EPs into a practical modeling tool. Constructing el- liptical distributions based on a normalizing flow provides high degree of flexibility without sacrificing computational tractability. This makes it possible to sidestep the "curse of Gaussianity", and adapt to heavy-tailed behavior when called for. We thus foresee many synergies between EPs and recently developed GP methods. We make first exploration of these, and simultaneously demonstrate the versatility of the clliptical process as a model for the prior and/or the likelihood, or as the variational postcrior. In more detail, our contributions are:
construction of the elliptical process and the elliptical likelihood as continuous scale-mixture of Gaussian processes parameterized by normalizing flow;
à variational approximation that can either learn an elliptical likelihood or handle known non- Gaussian likelihoods, such as in classification problems;
formulating sparse variational approximation for large-scale problems; and
describing extensions to heteroscedasttic data.
2

---

Published in Transactions on Machine Learning Research (08/2023)
2 Background
This section provides the necessary background on elliptical distributions, elliptical processes, and normal- izing flow models. Throughout, we consider a regression problem, where we are given a set of N scalar observations, (u......NN' at the locations X = [#1...., xN]", where æ, is D-dimensional. The measurement Yi is assumed to be noisy measurement, such that,
V/2 f(x;) +
(1)
where 6; is zero mean iidd. noise. The task is to infer the underlying function f: RD - R.
2.1 Elliptical distributions
Elliptical processes are based on elliptical distributions (Figure 2 which include Gaussian distributions
well as more heavy-tailed distributions, such as the Student's 's distribution and the Cauchy distribution.
The probability density of a random variable Y E R- that follows the elliptical distribution can be expressed as
p(u; 7) CN. Σ 1/2 11111//////
- gN m);
where u := (y EE- (1) is the squared Mahalanobis distance, μ the location vector, Σ is the non-negative definite scale matrix, and CN.7 is normalization constant. The density generator gNv (u: 7)) a non-negative function with finite integral parameterized by 7, which determines the shape of the distribution.
Elliptical distributions consistent, i.e., closed under marginalization, if and only if p(u; 17) is scale- mixture of Gaussian distributions [Kno 1994 The density can be expressed as
p(u; n)
e plE: me)dE,
(3)
2
using a mixing variable & ~ p(E; na). Any mixing distribution p(E: me) that is strictly positive can be used define consistent elliptical process. In particular, we recover the Gaussian distribution the mixing distribution is Dirac delta function and the Student's distribution if it is scaled inverse chi-square distribution. For more information the elliptical distribution, Appendix A
2.2 Elliptical processes
The elliptical process is defined analogously to Gaussian process as:
Definition 1 An elliptical process (EP) is a collection of random variables such that every finite subset has consistent elliplicul distribution, where the scale matrix is given by covariance kernel.
This means that EP is specified by a mean function ((x), scale matrix (a kernel) k(æ,, x') and the mixing distribution p(E; ng). Since the EP is built upon consistent elliptical distributions, it is closed under marginalization. The marginal mean le is the same as the mean for the Gaussian distribution, and the covariance Cov[Y] = E [E] Σ where Y is an elliptical random variable, Σ the covariance for a Gaussian distribution, and & the mixing variable.
Formally, stochastic process (X, I LE T} on probability space (Q,F, consists of random maps X, St, t E T, for measurable spaces (S,, t E T (Bhattacharya & Waymire 2007). We focus the setting where S = R and the index set T is subset of R in particular, the half-line [0, co). Due to Kolmogorov's extension theorem, we may construct the EP from the family of finite-dimensional consistent, elliptical distributions, which is due to the restriction to 5 = (which is a Polish space) and Kano's characterization above.
Note that the above definition of the clliptical process is essentially the same as that in Maume- Deschamps et al. (2017) but makes the use of covariance kernel explicit. (This definition also appears
3

---

Published in Transactions on Machine Learning Research (08/2023)
in Bånkestad et al. 20200) which is early preprint version of the present article.) Deceptively, there are also definitions of elliptical processes that differ in important ways from ours, notably ones based on Lévy processes [Bingham & Kiesel 2002) which assume independent increments.
Identifiability. When using GP for regression clas- sification, it is typically assumed that the data origi- nate from & single sample path, representing a realiza- tion from the GP. However, an elliptical process intro- duces a hierarchical model wherein & first sampled from the distribution p(:; ng). followed by sampling from GP(f: fu, Kg). In this structure, inferring the mixing dis- tribution p(:; me) from single path is impossible since we only have one observation from p(E; ne). In other words, to infer the mixing distribution p(E; mg), we must have multiple paths drawn from the same EP. Therefore, to infer an elliptical prior from data, the dataset would have contain multiple sample paths, such as multiple time series, all originating from the same prior.
E
ft
E
W
N
N
Ws
(b)
N
(a)
(c)
Figure 3: Graphical models of (a) the elliptical
likelihood, (b) the EP-prior, and (c) the EP with independent. elliptical noise, where w is sampled from the likelihood mixing distribution p(w; nw).
Prediction. To use the EP for predictions, we need the predictive mean and covariance of the corresponding elliptical distribution. The predictive discribution is guaranteed to be consistent elliptical distribution but not necessarily the same the original one the shape depends on the training samples. (Recall that consistency only concerns the marginal distribution.) The noise-free predictive distribution can be derived analytically (see Appendix B but to account for additive elliptical noise, we will instead solve it by approximating the posterior p(E) y; ne) with a variational distribution g(:: (P). The approximate inference framework also lets us incorporate (non- coon-Gussian)) noise according to the graphical models in Figure 3
We aim to model mixing distributions that can capture any shape of the elliptical noise in the data. Thus, to improve upon the piecewisc constant parameterization in carlier version of this work (Bänkestad et al. 2020), use normalizing flows: class of methods for learning complex probability distributions.
2.3 Flow based models
Normalizing flows are a family of generative models that map simple distributions to complex ones through a series of learned transformations (Papamakarios et al. 2021). Suppose we have a random variable æ that follows an unknown probability distribution (x). Then, the main idea of normalizing flow is express æ as transformation T, of variable z with known simple probability distribution P=(z). The transformation T, has to be bijective, and it can have learnable parameters 7. Both T and its inverse have to be diseerentiable. A change of variables obtains the probability density of æ:
(x( = (2) det (T,
(4)
dz
interested in modeling the mixing distribution. In particular, we use linear rational spline flows (Dolatabadi et al. 2020 Durkan et al. 2019), wherein the mapping T, is an elementwise, monotonic linear rational spline: a piecewise function where each piece is a linear rational function. The parameters are the number of pieces (bins) and the knot locations.
We focus one-dimensional flows since
3 Method
We propose the variational EP with elliptical noise, where the variational EP can learn any consistent ellip- tical process, and the elliptical noise capture any consistent elliptical noise. The key idea is to model the
4

---

Published in Transactions on Machine Learning Research (08/2023)
mixing distributions with normalizing flow. The joint probability distribution of the model (see Figure 3cl is
N
p(u, f,W, w.w. n) =(ffe: My (e(: ne)
II
p(yel (iww:) nw).
likelihood
Here, p(f m5) ~ (0, K is regular GP prior with the covariance kernel K containing the parameters nf p(:: n2) is the process mixing distribution, and p(w; nw) is the noise mixing distribution.
To learn the mixing distributions plš: ni) and p(w: nw) by gradient- based optimization, they need to be differentiable with respect to the parameters πソ and Tw in addition to being flexible and computationally efficient evaluate and sample from. Based on these critcria, spline flow (Section 2.3 natural fit.
We construct the mixing distributions by transforming a sample from a slandard normal distribution with spline flow. The output of the spline flow is then projected onto the positive real axis using differentiable function such as Softplus or Sigmoid.
In the following sections, we detail the construction of the model and show how to train it using variational inference. For clarity, we describe the likelihood first before combining it with the prior and describing a (computationally efficient) sparse approximation.
3.1 Likelihood
describes the mcasurement noise 6 (Equation 1 distribution of the independent elliptical likelihood is
By definition, the likelihood (Figure 3a
The probability
ple; o, nw) mw)== N(e; / wi)) n
(6)
where a > 0, and can be set to unity without loss of generality. In other words, the likelihood a continuous mixture of Gaussian distributions where, e.g,, E, follows Student's t distribution if w is scaled inverse chi-squared distributed.
Parameterization. We parameterize p(w; nw) as spline flow,
ØT mw)
p(w; nw) p(()
(7)
ac
although it could be, in principle, any positive, finite probability distribution. Here, p(S) ~ W0 1) is the base distribution and w = mww: represents the spline flow transformation followed by a Softplus transformation to guarantce positivity of w. The ficxibility of this blow- construction lets us capture a broad range of elliptical likelihoods, but we could also specify an appropriate likelihood ourselves. For instance, using categorical likelihood enables EP classification; see Section 4.5
Training objective. Now, suppose that we observe N independent and identically distributed residuals 6 = - fi between the observations and some function f. We are primarily interested in estimating the noise for the purpose of "denoising" the measurements. Hence, we fit an elliptical distribution to the residuals by maximizing the (log) marginal likelihood with respect to the parameters nw: that
log p(e: nw) = Σ log I NN (EE: (, plws: n )dws.
(8)
For general mixing distributions, this intogral lacks a closod-form expression, but since it ono-dimensional we can approximate it efficiently by numerical integration (for example, using the trapezoidal rule).
Ultimately, arrive at the likelihood
N
p(ylf)= = II /
N (ys: fr.we) plwis n ddw.
(9)
5

---

Published in Transactions on Machine Learning Research (08/2023)
3.2 Prior
Recall that our main objective is to infer the latent function f" f(x") at arbilrary locations xt E RRP given a finite set of noisy observations y. In probabilistic machine learning, the mapping + f is often defined by the posterior predictive distribution
p( - 11000-MM110N M
(10)
which turns modeling into a search for suitable choices of Iff and p(fly). Accordingly, the noise estimation described in the previous section is only done in pursuit of this higher purpose.
Sparse formulation. For an elliptical process (EP) we can rewrite the posterior predictive distribution as
p(f -( ( u, fly))
(11)
where we marginalizing not only over the mixing variable & and the function values f (at the given inputs X), but also over the function values at the so-called M inducing inputs Z. Introducing inducing points lets us derive a sparse variational EP-a computationally scalable version of the EP similar to the sparse variational GP Titssas 2009
The need for approximation arises because of the intractable second factor p(f, u, Ely) Equation 11 (The first factor fff*(f()) is simply a Normal distribution.) We summarize the sparse variational EP below and refer to Appendices D and E for additional details.
Variational approximation. We make the variational ansatz p(f. (ly))) = p(flu, ()q(ulE) g(E), where u is the latent. function at the induced input locations Z, and parameterize this variational posterior as an elliptical distribution. We do so for two reasons: first, this makes the variational postcrior similar to the true posterior, and second, we can then use the conditional distribution to make predictions. In full detail, we factorize the posterior as
u. sp) p(flu, ( 42;);
(12)
where 4 (ep,, Pu; (PE) are the variational parameters, quulE; (u,) = N(m, St) is Gaussian distribution with variational parameters m and S, and the mixing distribution { ~~ glE; (P:). Again, q(; (pe) could be any positive finite distribution, but since we only know that the posterior is clliptical with a data-dependent mixing distribution, and since it is impossible to "overfit" a variational distribution (Bishop Nasrabadi 2006), we choose a very flexible yet tractable model, namely spline flow.
Note that, because of the conditioning pair in u. Thus, marginalization over results in a Gaussian distribution, for which the marginals of l only depend on the corresponding input x' (Salimbeni et al. 2019
the first two factors in Equation 12 are a Gaussian conjugate
= -N (0y,a""))))))))))))))
(13)
where
ux(x²) K.- Kulm.
(14)
af(2") kin. k K K=l k.,
(15)
un
and ki. "(("xx")) k. k< )) and K... k).
Predictions on unseen data points x= are then computed according to (sce Appendix E
a(s"ly:"") = Ey(e:ve)) (x""),)( ay(()))
(16)
For training, we use variational inference (VT), i.c., maximizing the evidence lower bound (ELBO) to indi- rectly maximize the marginal likelihood. We train the model using stochastic gradient descent and black-box
variational inference (Bingham al.] 2019 Wingate & Weber 2013 Ranganath al. 2014).

---

Published in Transactions on Machine Learning Research (08/2023)
VI training. The marginal likelihood is
ply: Mf. me) - f ply, f. u, nu, me ne) / p(yllfp(fluu ny ny)p(u, Tur me f duddg. (17)
This intogral intractable since p(E: my) is parameterized by a spline flow. To overcome this we approximate the marginal likelihood with the ELBO
CELBO Th, ne, Ppr Pu, PP) = (()))))))) [logp(y]f) - DRr. (a(u, 3s xp) p(u, E))
(18)
Σ Bll (lag pog DEL (g(u, E: pp.u
Had the likelihood been Gaussian, the expectation (e((ieepp) [log fi; my)] could have been computed analytically. In our case, however, it is elliptical, and we use a Monte Carlo estimate instead. Inserting the elliptical likelihood (Equation (9) from the previous section yiclds
N
C(n. 4) = B( [log (N (yis fi, plws, nw)) DKL (a(u, ()))
(19)
3.3 Extension to heteroscedastic noise
We now extend the elliptical likelihood to capture heteroscedastic (input-dependent) noise. The main idea is to let the parameters nJ of the likelihood's mixing distribution depend the input location z;-
In heteroscedastic regression, the noise depends on the input location æ;, For example, heteroscedastic elliptical noise can be useful in time series where the noise variance and tail-heaviness change over time. Examples of this can be found in statistical finance (Liu et al. 2020) and robotics [Kersting et al., 2007). To model this, we use a neural network 97 (+) with parameters Yo to represent the mapping from input location
to spline flow parameters, 9nw
nw:
We train the model by maximizing the log- ligeliill
((ww) = Σ log f pluilSa wi)p(ws nw: = 9r. (x:))
N
(20)
Additional information, such time of the day or season for time series data, can be incorporated by simply passing it as extra inputs to the neural network (-).
4 Experiments
We examined the variational elliptical processes using four different experiments. In the first experiment, we investigated how well the elliptical likelihood (Section 3.1 recovers known clliptical noise in synthetic data. In the second experiment, we demonstrated the use of the heleroscedasticc EP on synthetic dataset. We evaluated regression performance on seven standard benchmarks in the third experiment where we compared the sparse and the heteroscedastic EP formulations with both sparse GP (SVGP) Hensman et al. 2013) and full GP baselines. Finally, in the fourth experiment, examined if using an EP is beneficial in classification tasks.
The mixing distribution of the variational EP used a quadratic rational spline flow, where we transformed the likelihood flow p(w) using Softplus and the posterior flow p(E) using arctan to ensure that they were positive (remember that a mixing distribution must be positive, see Equation (3) We used a squared exponential kernel with independent length scales in all experiments. See Appendix for further implementation details. Accompanying code is found at
Implementation.
7

---

Published in Transactions on Machine Learning Research (08/2023)
95% crediblitt
Triinin dita
Learned pla)
90% credibility region
mean
True p(w)
20
-20
-15
-1.
-05
0.0
1.0
1.5
2.0
X
(a)
(b)
p(w)
2.0
10
0.5
0.0
1.0
1.5
2.0
0.8
x
(d)
4
p(c))
-20
-15
-10
-05
0.0
1.0 1.5
2.0
0.4
0.8
x
(e)
(f)
Figure The posterior predictive distribution when using GP with elliptical noise modeled by a spline flow. Fach row represents synthetic dataset with different. noise. The first adds Gaussian noise 4), and the third row adds 1). The shaded areas show the latent function posterior f" and the noisy posterior y" 90% credibility areas. The histograms show the learned and the true noise mixing distribution.
d(w 0.04) the second row adds Student's t noise ~Scale-Invxx Cauchy noise w ~Scale-Ime Live-In-x'" (vv-xxxxx111 Inv-x²(v
4.1 Noise identification
To examine how well the elliptical likelihood, described in Section 3.1 captures different lypes of elliptical noise, we created three synthetic datasets with the same latent function fi = sin(3x;)/2 sampled indepen- dently at N = 200 locations 24 ~ 2(- Further, each dataset was contaminated with independent elliptical noise e, that was added to the latent function, Y: = fi + 4,. The added noise varied for the three datasets. The first. was Gaussian distributed which is the The second was Student's t which means that follows the scaled inverse chi-squared distribution ~Saale- Inv-x²(v 4). The third was Cauchy distributed w ~Scale-Invxx GP for each dataset with an elliptical likelihood.
w being Delta. distributed d(w 0.04).
1). We trained sparse variational
Figure ! illustrates the results from the experiments. The histograms in the right column compare the learned mixing distribution p(w; m.) to the true mixing distribution (the red curve) from which the noise e, originated. The learned distributions follow the shape of the true mixing distribution reasonably well, considering the small number of samples, indicating that can learn the correct likclihood regardless of the noise variance. Furthermore, if the noise is actually Gaussian, likelihood. The left column shows the predictive posterior of the final models, demonstrating that the models managed to learn suitable kernel parameters jointly with the likelihood.
at -0.7, then so is the learned
8

---

Published in Transactions on Machine Learning Research (08/2023)
Robust regression on synthetic data. An elliptical likelihood is better at handling outliers and non- non-Gaussian noise than a Gaussian likelihood because it can better match the whole distribution of the noise rather than just a single variance. This is shown in Figure 5 where a GP with an clliptical and a Gaussian likelihood were trained on a small synthetic dataset with additive Student's t noise, with 7 4. The Gaussian likelihood approximates the mixing distribution with single variance at approx- imately w = 0.4, while the elliptical likelihood fits the entire mixing distribution. As a result, the GP with the Gaussian likelihood needs to use a shorter length scale to compensate for the thin tail of the likelihood. In con- trast, the GP with the elliptical likelihood can focus on the slower variations, thus producing better fit to data.
Tra fune Fred mak
2.5
Leamidd plao
plas
2.0
Density
1.5
1.0
0.5
0.0
1.0
(b)
3.0
Leamodd plus)
2.5
True
Density
2.0
1.5
1.0
0.5
4.2 Elliptic heteroskedastic noise
0.0
1.0
W
In this experiment, we aimed to exemplify the benefits of using heteroscedastic elliptical noise as described in Sec- tion [3.]] To this end, we created the synthetic dataset shown in Figure 0 It consisted of 150 samples gen- erated by adding heteroscedastie noise to the function (x;) sin(5x;;) + where 2; (0, 4). Specifically, we added Student's noise, c(x;) ~ St (v(zz), (zz)), where the noise scale followed v(x,) = 25 - 11 + 110.9 and the standard deviation by (xx) = 0.5|2z + +|| 1| + 0.001. We
(c)
(d)
Figure 5: The predictive posterior after training on small synthetic dataset. The shaded areas show the 95 % credibility area of the latent func- tion posterior I' (blue) and the noisy posterior y" (magenta) when using (top row) a GP with elliptical noise modeled by a spline flow and a (bottom row) a GP with Gaussian noise. The histograms show the learned and the true noise mixing distribution.
used a variational sparse GP with heteroscedastic noise as described in Section 3.3
The experimental results, depicted in Figure 0 show that, qualitatively, even though the rapid change in the noise distribution and the low number of training samples, the model captures the varying noise in terms of the scale and the increasing heaviness of the tail. Remember that a single spike in the mixing distribution, as at x = 0.7, indicates that the noise is Gaussian, and the wider the mixing distribution is, as at x = -0.7, the heavier-tailed the noise is. When the synthetic data has Gaussian noise, then so has also the learned elliptical noise. Similarly, when the synthetic noise is heavier-tailed, so the learned mixing distribution. This indicates that this model could be helpful for data with varying elliptical noise.
4.3 Regression
We conducted experiments on seven datasets from the UCI repository (Dua & Grafl 2017) to study the impact of the elliptical noise, elliptical posterior, and heteroscedastic noise. We used various regression models based on a Gaussian Process (GP) prior; see Table 1 for a summary. As baselines, we compare with the sparse variational GP model of Hensman et al. (2013). which we call SVGP, and an exact GP.
Models evaluated. We used GP model with elliptical noise (EP-GP) to compare its performance to the traditional GP model with Gaussian noise. Theoretically, an elliptical posterior should result from combining Gaussian prior and an elliptical likelihood, but in this case, we approximated the posterior with a Gaussian. We also included a model that used an elliptical posterior (EP - EP) to explore the potential benefits of using the theoretically more accurate elliptical posterior. Additionally, we tested a heteroscedastic elliptical noise model (Het-EP) and a heteroscedastic Gaussian noise model (Het-gP) to compare their performance. The difference between these two is that in Het-GP the neural network only predicts the noise variance, whereas the Het-EP model predicts the 26 parameters corresponding to nine bins of the spline flow.
9

---

Published in Transactions on Machine Learning Research (08/2023)
0.4
0.0
0.4
0.7
ploa
16
23
20
0.0
0.2
0.4
0.0
0.2
0.4
0.0
0.5
1.0
2.0
pla) synth
230
000
10
0.0
0.0
0.4
1.0
Figure 6: The results from training GP with heteroscedastic elliptical noise on a synthetic dataset. The top row shows the posterior distribution where the shaded arcas the 95 % credibility areas of the latent posterior l (magenta) and the noisy posterior y" (blue). The histograms in the middle row show the noise mixing distributions at the different z-values indicated by the vertical dashed lines in the top plot. The bottom row shows the mixing distribution used when creating the synthetic data.
MPG
Concrete
Elevateer
Bike 10427,
Calfomia
Kind0s 24000,
Frotein
9963,
Het EP
Het GP
EP-EP
EP-GP
0.4
09
0.56
0.40
0.45
-20
-1.
-10
9.6
08
-0.4
-02
0.95
1.90
1/55
Neg LL
MPG
Concrete
Ewalor nurrin 9560, 10
Bike
Califorie n_tail 12384,
KindOk r_tran 24000, d
Proten n 27438, d9
n_ " 294, d
rutain " 778.d=8
1,30h 10427, d=
EP-GP
SVGP
Exact GP
Figure Predictive negative log-likelihood (LL) (top) and mean squared error (MSE) (bottom) on held- out test data from the regression benchmarks (smaller is better). We show the average of the ten splits as a dot and the standard deviation of the error varies wetween the datasets. For example, the MSE error for the Bike dataset is low for all six models. Overall, the EP posterior outperforms the GP posterior, regarding the log-likelihood, for the five larger datasets
line. The models with bold fonts
models. Note that the spread
First, summarize the results in Figure 7 and then discuss the results of each method in more detail. The figure displays the mean and standard deviation of ten randomly chosen training, validation, and test data splits. The training procedure for all models optimizes directly or indirectly the log-likelihood. Therefore, the most relevant figure of merit is the negative test lig-likelihood (LL), shown in the top row of Figure 7 We stress that log-likelihood is more critical than MSE bocause mean squared error (MSE) does not consider the predictive variance. However, we show the MSE on held-out test sets in the bottom row for completeness. Additional details the experiments be found in Appendix G
10

---

Published in Transactions on Machine Learning Research (08/2023)
Table
The different types of models we trained on the regression datasets.
NAME
APPROXIMATION
LOSS
LIKELIHOOD POSTERIOR
Exact GP
Exact
Marginal likelihood
Gaussian
Gaussian
SVGP
Variational
ELBO
Gaussian
Gaussian
EP.GP
Variational
ELBO
Elliptic
Gaussian
EP-EP
Variational
ELBO
Elliptic
Elliptic
Hel-GP
Varialional
ELBO
Gaussian
Gaussian
Het-EP
Variational
ELBO
Elliptic
Gaussian
GP baseline. To assess the quality of the approximations introduced, we first establish an exact GP baseline that made predictions without any approximation. We trained the GP hyperparamelers using L-- BFGS and early stopping
validation dataset. For this to be feasible on large datasets, used the Blackbox Matrix-Matrix multiplication inference procodure (Gardner et al. 2018 Wang et al. 2019). In the following sections, we discuss each method in detail.
Variational GP approximation. First, compare the exact GP bascline with its variational approx- imation, i.e., SVGP. First, consider the results on MPG and Concrete, where SVGP did not make use of inducing points due to the small sample size. Consequently, SVGP is worse performance on Concrete is only due to the change of inference method. For the other datasets, we investigated the effect of the number of inducing points on the predictive log-likelihood; see Figure [ in Appendix H The dependence is very similar for all methods. In particular, the performance saturates at roughly 500 inducing points on all datasets except Kin40k, which continues to improve. However, the relative performance of the different methods Kin40k is fairly stable.
Elliptical likelihood. Next, we consider whether is advantageous to use an elliptical likelihood instead of Gaussian. To this end, we compare the performance of SVGP and EP-GP, which only differ in this respect. The results show that switching to an elliptical likelihood improves the log-likelihood on most datasets, would be expected theoretically.
Elliptical posterior. We now compare EP-GP to EP-EP analyze the potential benefit of having an elliptical posterior. On three of the datasets (Bike, Kin40k, and Protein), which are all relatively large, the elliptical posterior produces clear improvement in log-likelihood. In contrast, on the other datasets, it is similar, possibly because the posterior is well-approximated by Gaussian. Regardless, we conclude that, when using an elliptical likelihood, an elliptical posterior is preferable over Gaussian one.
Heteroscedastic models. Is there additional benefit of having heteroscodastic noise? On the two smallest datasets (MPG and Concrete), the answer is clearly no: the heteroscedastic models perform worse than SVGP and EP-GP in terms of both log-likelihood and mean squared error, indicating potential over- fitting and that regularization may be warranted. (Note that the most relevant comparisons are Het-GP vs. SVGP and Het-EP vs. EP-GP.)
On the remaining datasets, however, the heteroscedastic models clearly outperform SVGP and EP-GP terms of log-likeldhood. On the other hand, they perform poorly in terms of mean squared error; in fact, worse than SVGP on all datasets. Hypothetically, this is because the heteroscedastiec models attribute too much variation to the likelihood, thus sacrificing the mean-function prediction. This could potentially be mitigated by decoupling the mean and the covariance models (Salimbeni et al. 2018 Jankowiak et al. 2020 Another oot would be to increase the weight of the KL--divrgence term in the ELBO (Higgins et al. 2017 We expect such improvements to be more critical for the Het-EP model, which has a more flexible likelihood than Hct-GP. Still, Het-EP already performs slightly better than Het-gP on the three largest datasets. Note, however, the EP-EP model often achieves both a log- likelihoohd similar to the heteroscedastic models and mean squared error similar to SVGP.
11

---

Published in Transactions on Machine Learning Research (08/2023)
900
95% confidence region
Training data
Test data
confidence region Data points
Wind power production (GWh)
Predictive mean
2010
2011
2012
2013
2014
2015
2016
2017
Year
Figure 9: The predictive posterior after training the elliptical process on the wind power dataset. The shaded areas show the 95 % credibility area of the latent function posterior f" and the noisy posterior y".
Computational considerations. Empirically, we found that replacing the Gaussian likelihood with the elliptical likelihood had a minor impact on the computational demand. Further changing to an elliptical posterior increased the computational time per iteration, but a faster convergence partially compensated for this. Finally, modeling heteroscedastiec noise with a neural network adds significant complexity, but this was offset by running it on GPU.
Prediction accuracy. In summary, the results show that an elliptical likelihood results in better or equal predictive log- logllikels than Gaussian likelihood. However, the advantage is less significant small datasets. Similarly, the more flexible elliptical posterior tends to produce better results. However, when looking at the mean squared error (MSE), the exact GP outperforms the other models. Thus, if predictive performance is the main objective, an exact GP (or, even better, a neural network) may be the best choice. However, the well-known scalability issues of exact GP elearly limit its applicability. In such scenarios, our results suggest that EPEEP better choice than SVgP.
4.4 Application: Forecasting wind power production
Wind power stands for a significant and increasing share of global power pro- duction. However, since wind power generation inherently stochastic and hard to control it poses severe challenges to power system balance and stability (Impram al. 2020). In principle, these could be addressed via sophisticated control theoretical tools such as stochastic model predictive control (Schildbachi et al. 2014), but that requires accurate and reliable wind power forecasts. In this section, we illustrate the applicability of the elliptical process to time series forecasting. Specifically, we consider the task of making a probabilistic forecast of German cointry-wide totals of wind power production at the end of 2016 based on data from the preceding seven years. The data comes from Wiese
((w)
0.0
0.5
1.0
w
Figure 8: The mixing dis- tribution p(w) of the ellip- tical noise.
et al. 2019).
Since the data is strictly positive, model the time series in the log domain, using an elliptical process with elliptical noise (EP- EP). To capture the periodicity in the data use sum of periodic kernel and a linear kernel. See Appendix []] for more details.
Figure presents the predictive distribution from the elliptical process. The model successfully captures the underlying annual periodicity and slowly increasing trend while producing qualitatively reasonable credibility regions. Figure shows the mixing distribution of the trained elliptical likelihood, which is rather broad and hence non-Gaussian. The two modes discernible in the mixing distribution reflect the scasonal changes in the underlying data. For comparison, Appendix shows the corresponding results when using full GP and
12

---

Published in Transactions on Machine Learning Research (08/2023)
heteroscedastie EP, where the full GP scems to be more likely to overfit to short-scale trends in the data and produces credible regions that are too wide. The heteroscedastic EP, on the other hand, produces a fit that is overall on par with the EP model but disentangles the seasonal variations.
4.5 Binary classification
To evaluate the EP on classification tasks, we perform variational EP and GP classification by simply replacing the likelihood with a binary one. To derive the expectation in Equation (18) we first sample f; ~
Heart
Monmograzhic
EPonor, pPot
gi-ptor, Eppport
( and then derive the likelihood as Ber (f)SSigoii(()))..
AUC
Heart
Pina
Manmoggaphhe
This realization is interesting since do not have likelihood that captures the noise in the data; instead, the process itself has to do it. Therefore, we can indicate the value of the el- liptical process itself without the elliptical noise. We compare two variational EP models with a variational GP model. The two EPs differ in the prior mixing distributions, where the first model has GP prior and a EP posterior. For the sec- ond model, we replace the GP prior to an ellip- tical one. We can see the trainable prior mixing distribution as using a contimuously scaled mix- ture of Gaussian processes, which can be more expressive than single GP.
EPoroor, EP-post
(Piptor, EP-post
Accumcy
Heart
Pira
Nanmogaphic
EPiprir, EPooos
gaopior, EPopos
GP
-066
-0.0
-1.55
-132
-45
-0.0
-0.0
-0.5
-0.0
Logikklinood
The classification AUC (Area Under the Curve), accuracy, and predictive log-likelihood from the ten-fold rross-validation (higher is better). We show the average of the ten splits a dot and the standard devi- ation as a line.
Figure 10:
To evaluate the models, we performed a ten-fold cross-validation where we trained the models on three classification datasets, described in Appendix J Figure
10 presents the results from ten folds. From the area under the curve and the log-likelihood score, we see that the EP separates the two classes better, especially using the sparse models. The variational elliptical distribution is sufficient to get a high log- ligliiillhood while training the mixing distribution of the EP did not further improve the score.
5 Related work
In general, altempts at modeling heavy-tailed stochastic processes modify either the likelihood or the stochas- tie process
both. Approximate inference typically needed when going beyond Gaussian likelihoods (Neal 1997 Jylänki et al., 2011)), e.g., for robust regression, but approximations that preserve analytical tractability have been proposed (Shah et al.) 2014).
Elliptical processes gain flexibility by learning the mixing distribution, which makes training and inference reasonably efficient and the resulting predictions interpretable. It is, however, also possible to construet deep process priors by composing several stochastic process layers (Damianou & Lawrence 2013) and including non-linear transformations as activation functions (Aitchison et al.] (201).. The tractability of such deep processes depends on the specifics of the distributions involved. For instance, the deep inverse Wishart process (Aitchison et al., 2021) uses an inverse Wishart prior over the kernel just the Student's L process (Shah et al. 2014). This suggests it might be possible to generalize this approach to deep elliptical processes. While intriguing. we leave this to future work.
Other attempts at creating more expressive GP priors include Maroñas et. al. 2021). who used a GP in combination with a normalizing flow, and Luo & Sun (2017), who used a discrete mixture of Gaussian processes. Similar idcas combining mixtures and normalizing flows have also been proposed
create
more expressive likelihoods (Abdelhamed al. 2019 Daemi et al. 2019 Winkler et al. 2019 Rivero &
Dvorkin 2020 and variational posteriors (Nguyen & Bonilla 2014). Non-stationary extensions of Gaussian
13

---

Published in Transactions on Machine Learning Research (08/2023)
processes, such when modeling heteroscedastic noise, are somewhat rare. Examples include Zhao et al. (2021) who propose a hierarchical model in parameter space, the mixture model of Li et al. (0021) and the variational model of Lázaro-Gredilla & Titsias (2011).
Deep kernel learning (Calandra et al. 2016 Wilson et al. 2016b aj is another class of deep GPs that uses a neural network to learn the input features
GP. similar approach was taken by Ma al. (2019), who describe a class of stochastic processes where the fuitee-dimensional distributions are only defined implicitly parameterized transformation of some base distribution, thereby generalizing carlier work on warped Gaussian processes (Snelson et al. 2003 Rios & Tobar 2019). However, the price of this generality is that standard variational inference is no longer possible. Based on an assumption of Gaussian likelihood, they describe an alternative based on the wake-sleep algorithm by Hinton et al. (1955)
In the statistics literature, it is weoll that elliptical processes can be defined as scale mixtures of Gaus- sian processes (Huang Cambanis 1979 O'Hagan 1991 O'Hagan al.] 1999). However, unlike in machine learning, little emphasis is placed on building the models from data (i.e., training). These models have found applications in environmental statistics because of the field's inherent interest in modeling spatial extremes (Davison et al. 2012). Like us, several works take the mixing distribution as the starting point and make localized predictions of quantiles Jlaume-Deschamps et al. 2017) or other tail-risk measures (Opitz 2016
6 Conclusions
The Gaussian distribution is the default choice in statistical modeling for good reasons. Even so, far from everything is Gaussian Caussian-casuallyy pretending it is, comes at risk. The elliptical distribution offers a com- putationally tractable alternative that can capture heavy-tailed distributions. The same reasoning applics when comparing the Gaussian process to the elliptical process. sensible approach in many applications would be to start from the weaker assumptions of the elliptical process and let the data decide whether the evidence supports Gaussianity.
We constructed the elliptical process as scale mixture of Gaussian distributions. By parameterizing the mixing distribution using normalizing flow, we showed how corresponding elliptical process can be trained using variational inference. The variational approximation we propose enables us to capture heavy-tailed posteriors and makes it straightforward to create a sparse variational elliptical process that scales to large datasets.
We performed both experiments on regression and classification. In particular, we investigated the benefits of various combinations of elliptical posterior and elliptical likelihood and their heteroscedastic counterparts. We concluded that using elliptical likelihood and an elliptical posterior often achieves better log- likelihood and similar mean squared error as the sparse variational GP.
The added flexibility of elliptical processes could benefit a range of classical and new applications. However, advanced statistical models are not cure-all, and one needs to avoid overreliance on such models, especially in safety-critical applications.
7 Acknowledgments
We thank Sebastian Mair and Zheng Zhao for providing leedback on the manuscript. We also want thank the reviewers whose comments significantly improved the quality and clarity of our paper. This work was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Kmut and Alice Wallenberg Foundation, the Swedish Foundation for Strategic Research grant SM19-0029, and by Kjell och Märta Beijer Foundation.
14

---

Published in Transactions on Machine Learning Research (08/2023)
References
Abdelrahman Abdelhamed, Marcus A Brubaker, and Michael S Brown. Noise flow: Noise modeling with conditional normalizing flows. In Proceedings of the IEEE/CVP International Conference on Computer Vision, pp. 3165-3173, 2019.
Laurence Aitchison, Adam Yang, and Sebastian W Ober. Deep kernel processes. In International Conference on Mochine Learning (ICML), pp. 130-140. PMLR, 2021.
Jesús Alealé-Fdez, Alberto Fernández, Julián Luengo, Joaquín Derrac, Salvador Garcia, Luciano Sánchez, and Francisco Herrera. Keel data-mining software tool: data set repository, integration of algorithms and experimental analysis framework. Journal of Multiple- Valued Logic & Soft Computing, 17, 2011.
Maria Bånkestad, Jens Sjölund, Jalil Taghia, and Thomas Schön. The elliptical processes: family of fat-tailed stochastic processes. arXiv preprint arXiv:203..0200, 2020.
Rabindra Nath Bhattacharya and Edward Waymire. A basic course in probability theory, volume 69.
Springer, 2007.
Eli Bingham, Jonathan P Chen, Martin Jankowiak, Fritz Obermeyer, Necraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D Goodman. Pyro: Deep universal probabilistic programming. Journal of Machine Learning Research, 20(1):973-99) 978, 2019.
Nicholas H Bingham and Rüdiger Kiesel. Semi-parametric modelling in finance: thcoretical foundations.
Quantitative Finance, 2(4)2241, 2002.
Christopher M Bishop and Nasser M Nasrabadi. Pattern recognition and machine learning, volume 4. Springer, 2006.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. Journal of the Americon statistical Association, 112(518):859-877 877, 2017.
Roberto Calandra, Jan Peters, Carl Edward Rasmussen, and Marc Peter Deisenroth. Manifold Gaussian processes for regression. In 2016 International joint conference neural networks (IJCNN), pp. 3338 3345. IEEE, 2016.
Atefeh Daemi, Hariprasad Kodamana, and Biao Huang. Gaussian process modelling with Gaussian mixture likelihood. Journal of Process Control, 81:209- 220, 2019.
Andreas Damianou and Neil D Lawrence. Deep Gaussian processes. In International Conference on Artificial Intelligence and Statistics (AISTATS), pp. 207-215, 2013.
Anthony Davison, Simone A Padoan, Mathicu Ribatet, et al. Statistical modeling of spatial extremos. Statislical science, 27(2):161-188 186, 2012.
Peter J Diggle, Jonathan A Tawn, and Rana A Moyeed. Model-based geostatistics. Journal of the Royal Statistical Sociely: Series C (Applied Statistios), 47 (3)::299-33))) 1998.
Hadi Mohaghegh Dolatabadi, Sarah Erfani, and Christopher Leckie. Invertible generative modeling using linear rational splines. In International Conference on Artificial Intelligence and Statistics (ATSTATS),
pp. 4236 4246, 2020.
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http:///rrchive.ic. icc. uci. edu/ m1
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. In H. Wal- lach, H. Larochelle, A. Beygelzimer, F. d' d'Alhé--Bue, E. Fox, and R. Garnett (eds.), Advances in Neural Information. Processing Systems (NeurIPS), volume 32, 2019.
Hadi Fanace-T T and Joao Gama. Event labeling combining ensemble detectors and background knowledge. Progress in Artificial Intelligence, ((2):113-122) 2014.
15

---

Published in Transactions on Machine Learning Research (08/2023)
Kai-Ta Fang, Samuel Kotz, and Kai Wang Ng. Symmetric multivariate and related distributions. Chapman and Hall, 1990.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. Gpytorch: Black- box matrix-matrix gaussian process inference with gpu acceleration. In Advances in Neural Information Processing Systems (NeurIPS), volume 31, 2018.
James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data, In Ann Nicholson and Padhraic Smyth (ods.), Uncertainty in Artificial Intelligence (UAT), volume 29. AUAI Press, 2013.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification. In International Conference on Artificial Intelligence and Statistics (AISTATS), pp. 351 360, 2015.
Irina Higgins, Loic Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick, Shakir Mohamed, and Alexander Lerchner. 3-VAE: Learning basic visual concepts with a constrained variational framework. In International conference learning representations (TCLR), 2017.
Geoffrey E Hinton, Peter Dayan, Brendan J Frey, and Radford M Neal. The wake-sleep algorithm for unsupervised neural networks. Science, 268(5214):1158 1161, 1995.
Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley. Stochastic variational inference. Journal of Machine Learning Research, 2013.
Steel T Huang and Stamatis Cambanis. Spherically invariant processes: Their nonlinear structure, discrim- ination, and ostimation. Journal of Multivariate Analysis, 911):59-83, 1979.
Semich Impram, Secil Varbak Nese, and Bülent Oral. Challenges of renewable energy penetration on power system flexibility: A survey. Energy Strategy Reviews, 31:100539, 2020.
Martin Jankowiak, Geoff Pleiss, and Jacob Gardner. Parametric Gaussian process regressors. In International Conference Machine Learning (ICML), pp. 4702 4712. PMLR, 2020.
Pasi Jylänki, Jarno Vanhatalo, and Aki Vehtari. Robust Gaussian process regression with Student-t likelihood. Journal of Machine Learning Research, 12(Nov):3227-3255 2011.
Yutaka Kano. Consistency property of elliptic probability density functions. Journal of Multivariate Analysis, 51(1):139-144 147, 1994.
Kristian Kersting, Christian Plagemann, Patrick Pfaff, and Wolfram Burgard. Most likely heteroscedastic Gaussian process regression. International Conference on Machine Learning (ICML), pp. 393 400, 2007.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Confer- ence on Learning Representations (ICLR), 2015.
Diederik P. Kingma and Max Welling. Auto-encoding variational bayes. In International Conference on Learning Representations (ICLR), 2014.
Miguel Távaro-Gredilla and Michalis K Titsias. Variational heteroscedastic Gaussian process regression. Tn International Conference on Machine Learning (ICML), 2011.
Tao Li, Di Wu, and Jinwen Ma. Mixture of robust Gaussian processes and its hard-cut EM algorithm with variational bounding approxination. Neurocomputing, 452:224 238, 2021.
Bingqing Liu, Ivan Kiskin, and Stephen Roberts. An overview of Gaussian process regression for volatility forecasting. In International Conference on Artificial Intelligence in Information and Communication (ICAIIC), pp. 681-686, 2020.
Chen Luo and Shiliang Sun. Variational mixtures of Gaussian processes for classification. In International Joint Conferences Artificial Intelligence (IJCAT), volume 357, pp. 4603 4603-4609, 4609, 2017.
16

---

Published in Transactions on Machine Learning Research (08/2023)
Chao Ma, Yingzhen Li, and José Miguel Hernández-Lobato Variational implicit processes. In International Conference on Machine Learning (TCML), pp. 1222-4233, 2019.
Benoit. Mandelbrot. The variation of certain speculative prices. The Journal of Business, 36(4):394 419, 1963.
Juan Maroñas, Oliver Hamelijnck, Jeremias Knoblauch, and Theodoros Damoulas. Transforming gaussian processes with normalizing flows. In Arindam Bancrjec and Kenji Fukumizu (eds.), International Confer- ence on Artificial Intelligence and Statistics (AISTATS), 2021.
Véronique Maume- Deschamps, Didier Rullière, and Antoine Usseglio- Carleve. Quantile predictions for ellip- tical random fields. Journal of Mullivariale Analysis, 159:1-17, 2017.
Radford M Neal. Monte Carlo implementation of Gaussian process models for Bayesian regression and classification. Technical Report 9702, Department of Statistics, University of Torouto, 1997.
Trung V Nguyen and Edwin V Bonilla. Automated variational inference for Gaussian process models. Advances in Newral Information Processing Systems (NeurIPS), 27, 2014.
Scbastian W Ober and Laurence Aitchison. Global inducing point variational posteriors for Bayesian neural networks and deep Gaussian processes. In International Conference Machine Learning (ICML), pp. 8248 8259, 2021.
Anthony O'Hagan. Bayes Hermile quadrature. Journal of statistical planning and inference, 29(3):245 260, 1991.
Anthony O'Hagan, Mare C Kennedy, and Jeremy E Oakley. Uncertainty analysis and other inference tools for complex computer codes. In Bayesian statistics 6, pp. 503-524. Oxford University Press, 1999.
Thomas Opitz. Modeling asymptotically independent spatial extremes based Laplace random fields.
Spatial Statistics, 16:1 18, 2016.
R. Kelley Pace and Ronald Barry. Sparse spatial autoregressions. Statistics & Probability Letters, 33(3): 291 297, 1997.
George Papamakarios, Eric Nalisuick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshmi- narayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1-64 2021.
Rajesh Ranganath, Sean Gerrish, and David Blei. Black Box Variational Inference. In International Con- ference on Artificial Intelligence and Statistics (AISTATS), pp. 814 822, 2014.
Carl Edward Rasmussen and Christopher K 1 Williams. Gaussian processes for machine learning. The MIT Press, 2006. p. 194.
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International Con- ference on Machine Learning (ICML), pp. 1530-1538, 2015.
Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic backpropagation and approximate inference in deep generative models. Tn International Conference on Machine Learning (ICML), pp. 1278- 1286, 2014.
Gonzalo Rios and Felipe Tobar. Compositionally-warpee Gaussian processes. Neural Networks, 118:235 246, 2019.
Ana Diaz Rivero and Cora Dvorkin. Flow-bascd likelihoods for non-Gaussian inference. Physical Review D, 102(10):003507, 2020.
Hugh Salimbeni, Ching-An Cheng, Byron Boots, and Mare Deisenroth. Orthogonally decoupled variational Gaussian processes. Advances in neuraa information processing systems (NeurIPS), 31, 2018.
17

---

Published in Transactions on Machine Learning Research (08/2023)
Hugh Salimbeni, Vincent Dutordoir, James Hensman, and Mare Deisenroth. Deep Gaussian processes with importance-wer weighted variational inference. In International Conference on Machine Learning (TCML., pp. 5589 5598, 2019.
Georg Schildbach, Lorenzo Fagiano, Christoph Frei, and Manfred Morari. The scenario approach for stochas- tic model prodictive control with bounds on closod-loop constraint violations. Automatica, 50(12):3009- 3018, 2014.
Amar Shah, Andrew Wilson, and Zoubin Ghahramani. Student-t processes as alternatives to Gaussian processes. In International Conference on Artificial Intelligence and Stalislics (AISTATS), pp. 877-885, 2014.
Jack W Smith, James Everhart, WC Dickson, William Knowler, and Robert Scott Johannes. Using the adap learning algorithm to forecast the onset of diabetes mellitus. In Proceedings of the annual symposium on compuler application in medical care, pp. 261. American Medical Informatics Association, 1988.
Edward Snelson, Zoubin Ghahramani, and Carl Rasmussen. Warped gaussian processes. In Adoances in
Newral Information Processing Systems (NeurIPS), volume 16, 2003.
Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In International Conference on Artificial Intelligence and Statistics (AISTATS), pp. 567-574, 2009.
Dustin Tran, Rajesh Ranganath, and David M Blei. Variational Gaussian process. In International Confer- ence on Learning Representations (ICLR), 2016.
Martin J Wainwright, Michael I Jordan, et al. Graphical models, exponential families, and variational inference. Foundations and Trends in Machine Learning, 1(1 2):1 305, 2008
Ke Wang, Geoff Pleiss, Jacob Gardner, Stephen Tyree, Kilian Q Weinberger, and Andrew Gordon Wilson. Exact gaussian processes on a million data points. In Advances in Neural Information Processing Systems (NeurIPS), volume 32, 2019.
Frauke Wiese, Ingmar Schlecht, Wolf-Dieter Bunke, Clemens Gerbaulet, Lion Hirth, Martin Jahn, Friedrich Kunz, Casimir Torenz, Jonathan Mühlenpfordt, Juliane Reimann, et al. Open power system data- frictionless data for electricity system modelling. Applied Eneryy, 236:401 409, 2019.
Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing. Deep kernel learning. Tn International Conference on Artificial Intelligence and Statistics (AISTATS), pp. 370-378. PMLR, 2016a.
Andrew Gordon Wilson Zhiting Hu, Ruslan Salakhutdinov, and Eric P. Xing. Deep kernel learning. In International Conference on Artificial Intelligence and Statistics (AISTATS), 2016b.
David Wingate and Thcophanc Weber. Automated variational inference in probabilistic programming. arXiv preprinl ar Xiv: 1301. 1299, 2013.
Christina Winkler, Daniel Worrall, Emiel Hoogeboom, and Max Welling. Learning likelihoods with condi- tional normalixing flows. arXie preprint or
2019.
I-C Yeh. Modeling of strength of high--peformancc concrete using artificial neural networks. Cement and Concrete research, 28(12):1797-188) 1808, 1998.
Zheng Zhao, Muhammad Emzir, and Simo Särkkä. Deep state-space Gaussian processes. Statistics and Computing, 31:1-26, 2021.
Abdelhak M Zoubir, Visa Koivunen, Yacine Chakhchoukh, and Michael Muma. Robust estimation in signal processing: A tutorial-stylc treatment of fundamental concopts. IEEE Signal Processing Magazine, 29(4): 61 80, 2012.
18

---

Published in Transactions on Machine Learning Research (08/2023)
A The elliptical distribution
The Gaussian distribution-thh basic building block of Gaussian processes- several attractive properties that we wish the elliptical process to inherit, namely (i) closure under marginalization, (ii) closure under conditioning, and (iii) straightforward sampling. This leads us to consider the family of consistent elliptical distributions. Following Kano (1994), we say that a family elliptical distributions {p(u(yn); m) | E N} is consistent if and only if
p(u(yN+1): n) dyN-1 (u(yN): 7))
(21)
In other words, a consistent clliptical distribution closed under marginalization.
Far from all elliptical distributions consistent, but the complete characterization of those that is provided by the following theorem (Kano 1994
Theorem An elliptical distribution is consistent if and only if it originates from the integral
p(u; n) = 12
E E pl:: nz)dE,
(22)
E
where & is a mixing variable with the corresponding, strictly positive finite, mixing distribution p(E: n), that is independent of N.
This shows that consistent elliptical distributions p(u; 17) are ccale-mixtures of Gaussian distributions, with mixing variable p(E: 7).. Note that any mixing distribution fulfilling Theorem can be used to define a consistent elliptical process. We recover the Gaussian distribution if the mixing distribution is Dirac delta function and the Student's distribution if it is scaled inverse chi-square distribution.
If p(u; 7) is a scale-mixture of normal distributions, it has the stochastic representation
Y| ~ (M, ΣE), E~ <<: n) By using the following representation of the clliptical distribution,
(23)
= 1 1/2 1///2/1/2
(24)
where Z follows the standard normal distribution, get the mean
E|Y] = 1/ |z] E μ
(25)
and the covariance
Cov(Y) = (Y) - =
=E 1E (E)E
= =E ( EE(
zz (x113)
= E
(26)
The variance is a scale factor of the scale matrix Σ. To get the variance we have to derive [g]. Note that if follows the scaled inverse chi-square distribution, E[g] = v/(v - 2). We recognize this from the Student's distribution, where Cov(Y) v/(v 2 2)E.
B Conditional distribution
To use the EP for predictions, we need the conditional mean and covariance of the corresponding elliptical distribution, which we derive next. We partition the data as y [y1, where y1 are the N₂ observed data points, y2 the next N₂ data points to predict, and N₁ No N. We have the following result:
19

---

Published in Transactions on Machine Learning Research (08/2023)
Proposition 1 If the data y 1y1. y2l originate from the consistent elliptical distribution in Equation the conditional distribulion originates from the distribulion.
3
Pyzlun (32)
CN
E se (u2) (0g)) p(E: v)d
(27)
Eeat1
(2x)
with the conditional mean E E 12l1 and the conditional covariance
Cov|Y2|Y, y2l E2() & E~ Ely
(28)
where = (y) -1) constant. The conditional scale matriz Σ2nt1 and the conditional meam veclor Paj1 are the same as the mean and the covariance matriz for Gaussian distribution.
- En pI), W2/1 = (y2 - (2l1 2 (32 (2n1), and CNim is normalization
The conditional distribution is guaranteed to be consistent elliptical distribution but not necessarily the same as the original one-the shape depending on the training samples. (Recall that consistency only concerns the marginal distribution.)
Proof of proposition 1
The joint distribution of y,, y2) is p(y (2( 7) and the conditional distri- bution of y2s given yı is (lgelli.E)p(Elbl {)p(Ely1: 77).
For given E, () is the conditional normal distribution and so
p(uzly1. E) ~ N( E22pE), E~ p(Elp): n)
(29)
where,
M2j1 = p12 + E21 Eji (y) (y) - 11)
(30)
Σ221 = Σ22 Σ21 Eji ET Σ21,
(31)
the same as for the conditional Gaussian distribution. We obtain the conditional distribution p(Elu1 m) by remembering that.
p(v1k)) ~ NM(1, E)1))
(32)
Using Bayes' Theorem get
p(ll: n) p(y)15))(E: n)
|nn 1/2 exp
u1 p(E: n)
2 f
x N212 exp
{ -5 U1
p(E: n))
(33)
Recall that u1 (y - -11) E 11)). We normalize the distribution by
N l0 5- E-N//2 exp
{
U1 p(e: n)d
(34)
2
The conditional mixing distribution is
p(C yı; 77) oNi.ns
N2.2 exp
p(E: n7).
(35)
2
The conditional distribution of y2 given is derived by using the consistency formula
p(y2lv1) = Ex1 (12 No/2 N2/2 g g-Na/2 exp p(Elo))dE
U2j1
(36)
2
where U2j1 = (y2 )2)1a Exg (y2 - (211). Using Equation 35
we get
p(y2111)
CN.
e- n/2 (wap (wpp ( p(E: n)dg.
(37)
1/2 (2 No12
20

---

Published in Transactions on Machine Learning Research (08/2023)
C Derivation of the credible intervals of the elliptical process
We derive the credible interval of the elliptical process, by using the Monte Carlo approximation of the integral, as
1
{ -/122 =2/(200)) ( p(E)dEdx
p < < 20)
(38)
OV 2T
1
20
1
m
1/2 =x =/2(00)))) dx
(39)
0 V
2m
M
(=1
50
=f =22(02)0)) dx
(40)
omV
2n
1=1
V2Ei
e
du
(41)
MV
π
erf
(42)
m
22
For every mixing distribution, can derive the credibility of the prediction. It is the number of samples we take that decides the accuracy of the credible interval.
D Details on the non-sparse variational elliptical process
For a Gaussian process, the posterior of the latent variables f is p(fl) x plylf)((f)
(43)
Here, the prior x(f ~ No0 K) is a Gaussian process with kernel K and the likelihood p(ylf) ~ ~NI, I2) is Gaussian. The postcrior derives to
p(flv) ~ N
y. (K to to 2))
(44)
and the predictive distribution of an arbitrary input location x² is
p( f MFIFF((W)
(45)
where pf, x²) the conditional distribution, which again Gaussian with
N (f
kk, (ke - k. [K + o2
(46)
Going back to the elliptical process, we want to derive the predictive distribution. The problem, though, is that the posterior now intractable. In order to get tractable posterior, we train the model using variational inference, where we approximate the intractable posterior with a tractable one,
p(f. Ely: ny, ne) g(f, E: &p = a(fls; (P))((; spg).
(47)
Here of are the parameters of the prior GP process (such as the kernel parameters). ns are the parameters of the mixing distribution, and q(flE: ppy) N(m,, Sy where my and Sy are variational parameters. The postcrior 45; (p&) is parameterized with any positive distribution such normalizing flow. We use this approximation when we derive the predictive distribution
p(f*))) = / p( mr)p(f.el) ns de)
(48)
/ piS'Sss; mf)e(f Elu: nr n dffd
(49)
- ns)a(fls: sp))((E: se)ddf
(50)
21

---

Published in Transactions on Machine Learning Research (08/2023)
By first taking look at the prior distribution p(f f(5) when is constant,
f- & N (0,
Kin k__
k. K {)
(51)
we arrive at the conditional distribution
pff &: ns) = KK (Riz. - 1.K
(52)
We use this expression together with the variational approximation to derive the posterior predictive distri- bution
/ pPSMIS ns)a(fle sp))(E;; speddfde
p(f'lz))
(53)
I purMIss:m vy)df
(54)
/ f. (Kes - k. K KE_KK'K, kk)))) N (E)
= =BulE: vP&)
(55)
(56)
where
My(x²) = k2_ K--K-m 0f(x²) hina k. (K -
(57)
K-SSK" 1) k..
(58)
We get the variance by [x())
Optimizing the ELBO We train the model by optimizing the evidence lower bound (ELBO) given by
(ELBO ne; s P5) Bychsss)) [logrlulf) - DKL (a(S, 9p, (2) 1lp)) ng))
(59)
E
Details on sparse elliptical processes
With the variational inference framework, we create a sparse version of the model
f plf u. m)ds = / p(flu,,; m))p(uls; nw)p(E: ng)dg,
(60)
where u are outputs of the elliptical process located at the inducing inputs Z. We approximate the posterior with
pff. u. Ely: 7) p(f u. E: ny ny)q(ul;; pu))(E: sp²). The posterior predictive distribution then given by
(61)
p(f" 8) f p( u, wp(f u, Ely: dddfdd
f plsissu u. m)(( E; ng)g (u)g(E; d Rd
/ / / u m)((su, &: ny)df g(ulf; (P))(EE spd)dudg.
(62)
We can simplify the inner expression by using the fact that the elliptical distribution is consistent,
f [Pf u, m)p(fu. n)df / flu, n)df pu n)
(63)
22

---

Published in Transactions on Machine Learning Research (08/2023)
Hence, Equation
62 is simplifies to
p( - f (l m)a(uls pu)a(Ei Pe)duds,
(64)
where q(ul:: Pu) (m.. S,4) with the variational parameters my and Su, and & is parameterized, e.g., by a normalizing flow.
Finally, we obtain the posterior p( = (eterve) (M(M))gga"")) )))) where
pg(x²) = -k_ KuKK Kuum
(65)
,f(x²) hin. - k. (Kuu
K Kun Skuk) k..
(66)
Here, k, k(x', Z), k.+ k( x²), and Kuu k(Z, Z).
F Implementation: variational inference
We used the Pyro library (Bingham et al. 2019 universal probabilistic programming language (PPL) written in Python and supported by PyTorch on the backend.
In Pyro, we trained a model with variational inference (Kingma & Welling 2014 by creating "stochas- tic functions' called model and a guide, where the model samples from the prior latent distributions p(f, w, v), and the observed distribution p(ylf,w), and the guide samples the approximate posterior qf E: (p))((E; pg). We then trained the model by maximizing the evidence lower bound (ELBO), where simultancously optimized the model parameters n and the variational parameters s. (See more details
here
To implement the model in Pyro, we created the guide and the model (see Algorithm 1) by building upon the already implemented variational Gaussian process. We used the guide and the model to derive the ELBO. which we then optimized with stochastic gradient descent using the Adam optimizer (Kingma & Ba, 2015
Algorithm Py PyTorch implementation of the variational sparse elliptical process (VI-EP-EP).
1: procedure MODEL(X, y)
2: Sample { from p(E: ne)( Normalizing flow) Sample from NN )u))
Take sample from the latent and II -N b During training is sampled from the 1=1
4: Derive the variational posterior
posterior/guide
5: Take Monte-Carke sample f. from each ((filE:p))
For or each yi approximate by, = log JN (yilfi, p(u; m mow using the trapezoid rule.
7: Get the log probability of y by Ei Eye-
end procedure
9: procedure GUIDE
10: Sample { from ((::4p)) Normalizing flow
11: Sample from ((m,
12: end procedure
G Regression experiment setup
In the regression experiments in Section 4.3 we ran all experiments using the Adam optimizer (Knggaa & Ba 2015 with learning rate of 0.01. For the full GP, we used the L-BFGS optimizer to train the hyperparameters. Here, we, in the same way as for the other models, used early stopping on a validation dataset, which operated by saving the model with the lowest validation log predictive likelihood.
For all experiments, we created ten random train/val/test splits with the proportions 0.6/0.2/0.2, except for the two smallest. datasets (mpg and concrete), where we neglected the validation dataset and used
23

---

Published in Transactions on Machine Learning Research (08/2023)
train/test proportions of 0.7/0.3. For the test set evaluation, used the model with the highest predictive probability on the validation set. For the large datasets (N > 1000), we used 500 inducing points. We did not use a sparse version of the model for the small datasets but instead used Z = Xtrain and kept them fixed during the training. We run the optimizer for the large dataset for 500 epochs and the small dataset for 2000 epochs.
Elliptical process setup. The likelihood mixing distribution uses spline flow with nine bins and Softplus its output transformation. The elliptic posterior mixing distribution Sigmoid output transformation. These parameters were not tuned, and fixed for all experiments.
spline flow with five bins and
For the beteroscedastic noise models, we used two-layer neural networks with hidden dimensions of 128. For the elliptical noise, we learned spline flow with nine bins which results in 26 hyperparameters to learn while for the heteroscodastic Gaussian likelihood we learned the variance solely.
H Results
The regression results from Figure 7 are presented in Tables 2 and 3 Figure 11 presents the outcome for different numbers of inducing points. We see that the results have stabilized for almost all datasets at 500 inducing points. We also notice that the relative log-likelihood between the models stays constant after 400-500 inducing points.
Table 2: Predictive mean squared error (MSE) on the hold-out sets from the experiments. We show the average of the ten random splits and one standard deviation in parenthesis.
MPG
CONCRETE
ELEVATORSS
BIKE
CALIFORNIA
KIN40K
PROTEIN
Het-E'P
0,144 (0.018)
0.127 (0.014)
0.149 (0.005)
0.018 (0.002)
0.223 (0.012)
0.141 (0.008)
0.531 (0.009)
Het-GP
0.142 (0.020)
0.178 (0.022)
0.148 (0.005)
0.020 (0.002)
0.230 (0.011)
0.112 (0.007)
0.508 (0.007)
EP-EP
(0.027)
(0.022)
(0.001)
(0.012)
(0.002)
(0.008)
EP-97
0.121 (0.026)
0.128 (0.013)
0.145 (0.005)
0.011 (0.001)
0.226 (0.013)
0.056 (0.003)
0.481 (0.007)
SVGP
(0.018)
(0.011)
(0.004)
(0.001)
(0.012)
(0.001)
(0.007)
Exact GP
0.135 (35)
0.103 103)
0.134 (0.134)
0.002 (0.002)
0.134 (0.134)
0.006 (0.000)
0.357 (0.006)
Table 3: Negative log likelihood (Neg LL) on the hold-out test sets from the experiments. We show the average of the ten random splits and one standard deviation in parenthesis.
MPG
CONCRETE
ELEVATORS
BIKE
CALIFORNIA
KIN40K
PROTEIN
Het-EP
0.463 (0.089)
(0.033)
0.400 (0.018)
-1.327 (0.020)
(0.022)
-0.397 (0.012)
(0.017)
Het-GP
0.530 (0.155)
0.456 (0.074)
0.399 (0.017)
-1.532 (0.047)
0.376 (0.021)
-0.316 (0.010)
0.986 (0.013)
EP-EP
0.266 (0.074)
0.443 (0.083)
0.425 (0.015)
-1.406 (0.028)
0.506 (0.022)
-0.246 (0.020)
0.976 (0.008)
EP-GP
0.268 (0.071)
0.344 (0.031)
0.427 (0.014)
-1.304 (0.023)
0.515 (0.021)
-0.053 (0.049)
1.056 (0.006)
SVGP
0.352 (0.062)
0.382 (0.030)
0.116 (0.013)
-0.865 (0.016)
0.619 (0.022)
-0.028 (0.004)
1.056 (0.006)
Exact GP
(0.387)
(0.206)
(0.463)
(-1.103)
(0.463)
(0.097)
(0.005)
24

---

Published in Transactions on Machine Learning Research (08/2023)
Train, Elevators
Val, Elevators
Train, Bike
Val, Bike
x
x
-0.8
x
x
SVGP
0.44
X
-1.0
x
x
x
x
X
X
x
EPGP
EPEP
042
-1.2
Gr-hetero
ETheeter
Neg
Neg
1.4
0.40
1.6
1.8
2.0
200
400
600
200
400
800
200
400
600
200
400
600
inducing
inducing
inducing
inducing
Train, Calefomia
Val. Calefornia
Train, Kin40k
Val,, Kir40k
0.85
0.2
X
0.1
Neg 11
0.55
-0.1
Neg
0.50
-0.2
0.45
-0.4
0.40
200
800
200
400
200
400
600
Num inducing
Num inducing
Num inducing
Num inducing
Train, Protein
Val, Prolein
X
x
1.10
Neg LL
9
X
1.05
0.85
200
400
600
200
400
800
Num inducing
inducing
Figure 11: The train and validation negative log- log-likelihood for the datasets using a varying number of inducing points.
Application: Wind power production
In the wind power production experiment, we used an EP process with elliptical likelihood, elliptical poste- rior, and 300 inducing points; heteroscedastic EP with heteroscedastic elliptical likelihood, and elliptical posterior; and a full GP with Gaussian likelihood. They all used kernel that was sum of periodic kernel and a linear kernel. Before training, we transformed to data to log scale and then normalized it. See Figure 12 for a visualization ion of the data before and after the transformation. This transformation made the data more symmetric around its mean and removed the nonnegative constraints.
Figure Section 4.4 shows the results when training the data with EP process with elliptical likelihood. For comparison, here plot the results when using full GP (Figure 13 and heteroscedastic EP (Figure 14
25

---

Published in Transactions on Machine Learning Research (08/2023)
Before transformation
After transformation
800
Wind power
600
400
200
-2
2010
2012
2014
2016
2018
2010
2012
2014
2016
2018
Year
Year
Figure 12: The data used in the wind power production experiments, where the left plot illustrates the data before is transformed and the right plot shows the data after it transformed, which is what we used as input to the model.
The main takeaway is that all three models able capture the periodicity in the data. The full GP seems to have lendency to overfit the data and also to create wide credibility regions. This probably because the data is non-Gaussian and the thin tails of the GP forces variances to be high. The heteroscedastic EP fits the data well, although the difference from using a non-heteroscedastic likelihood is small. Figure 15 illustrates the mixing distribution during some of the months of the year where we clearly can see that we have a seasonal change in the noise.
900
95% confidence region
Training data
Test data
Wind power production (GWh)
800
y" 95% confidence region
700
Predictive mean
600
500
400
300
200
100
2012
2013
2014
2015
Year
Figure 13: Wind power forecast using an exact GP The plot shows the predictive posterior where the shaded areas show the 95 % credibility area of the latent function posterior f and the noisy posterior y*:
J Datasets
Bike dataset Fanaee-T & Gama 2014) is obtained from bike sharing data, especially it contains the hourly and daily count of rental bikes between the years 2011 and 2012 with the corresponding weather and seasonal information.
Elevators dataset Duua Grail 2017 is obtained from the task of controlling an F16 aircraft, and the objective related to an action taken on the elevators of the aircraft according the status attributes of the airplane.
Physicochemical properties of protein tertiary structure dataset The dataset is taken from CASP 5-9. There 45730 decoys and sizes varying from 0 to 21 Armstrong.
26

---

Published in Transactions on Machine Learning Research (08/2023)
900
95% confidence region
Training
confidence region Data points
Wind power production (GWh)
Predictive mean
2010
2011
2012
2013
2014
2015
2016
2017
Year
Figure 14: The predictive posterior after training the elliptical process with heteroscedastic noise the wind power dataset. The shaded arcas show the 95 % credibility arca of the latent function posterior f* and the noisy posterior y"
Mars
June
September
December
20
10
0.7
0.5
0.5
0.6
0.7
0.6
0.7
W
W
(
The mixing distribution of the elliptical heteroscedastic likelihood sampled during different months during the year.
Figure 15:
California housing dataset was originally published by Pace & Barry 1997). There are 20 640 samples and 9 feature variables in this dataset. The targets are the prices of houses in the California area.
The Concrete dataset Yeh. 1998) has 8 input variables and 1030 observations. The target variables are the concrete compressive strength.
Auto MPG dataset (Aleal--Fdez al. 2011 originally from the StatLib library which is maintained Carnegie Mellon University. The data concerns city-cycle fuel consumption in miles per gallon and consists of 392 samples with five features each.
Pima Indians Diabetes Database Smith et al. 1988) originally from the National Institute of Diabetes and Digestive and Kidney Diseases. The objective of the dataset is to predict diagnostically whether or not a patient has diabetes based on certain diagnostic measurements included in the dataset. The dataset consists of 768 samples with eight attributes.
The Cleveland Heart Disease dataset consists of 13 input variables and 270 samples. The target classifies whether a person is suffering from heart disease not.
The Mammography Mass dataset predicts the severity (benign malignant) of mammography mass lesion from BI-RADS attributes and the patient's age. This dataset consists of 961 with six attributes.
27


Reviews:
[{"review_id": "4DtxMQqXXp", "reviewer": "Reviewer_bt3E", "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": "The paper brings an interesting topic to the audience of GP models, which is the problem of correctly calibrating the uncertainty in probabilistic non-linear regression and classification tasks. The use of more heavy-tailed distributions instead of just Gaussian densities makes sense to me, and authors motivated it better in this new version of the manuscript. The general presentation of the method and the technical details are clear, and easy to follow, apart from some details around the concept of identifiability and the theoretical conditions to be a stochastic process in Section 2.2, which are a bit short and unclear to me. (More details around these two points would be highly appreciated). On the modelling side, I saw the main point of the model, playing around with the parametric mixing distribution to control the elliptical model in both the likelihood and prior functions of the GP (or EP). The application of the proposed model to several datasets and comparison with SOTA methods is also a valuable point that I remark.", "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 proposes elliptical processes (EPs), a new family of stochastic processes inspired on GPs (Rasmussen & Williams, 2006) and Student-T processes (SPs) (Shah, ICML 2014). The underlying principle of these new EPs is that they are based on elliptical distributions, which are a more general family of distributions, typically more heavy-tailed, that includes Student-T, Cauchy as well as the Gaussian density. The key point of the elliptical distributions is their density generator, which in the end, depends of a parametric mixing variable \\xi. The authors consider this mechanism for both the function prior and the likelihood distribution, as they hope to better capture the uncertainty with the parametric heavy-tailed distribution. For the mixing distribution on the likelihood function, they consider a normalising flow based on splines. Experimental results show some evidence of the performance in regression, classification and heteroscedastic tasks, mainly in the well-known UCI datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Authors followed the requested changes {1,2} in the previous revision (see my previous comments as Reviewer CEkx), which have now improved quite a lot the quality of the manuscript. Additionally, removing the multi-task details made the manuscript much clearer and detailed, allowing the reader to focus on the main contributions and the technical utility of the proposed method.\n\nTo me, even if some extra experiments have been added, I still have a general concern around the limitation of the empirical results. To give an example, the performance of HetEP, EP-EP and EP-GP is kind of in-between the rest of SOTA metrics in Figure 7. Even if some analysis/conclusions are included in the end of Section 4.3, this is not entirely sufficient to me. I remark that I am not looking for metrics or datasets that indicate a better performance for the proposed model, but I would like to see that there is a clear utility on using the model indicated by some empirical results.\n\n> One advise would be to think about scenarios, datasets or problems (i.e. regression) where the use of Gaussians is predominantly worse than other heavy-tailed distributions. The authors indicated \"finance\" as one of the applications for more heavy-tailed stochastic processes, perhaps that could be a good one. But at least, I would like to see a real example of this effect and how the proposed EP correctly captures the whole uncertainty, improving standard GPs.\n\n> Other extra (minor) points could be to improve the explanations around the effect of using the elliptical distribution in the likelihood and the function prior, for example. Also the presentation of the sparse/inducing points could be fixed a bit, particularly the use of X_u and x_n, which seemed a bit odd to me.\n\n", "broader_impact_concerns": "I do not detect particular ethical implications of the work. ", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "8xJVvnPFBI", "reviewer": "Reviewer_j7Xo", "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\nSynthetic evidence : \n\nI liked the experiment of Figure 5 where the consequence of utilizing the more straightforward GP w/ Gaussian likelihood results in an overfit latent function with smaller estimated lengthscales. This is among the clearest examples provided for how using an elliptical distribution for noise can benefit modeling. The heteroscedastic experiment in Section 4.2 is also nicely illustrative.\n\n\nWeaknesses\n\nClarity : \n\nFor the use case where an elliptical distribution is used as observation likelihood, there is either reader confusion in interpreting the joint probability (5) and noise likelihood (6) or typos in these equations : The mixing distribution $p(\\omega; \\eta_\\omega)$ for the likelihood can be defined either per observation or hierarchically as (5) and (6) specify. In the former case, the mixing distribution can be estimated from observations IF the same elliptical distribution is used across observations and IF it happens to be true that the equivalence of two univariate elliptical distributions implies equivalence of their mixing distributions. In the latter (hierarchical) case, the mixing distribution cannot be estimated from observations for the same reason that it cannot be estimated when an elliptical distribution is used to model a random process, i.e., only a single draw of the mixing variable is observed. Figure 3a,c and the experiments of Figure 4 imply the paper is assuming the former model, but the a priori and a posteriori model specifications are for the latter model.\n\nFor the use case where an elliptical distribution is used within a random process prior, there are some inconsistencies in the exposition that create confusion. In Section 2.2, \"Prediction\" paragraph, in the statement \"The conditional distribution can be derived analytically (see Appendix B), but we will instead solve it by approximating the posterior $p(\\xi|y; \\eta_\\xi)\" what does \"conditional distribution\" refer to? The results of Appendix B refer to a pure elliptical distribution while the statement of Section 2.2 refers to conditioning on observations, i.e., drawn from a pure elliptical distribution plus i.i.d. noise.\n\nExperimental methodology : \n\nThe paper does not provide a procedure for how the number of bins in an a priori or a posteriori mixing distribution spine flow should be set and how an output transformation (softplus, sigmoid, arctan) after the flow should be specified. In other words, it should be clear to a reader whether or not the settings described in appendix G were tuned for the experiments (preferably, not tuned). Similarly, for implementing a heteroscedastic elliptical distribution noise model, how would one select an architecture for the input-varying network model $g_{\\gamma_\\omega}( )$ of Section 3.3 in a novel problem?\n\nRegression experiments :\n\nThe EP-EP results w.r.t. MSE plotted in Figure 7 are worse than the exact GP results in 5/7 datasets and not very different in the remaining 2. The claim that \"...log-likelihood is more critical than MSE because mean-squared error (MSE) does not consider the predictive variance\" is not unconditionally true. One would have to show that, for a given value of MSE, no value of variance could yield a given log likelihood. In other words, both metrics generally provide information. In the 5/7 cases where the exact GP has worse log likelihood, it could be argued that the exact GP is somehow poorly calibrated. The EP-GP method does not seem to perform any better than the EP-EP method across the metrics and datasets. The heteroscedastic EP method appears to be trading log likelihood for MSE even more than the vanilla EP-EP method. Finally, it isn't clear what additional value SVGP provides for interpreting the proposed method's performance.\n\nClassification experiments : \n\nThe evidence among the non-sparse variants in support of the proposed method seems marginal. One dataset (mammographic) shows little difference among the methods in both performance metrics. Another dataset (heart) shows some improvement in one metric (accuracy) but not much in the other (AUC). The final dataset (pima) does show what appears to be statistically significant improvement in one metric (AUC). In total, one could say that the use of EP distributions does not hurt modeling binary classification and can sometimes provide a little benefit. The included sparse variants are not necessary since the datasets are all small (<1000 sample size) and obfuscate interpreting the results. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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 proposes the use of the elliptical distribution for a priori modeling and a posteriori approximation of random processes and observations. Noting that a (consistent) elliptical distribution only arises as a scale mixture of Gaussians, the paper proposes modeling elliptical distributions generically using a normalizing flow for mixing distribution (for both a priori and a posteriori cases). A sparse variational solution for the case of elliptical process prior is provided. Synthetic experiments in regression exploring elliptical noise and heteroscedastic noise settings is provided as well as experiments on real regression and classification datasets.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Section 2.2, \"Prediction\" paragraph : Please correct the statement \"The conditional distribution can be...\"\n\nSection 3 : Please correct the likelihood definitions in (5), (6), (8), and (9) to include an $i$ subscript on the mixing variable or clarify how the mixing distribution can be estimated as in the experiments of Figure 4 from a single draw of the mixing distribution for all observations. In other words, I believe the correct expression for the likelihood term in (5) is $\\prod_{i=1}^N p(y_i | f_i, \\omega_i) p(\\omega_i; \\eta_\\omega)$.\n\nSection 3.2 : Please provide some motivatation for selecting a spline flow for the variational posterior $q(\\xi; \\phi_\\xi)$. Since the prior model (5) selects a single draw from the mixing distribution $p(\\xi; \\eta_\\xi)$ it isn't clear why a complicated posterior is necessary.\n\nSection 4.4 : If the paper would like to include a comparison among sparse variants, please select a dataset large enough to warrant them, e.g., the airlines dataset of the sparse GP classification paper by Hensman (2015) et al.\n\nSection 4.4 : Please comment on the estimated test-set log predictive densities of the two (full) EP methods compared to the GP method (or include results and discussion).\n\nAppendix D : Please explain the distinction b/w $\\eta$ in (48) and $\\eta_f$ in (49).\n\nAppendix G : Please elaborate on the choices of output transformations utilized in the likelihood mixing distribution (softplus) and posterior mixing distribution (sigmoid). The included hypothesis that the posterior is inherently \"more difficult to learn\" should be backed up with experimental evidence.\n\nMinor:\n\nSection 3 : \"...is a regular EP prior with the covariance kernel...\" -> a regular GP prior (?)\n\nSection 3.2 : Please modify the title \"Prior\" of this section. It is not obvious why this title applies.\n\nSection 3.2 : Extra vertical bar in (17) in expression $p(u, \\xi |; \\eta_u, \\eta_\\xi)$.\n\nSection 4.2, 1st para : The expression for the true function $f( x_n )$ is missing the additive noise.\n\nSection 4.4, 2nd para : \"Therefore, we can indicate the value of the elliptical...\" What does it mean to \"indicate the value\" of the process? Please clarify or re-word.\n\nSection 4.4, 3rd para : \"performed a ton-fold cross-validation\" -> ten-fold\n\nAppendix B : At the end of the line starting with \"Proof of proposition 1.\", the variable $M$ in $p(\\xi|y_1 M ; \\eta)$ is undefined.\n\nAppendix F : There is a missing $\\xi$ in the expression in the line above (66).", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "04aZdoirAs", "reviewer": "Reviewer_g9nQ", "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- The paper is well-written and a pleasure to read. I believe it does an especially great job in communicating the theoretical contributions required to formulate $\\mathcal{EP}$s, while also including multiple illustrative examples that convey the utility of the proposed model in practical set-ups. Crucially, the paper comes across as fairly \"complete\", and is especially enhanced by the ablation studies and discussion contained in Section 4. I also appreciated that the experiments highlight both strengths and weaknesses of the proposed model (such as the greater risk of overfitting when accounting for heteroscedastic noise).\n\n- The problem itself is very interesting and relevant, and as the authors indicate, highly-relevant to practitioners working with data having heavy-tailed distributions (without compromising on the effectiveness of GPs for normal-distributed data). There would definitely be large value in having the model proposed here being integrated in a more generic toolkit such as $sklearn$ in order to improve usability.\n\n**Weaknesses**\n\n- One aspect that sometimes comes across as unclear is the core contribution of itself, i.e. is it the fundamental formulation of elliptical processes, or the incorporation of variational inference to enable tractability? The initial description of $\\mathcal{EP}s$ is provided in the *Background* section, but there is no reference to the unpublished arXiv manuscript on a similar topic (https://arxiv.org/abs/2003.07201).\n\n- A small suggestion for the Experiments section is to perhaps include a critical difference diagram that can more succinctly summarise how GPs and different variations of $\\mathcal{EP}s$ compare against each other. Beyond commenting on computational considerations, a more formal description of computational complexity would also be welcome.\n\n- In the related work section, although there are several references to alternative models permitting more complex and flexible posteriors than standard GPs (such as deep GPs), these methods are not featured in the experimental evaluation itself. While I wouldn't expect the evaluation section to be extended further, a little more commentary on what performance can be expected form these competing models when applied to heavily-tailed distributions could be just as illuminating.\n\n- [nitpick] The final sentence of the paper reads as unnecessarily pessimistic, and could with either a citation or a supporting argument.", "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": "In this work, the authors present a novel interpretation of elliptical processes that leverages variational inference in order to approximate intractable likelihoods and reduce computational complexity. This class of models is expected to outperform Gaussian processes (GPs) on datasets having heavier-tailed distributions, while also subsuming the former along with Student-$t$ processes. Via an extensive combination of experiments on synthetic and real-world datasets, the authors demonstrate how different aspects of variational $\\mathcal{EP}$s (i.e. their likelihood, posterior approximation, etc) contribute towards the overall performance of the model in comparison to the more widely-used GPs.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The main changes I would like to see in a future revision of the paper are closely tied to the weakness listed in the previous section. For clarity, these are:\n\n- *[major]* Clarifying whether the specification of elliptical processes is a core contribution of the paper itself, or whether the focus is on the variational approximation. The current placement of the introduction to $\\mathcal{EP}$s in the Background section makes this unclear.\n\n- *[major]* Expanding the Related Work section to more clearly express connections to alternative methods such as deep GPs and deep kernel learning methods.\n\n- *[minor]* Slight tweaks to Experiments section to make key takeaways more clear.\n\n- *[minor]* Tiny spotted typos:\n\n    - pg.1 - \"Gaussian distributions *is* attractive...\"\n    - pg.12 - \"a *ten*-fold cross-validation where *we*...\"", "broader_impact_concerns": "I do not foresee any broader impact concerns resulting from this work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Overall, this paper is of interest to the community and the clarity and focus of the paper are much improved from the previous submission. Indeed, those reviewers who reviewed both versions were very complementary about the improvements in this version.

While there are some concerns about the experimental results (which show marginal improvements over alternative methods), the reviewers all agree that EPs, and the proposed methodology, are of interest to the community. Therefore, I do not request any further experiments.

I *would* like to see clearer delineation of the contributions. Two out of three reviewers were of the impression that the Elliptical process was a new concept; however as the authors address in response to reviewer g9nQ, the idea exists already in the statistical literature (albeit, not as a practical modeling method). Making EPs a practical modeling tool is certainly an important contribution, and under "Our Contributions", the authors do indeed say that their contribution is a new construction of EPs, not the introduction of EPs. However, I would like to see this distinction made clearer in the Introduction and in Section 2.2 (where it is unclear whether Definition 1 is a new contributions since it has no attribution). 

In addition, the reviewers noted that, while the paper is much approved, it would benefit from a final proof read -- for example, in (5) the subscript should be $n$ not $i$
Decision: Accept (poster)