Computationally-efficient initialisation of GPs: The generalised variogram method
Authors: Felipe Tobar
Elsa Cazelles
Taco de Wolff
Venue: Accepted by TMLR ()

Abstract: We present a computationally-efficient strategy to initialise the hyperparameters of a Gaussian process (GP) avoiding the computation of the likelihood function. Our strategy can be used as a pretraining stage to find initial conditions for maximum-likelihood (ML) training, or as a standalone method to compute hyperparameters values to be plugged in directly into the GP model. Motivated by the fact that training a GP via  ML is equivalent (on average) to minimising the KL-divergence between the true and learnt model, we set to explore different metrics/divergences among GPs that are computationally inexpensive and provide hyperparameter values that are close to those found via ML. In practice, we identify the GP hyperparameters by projecting the empirical covariance or (Fourier) power spectrum onto a parametric family, thus proposing and studying various measures of discrepancy operating on the temporal and frequency domains. Our contribution extends the variogram method developed by the geostatistics literature and, accordingly, it is referred to as the generalised variogram method (GVM). In addition to the theoretical presentation of GVM, we provide experimental validation in terms of accuracy, consistency with ML and computational complexity for different kernels using synthetic and real-world data.

Published in Transactions on Machine Learning Research (04/2023)
Computationally--ffcceet initialisation The generalised variogram method
inittiacooooniiitiaaiitttiooo of GPs:
Felipe Tobar
ftobar@uchill... cl
Initiative for Data & AT Universidad de Chile
Elsa Cazelles
elsa.cazelles@hrrt.fh
CNRS, IRIT
Université de Toulouse
Taco de Wolff
tacodewolfi@gmace com
Inria Chile
Reviewed on OpenReview: https: // openreview. net/ forum? id= ad=slaae slaAHHP5h
Abstract
We present a computationalll--ffiiein strategy to initialise the hyperparameters of a Gaus- sian process (GP) avoiding the computation the likelihood function. Our strategy can be used as a pretraining stage to find initial conditions for maximun-likelihod (ML) training, or as standalone method compute hyperparameters values to be plugged in directly into the GP model. Motivated by the fact that training GP via ML is equivalent (on aver- age) to minimising the KL-divergence between the true and learnt model, set to explore different metrics/divergenn
among GPs that are computationall inexpensive and provide hyperparameter values that are close to those found via ML. In practice, we identify the GP hyperparamoters by projecting the empirical covariance or (Fourier) power spectrum onto a parametric family, thus proposing and studying various measures of discrepancy operating on the temporal and frequency domains. Our contribution extends the variogram mothod developed by the geostatistics literature and, accordingly, it is referred as the generalised variogram method (GVM). In addition to the theoretical presentation of GVM, pro- vide experimental validation in terms of accuracy, consistency with ML and computational complexity for different kernels using synthetic and real-world data.
1 Introduction
Gaussian processes (GPs) are Bayesian nonparametric models for time series praised by their interpretability and generality. Their implementation, however, is governed by two main challenges. First, the choice of the covariance kernel, which usually derived from first principles or expert knowledge and thus may result in complex structures that hinder hyperparameter learning. Second, the cubic computational cost of standard, maximun--ikclhhood-besse thousands observations. The GP community actively largets these issues though the development of robust, computationally--ffiin efficicien training strategies. Though these advances have facilitated the widespread use of GP models in roalistic settings, their success heavily depends on the initialisation of the hyperparameters.
training which renders the exact GP unfeasible for more than a few
In practice, initialisation either follows from expert knowledge or time-consuming stochastic search. This is in sharp contrast with the main selling point of GPs, that is, being agnostic to the problem and able freely learn from data. To provide rescarchers and practitioners with an automated, general-application and cost-efficient initialisation methodology, we propose learn the hyperparameters by approximating the empirical covariance by parametric function using divergences bet covariances that inexpensive

---

Published in Transactions on Machine Learning Research (04/2023)
to compute. This approach is inspired by the common practice, such when computes some statistics (e-g-, mean, variance, discrete Fourier transform) to determine initial values for the kernel hyperparameters In the geostatistics literature, method that follows this concept is the Variogram (Cressie) 1993 Chiles & Delfiner 1999). which is restricted to particular cases of kernels and divergences. Therefore, we refer to the proposed methodology as Generalised Variogram Method (GMV) the sense that it extends the application of the classic methodology to a broader scenario that includes general stationary kernels and metrics, in particular, Fourier-bassee metrics. Though our proposal is conceptually applicable to an arbitrary input/outpin dimension, focus the scalar-input scalar-output case and leave the extension to general dimensions as future work H
The contributions of our work are of theoretical and applied nature, and are summarised as follows:
novel, conputationalll--ffiiiini training objective for GPs as an alternative to maximum likeli- hood, based on projection of sample covariance (or power spectrum) onto parametric space;
a particular instance of the above objective, based on the Wasserstein distance applied to the power spectrum, that is convex and also admits a closed-form solution which can thus be found in a single step;
a study the computational cost of the proposed method and its relationship to maximum likeli- hood, both in uaucept and empirical terms;
experimental validation of the proposed method assessing its stability with respect to initial con- ditions, confirming its linear computational cost, ability train GPs with a large number of hyperparameters and its advantages initialiser for maximunn-likelhood training.
2
Preliminaries
2.1
Motivation
Let us consider ~ GP(O, K) and its realisation y [vi...., Vn] E IR" at times t == [ t,] R". The kernel K is usually learnt by choosing family beee and optimising the liglliillhood
1(0) - Tr (K, "y log Kg log 2nt
(1)
with respoct to 0 E O where we used the cyclic property of the trace, and defined K def Ko(t) according to [Kols Ko(t: i, nI.. Defining K def K(t) ByyT, we note that
EI(A) = - Tr (K₂ K) log K log 2n.
(2)
2
2
Notice that, up to terms independent of 6, equation is equivalent the negative Kullback-Leibler divergence (NKL.) between the (zero-mean) multivariate normal distributions NN K) and (O, Ko) given by
K K
DNKL(N (O)
2 Tr (K₂ K) + log
n
(3)
which, with a slight abuse of notation, can be expressod as a function of only the covariances as DNKL(K, Ko) def DN K)) Ka)).
This reveals that learning a GP by maximising I(0) in equation I can be understood (in expectation) as minimising the KT. between the t-marginalisations of the true process GP(O, K) and a candidate process GP(O, Ko). This motivates the following remark.
Remark 1. Since macimum--lkkelttun learning of GPs has a cubic computational cost but it (on average) equivalent to minimising a KL divergence, what other divergences or distances, computationally cheaper than the likelihood, can be considered for learning GPs?
For illustraionn and completeness, we incorporate toy experiment featuring 5-dimensional input data Sec. 5.7
2

---

Published in Transactions on Machine Learning Research (04/2023)
2.2 Divergences over covariance functions
We consider zero-mean stationary GPs 2 The stationary requirement allows us to i) aggregate observations in time when computing the covariance, and ii) compare covariances in terms of their (Fourier) spectral content. Both perspectives will be present throughout the text, thus, we consider two types of divergences over covariances: i) temporal ones, which operate directly to the covariances, and ii) spectral ones, which operate over the the Fourier transform of the covariances, i.c., the GP's power spectral density (PSD). As our work relies extensively on concepts of spectral analysis and signal processing, a brief introduction to the subject is included in Appendix
Though we can use most metrics (c.g., L1 L2) on both domains, the advantage of the spoctral perspective is that it admits -ensity- divergences as it customary in signal processing (Basseville 1989). Bregman divergences, which include the Euclidean, KL and Itakura-Saito (IS) (Itakura, 1968), are vertical measures, i., they integrate the point-wise discrepancy between densities across their support. We also consider horizontal spectral measures, based on the minimal-cost to transport the mass from one Iistribution-across the support recently been considered for comparing PSDs using, e.g., the 1- and 2- 2-Wasserstein distances, denoted W1 and W2 (Cazelles al. 2021 Henderson & Solomon 2019 Soe Appendix B for definitions of vertical and horizontal divergences.
another. This concept, known as optimal transport (OT) (Villani 2009) has only
2.3 Related work
A classic notion of kernel dissimilarity in the machine learning community is that of kernel alignment, concept introduced by Cristianini et al. (2001) which measures the agrecment between two kernels or a kernel and a task; this method has been mostly applied for kernel selection in SVM-based regression and classification. This notion of dissimilarity is based on the Frobenius inner product between the Gram matrices of each kernel given a dataset- see (Cristianini et al. 2001] Def. 1). Though, in spirit, this concept is related to ours in that the kernel is learnt by minimising a discrepancy measure, we take a signal-proccssing inspired perspective and exploit the fact that, for stationary GPs, kernels are single-input covariance functions and thus admit computationalll-effiiees discrepancy measures. In addition to the reference above, the interested reader is referred to (Cortes et al. 2012) for the centred kernel alignment method.
Within the GP community, two methodologies for accelerated training can be identified. The first one focuses directly on the optimisation procedure by, e.g., avoiding inverses (van der Wilk et al.] 2020), or derivatives (Rios & Tobar 2018), or even by parallelisation; combining these techniques has allowed to process even a million datapoints (Wang et al. 2019). A second perspective that of sparse GP approxima- tions using pseudo-inputs Quinonero-Candela & Rasmussen 2005), with particular emphasis on variational
methods (Titsias 2009 The rates of convergence of sparse GPs has been studiod by Burt et al. (0199)
and the hyperparameters in this approach have also been dealt with in Bayesian manner by Lalchand
et al. 20022) Sparse GPs have allowed for fitting GPs to large datasets (Hensman et al.) 2013 training
non- non-paramctric kernels Tobar et al.] 2015), and implementing deep GPs (Damianou & Lawrence 2013
Salimbeni & Deisenroth, 2017). Perhaps the work that is closest to ours in the GP literature is that of Liu
et al.] (2020), which trains a neural network to learn the mapping from datasets to GP hyperparameterss thus avoiding the computation of the likelihood during training. However, the authors state that " "training on very large datasets consumes too much GPU memory or becomes computationally expensive due to the kernel matrin size and the cost of önwerling the malria". This is because they still need to compute the kernel matrix during training of the net, while we bypass that calculation altogether.
The Wasserstein distance has been used to compare laws of GPs MMsaarotto et al. 2019 Mallasto & Feragen 2017). , and applied to kernel design, in particular to define GPs (Bachoc et al. 2018) and deep GPs (Popescu ct al. p over the space of probability distributions. Cazelles et al. (2021 proposed a distance between time series based on the Wasserstein, termed the Rosserstein-Pooier distance, and showed its application
We do so for convenience in computation since non-zcro-mean CPs can be understood mixture of a zero-mean CP and a parametrie regression model, where the GP models the residuals of the parametric part. See Section 2.7 in Rasmussen & Williams (2005) for discussion.
3

---

Published in Transactions on Machine Learning Research (04/2023)
to GPs. However, despite the connection between GPs and the Wasserstein distance established by these works, they have not been implemented to train (or initialise) GPs.
In geostatistics, the variogram function (Cressie, 1993 Chiles & Delfiner 1999 is defined as the variance of the difference of a process y at two locations f and t2, that is, 7(1 - 12) V[y(tz) - y(t2)] The variogram is computed by choosing a parametric form for (() and then fit to a cloud of points (sample variogram) using least squares. Common variogram functions in the literature include exponential and Gaussian ones, thus drawing natural connection with GP models. Furthermore, when the process is stationary and isotropic (or one-dimensional)
in the GP models considered here, the variogram and the covariance K(t) follow the relationship ((t) ((0) (t),, therefore, given a kernel function the corresponding variogram function can be clearly identified (and vice versa). The way in which the variogram is fitted in the geostatistics literature is what inspires the methodology proposed here: we match parametric forms of the covariance (or the PSD) to their corresponding sample approximations order to find appropriate values for the kernel hyperparameters. Also, as we explore different distances for the covariance and PSD beyond the Euclidean one (least squares) we refer to our proposal the Generalised Variogram Melhod (GVM).
GVM complements the literature in way that is orthogonal to the above developments. We find the hyperparameters of GP in likelihood-free manner by minimising a loss function operating directly on the sample covariance or its Fourier transform. As we will see, GVM is robust to empirical approximations of the covariance or PSDs, admits arbitrary distances and has a remarkably low computational complexity.
3 A likelihood-free iwariance-matchiin strategy for training GPs
Overview of the section. As introduced in Section 2.1 objective is to approximate the ground-truth kernel K with parametric kernel Kø-; to this end we will rely the kernel denoted Kn, We will proceed by matching K, with parametric form Ko; using metrics in the temporal domain, i.c., solving equation 7 or in the spectral domain matching the Fourier transform of Rn, denoted S2 F FK,, with a parametric form S i.e., solving equation [8] In the following, we present the Fourier pains Kg and Sa, the estimators Kn and Sn, and the metrics considered for the matching, We then present an explicit case for location-scall families, and we finally give theoretical arguments for the proposed learning method.
empirical data-driven estimate of
3.1 Fourier pairs Ko and S, and their respective estimators
Let us consider the zero-mean stationary process y GP(O, Ko) with covariance K2 and hyperparameter f E O. If the covariance Kg is integrable, Bochner's 's Theorem (Bochner 1959) states that Kg and the process' power spectral density (PSD), denoted Su, are Fourier pairs, that is,
def / Kolt)e (olt)e-
So(E) - FF
12nEt di,
(4)
where j is the imaginary unit and F {.} denotes the Fourier transform.
Since zero mean stationary GPs are uniquely determined by their PSDs, any distance defined on the space of PSDs can be "lifted" to the space covariances and then to that of GPs. Therefore, we aim at learning the hyperparameter 0 E θ by building statistics in either the temporal space (sample covariances) or in spectral space (sample PSDs).
First, we consider the following statistic Kn in order approximate the ground-truth kernel K.
Definition 1. Let E R be zero mean stochastic process over R with observations y [V....., Vnl E R. al limes f. = 14,,..., 1a] E R The empirical covariance of y is given by
kn(()
Σ
Viyl 1 Card( = =1 1,
(5)
=1=1
4

---

Published in Transactions on Machine Learning Research (04/2023)
The estimator (t) the kernel Kn in equation that is
the PSD S then simply corresponds to applying the Fourier transform to the empirical
Sn(E) / (t)e j2 dt.
(6)
Another usual choice for estimating the PSD is the Periodogram Sper (Schuster 1900). Though Sper is asymptotically unbiased (VV, Per S()),, it is inconsistent, i.c., its variance does not vanish when oo (Stoica & Moses 2005 [Sec. 2.4.2]. Luckily, the variance of Spee(E) can be reduced via windowing
and the Welch/Bartlett techniques which produce (asymplotically) consistent and unbiased estimates of S(E), vg Stoica & Moses
2005
3.2 Fourier-based covariance divergences
The proposed method builds on two types of (covariance) similarity. First, those operating directly on Rn and Kg, known as temporal divergences, which include the L1 and L2 distances. Second, those operating over the Fourier transforms of K,, and Kg. that S, and Sg, known as spectral divergences. In the temporal case, we aim to find the hyperparameters of y by projecting Kn (t) in equation 5 onto the parametrised family K (Ko, O 0}. That is, by finding 0* such that Ko(-) is as close as possible to kn(t), i..,
0 arg min D(Ra, Ko);
l
where the function D(:-)) is the chosen criterion for similarity. Akin to the strategy of learning the hy- perparameters of the GP by matching the covariance, can learn the hyperparameters by projecting estimator of the PSD, namely S, in equation onto a parametric family S = {So, 0 E O). that is,
0 = arg min
(8)
OE
where Dp(:,)) is a divergence operating on the space of PSDs.
Remark 2. Since the mop Kg + So is one-to-one, equation 7 and equation 8 are equivalent when Dr(s,.) D(F F()) and S [KC..
We will consider parametric families with explicit inverse Fourier transform, since this way parametrises both the kernel and the PSD and can be learnt in either domain. These families include the Dirac delta, Cosine, Square Exponential (SF),, Student's f, Sinc, Rectangle, and their mixtures; see Figure 10 in Appendix A for an illustration of some of these PSDs and their associated kernels.
Remark 3 (Recovering kernel paramcters from PSD parameters). For parametric Fourier pair (K, S), the Fourier transform induces map between the kernel parameler space and the PSD parameler space. For kernel/PSD families considered in this work this map will be bijective, which allows us to compute the estimated kernel parameters from the estimated PSD parameters (and vice versa); see Table 1 for two examples of parametric kernels and PSDs. Furthermore, based on this bijection we refer as O to both kernel and PSD parameters.
3.3 A particular case with explicit solution: location-scale family of PSDs and 2-Wasserstein distance
Of particular relevance to our work is the 2- 2-Wasserstein distance (W2) and location-see scale families of PSDs, for which the solution of oquation X is completely determined through first order conditions.
Definition 2 (Location scale). A fomily of one-dimensional integrable PSDs is said to be of location-scaa type if it is given by
1 S (p
(5-1)) , ..E ..ER,OER (
(9)
where μ E R. the location parameter, R4 the scale parameter and Su,1 is the prototype of the family.
For arbitrary prototypes So,1, location-scale families of PSDs are commonly found in the GP literature. For instance, the SE, Dirac delta, Student's i, Rectangular and Sine PSDs, correspond to the Exp-cos, Cosine,
5

---

Published in Transactions on Machine Learning Research (04/2023)
Table 1: Location-scale PSDs (left) and their covariance kernel (right) used in this work. We have denoted by t and & the time and frequency variables respectively.
PSD
(4))(()
Prototype (E)
Kernel
Kt,o(t) VTo exp(-o² x cos(2 cos=mt
Ko, v <rexp(==22 sinc(t)
Square-exp
exp
exp(-E")
Exp-cos
Rectangular
rect
rect(()
Sinc
osinc(ar) c0s(2
Laplace, Sine, and Rectangular kernels respectively. Location-scale families do not, however, include kernel mixtures, which are also relevant in our setting and will be dealt with separately. Though the protolype So,1 might also be parametrised (e.g., with a shape parameter), we consider those parameters to be fixed and only focus on == (μ, for the rest of this section. Table shows the two families of location-scale PSDs (and their respective kernels) that will be used throughout our work.
Remark 4. Let consider location-scale family of distributions with prototype So, and arbitrary member Sy Their quantile (t.e., inverse cumulattve) functions, denoted Q0 and Que respectivcly, obey
Quo(p) = = (00,
(10)
The linear expression in equation 10 is pertinent in our setting and motivates the choice of the 2-Wasserstein distance W2 to compare members of the location-scale family of PSDs. This is because for arbitrary one- dimensional distributions S₁ and S2, W3(S),, S2) can be expressed according to:
w3 (2) -
(11)
where Q1 and Q2 denote the quantiles of S1 and S2 respectively. We are now in position to state the first main contribution of our work.
Theorem 1. S is location-scale family with prototype Su,,, and S is arbitrary PSD, the minimiser (m",a") of W2(S, is unique and given by
n² -==
=== @(D)do
and a² So Qb L Q(1)
(12)
where Q is the quantile function of S.
Proof Sketch. The proof follows from the fact that (2(S, Sp.o) is convex both j and a, which is shown by noting that its Hessian is positive via Jensen's 's inequality. Then, the first order conditions give the solutions
in equation 12 The complete proof can be found in Appendix C
Remark 5. Although integrals of quantile functions are not usually available in closed-form, computing is straightforward. First, (p)(p)) is known for large class of normalised protolypes, including square exponential and reclangular functions. Second, p¹ Ex-s(a] and fo (r) Ex~s[=²), where x is random variable with probability density S. Third, both integrals are one-dimensional and supported [0, 1], thus numerical integration is inexpensive and precise, specially when S has compact support.
equation
As pointed out by Cazelles et al. (2021) W Sy) is in general non-convex in 0, however, for the particular case of the location-scale family of PSDs and parameters (), convexity holds. Since this family includes [ guarantees closed form solution for the optimisation problem in equation and thus can be instrumental for learning GPs without computing (the expensive) likelihood function.
usual kernel choices the GP literature, Theorem
3.4 Learning from data
Learning the hyperparameters through the optimisation objective in equation or equation g possible provided that the statistics K,, and S, converge to K and S respectively. We next provide theoretical

---

Published in Transactions on Machine Learning Research (04/2023)
results on the convergence of the optimal minimiser 8% The first result, Proposition 1 focuses the particular case of the 2- Wasseesstei distance and location-scale families, presented in Section 3.3
Proposition 1. Let So be location-scale family, S the ground truth PSD and DF - W Then (S, Sw) + 0 implies that the empirical minimiser 0% in equation converges to the true minimiser 0r = arg mingge (S, So), meaning that E - + -
Proof Sketch. First, in the location-scale family we have 0 = (1, σ). Then, the solutions in Theorem alllow us to compute upper bounds for 1,2 -1= and - via Jensen's 's and Hölder's inequalities, both of which converge to zero as (S?(S, S..) + The complete proof can be found Appondix C
The second result, Proposition 2 deals with the more general setting of arbitrary parametric families, the distances L1, L2, W1, W2 and cither temporal and frequency based estimators. To cover both cases, we will denote In to refer to either Kn Sn. Let also consider the parametric function so E {Solo O) and the ground-truth function f which denotes either the ground-truth covariance ground-truth PSD.
Proposition 2. For general parametric famálies {Sa θ), the empirical solution 0% conwerges the true solution H under the following sufficient and stronger conditions: (i) D = W, or L,,7 = 1,2, (ii) f) f)
to
0.8. 0; (iit) On + D(fe,, fo) O; and (tv) the parameter space θ is compact.
Proof. This result follows as a particular case from Theorem 2.1 in Bernton et al. 2019 where the authors study general r- Wasserstein distance estimators for paramotric familics of distributions for empirical mca- sures, under the notion of T-convergence or, equivalently, epi-convergence. The proof for the Le, 1, case is similar to that of W,.
4 Practical considerations
This section is dedicated to the implementation of GVM. We first discuss the calculation of the estimators Rn and Sn and their complexity, we then present the special case of location-scale family and spectral 2- Wasserstein loss. The last subsection dedicated to the noise variance parameter and the relationship between GVM and ML solutions.
4.1 Computational complexity of the estimators Rn and S..
First, for temporal divergences (operating on the covariance) we can modify the statistic in equation 5 using binming, that is, by averaging values for which the lags are similar; this process automatic the case of evenly-sampled data and widely used in discrete-time signal processing. This allows to reducing the amount of summands in K, from n(a+1) to an order n or even lower if the range of the data grows beyond the length of the correlations of interests.
2
Second, the Periodogram Sper can be computed at a cost of O(nk) using bins, where n the number of observations and k the number of frequency bins; in the evenly-sampled case, one could set k n and apply the fast Fourier transform at cost On logn). However, for applications where the number of datapoints greatly exceeds the required frequency resolution, & can be considered to be constant, which results in cost linear in the obsorvations O(n).
4.2 The location-scale family and 2- 2-Wasserstein distance
Algorithm 1 presents the implementation of GVM in the case of location-scale family of PSDs (Sg)ece and 2- Wasserstein loss presented in Section 3.3
7

---

Published in Transactions on Machine Learning Research (04/2023)
Algorithm 1 GVM Spectral loss W2 & ([So)oee is location-scale)
Require: 4,, .,....7 .
Require: location-scale PSD family (Se) (egg, Gaussians)
Define a grid over frequency space g (1o.
Compute S,, = (g), from eq. (6), over the chosen frequency grid Compute On the quantile function of S,
Compute 8r given by eq. 12 with Q = Qn
Linear computational complexity and quantiles. The cost. of the 2-Wasserstein loss is given by calculating i) the Periodogram Sper of the estimator S, in equation 6 ii) its corresponding quantile function, and iii) the integrals in equation 12 Though the quantile function is available in closed form for some famílics of PSDs (see Remark 5 they can also be calculated by quantisation in the general case: for a sample distribution (or histogram) one can compute the cumulative distribution and then invert. it to oblain the quantile. Additionally, the integrals of the quantile functions also have a linear cost but only in the number of frequency bins ((k) since they are frequency histograms, therefore, computing the solution has linear cost in the data.
4.3 Solving the general case: beyond the location-scale family
Recall that the proposed method has closed-form solution when considers PSDs in location-scale family and the W2 distance over the spectral domain. However, in the general case in equations 7 and S the minimisation is not convex in 0 and thus iterative/numeriaal optimisation is needed. In particular, for horizontal Fourier distances (c.g., Wasserstein distances) the derivative the loss depends on the derivative of the quantile function Qs with respect to 8, which might not even be known in closed form in general, specially for mixtures. However, most cases the lags of the empirical covariance or the frequencies of the Periodogram belong to a compact space and thus numerical computations are precise and inexpensive. Therefore, approximate derivatives can be considered (e.g., in conjunction with BFGS) though in our ex- periments the derivative-free Powell method also provide satisfactory results. Algorithms 2 and 3 present the application of GVM using the temporal and spectral metrics respectively. Note that we will usually consider L1, L2 as temporal metrics and also Wi, W2, Itakura- Saito and KL as spectral divergences Dr. Algorithms and are presented side to side for the convenience of the reader.
Algorithm GVM Temporal loss Require: L;, V/i, ,,
Algorithm GVM Spectral loss
Require: 4, Di, 1,...,
Require: parametric family {Ko]eee
Require: parametric family {Sy)eee
Require: temporal metric D)
Require: spectral metric Dp(:-))
Define temporal grid t (t1, (2.., tn
Define frequency grid g = (Eo. E............ {
Compute K₂ Ko(t)
Compute Sn So(g)
Compute Kn (t) in eq.
5
Compute S, S,)g), from eq.
6
Construct loss 0 - D(K,, Ko)
Construct loss 0 - (r(Ss Sn)
Find 0*, in eq. using BFGS or Powell
Find 0% in eq. using BFGS or Powell
Linear computational complexity. For the general case of spectral losses using numerical optimisation methods, we need to calculate Š, or its quantile (which are (n)) only once, to then compute the chosen distance Dr, which is Okk) for discrete measures defined on a k-point frequency grid as many times as the optimiser requires it. Therefore, the cost is Okk) but with a constant that depends on the complexity of the parametric family {S,, Q} and the optimiser of choice. For spatial losses, the complexity follows the same reasoning for the chosen spatial distance D and parametric family (Ke, O). and thus also linear in the datapoints.
8

---

Published in Transactions on Machine Learning Research (04/2023)
4.4 Noise variance and relationship to maximum likelihood
Following the assumplions of the Fourier transform, the spectral divergences considered apply only to Tebesgue-interrable PSDs, which rules out the relevant case of white-noise-corrpppee observations. This is because white noise, defined by a Dirac delta covariance, implies a PSD given by a constant, positive, infmite-support, spectral floor that is non-integrable. These cases can be addressed with temporal diver- gences, which are well suited (theoretically and in practice) to handle noise.
The proposed hyperparameter--earrd method intended both as a standalone likelihood-free GP learning technique and also routine. In this sense, we identify relationship between the ML estimator 0 and the proposed estimator 0" arg mingee (2(S, So). From equation 4 and Plancherel's theorem, we have || So ra = || Koll Ag Then, by definition of the estimators and Lemma 2 in Hoffman & Ma 2020 obtain the following inoquality
initialisation approach to feed initial conditions to maximum likelihood (ML)
1 | "SK-"1lall |aal|K ||
DKL. [K|| ) DKL Kg.)
Ka-
(13)
where || denotes the matrix Frobenius norm, DKI. (AIlBB denotes the KL divergence between zero-mean multivariate normal distributions with covariances A and B, and recall that K is the kernel of the ground truth GP. Denoting the ball centred at with radius M by BO0, M), we present the following remark.
Remark 6. The inequality equation 13 states that if the proposed estimator 4* is such that ||S ll,, e B(0, M), then DKL |K|| E B 2 K ||22||/7711a Therefore, under the reasonable assumplion that the function 0 + Kg only produces well-conditioned matrices, the factor ||K K ll2 is bounded and thus both balls have radius of the same order.
5 Experiments
This section illustrates different aspects of the proposed GVM through the following experiments (E):
El: shows that GVM, unlike ML, is robust different initialisation values and subsets of observations,
E2: studies the sensibility of GVM to the calculation of Sn, using Periodogram, Welch and Bartlett,
E3: validates the linear computational complexity of GVM in comparison to the full and sparse GPs,
E4: compares spectral against temporal implementations of GVM on an audio time series,
E5: exhibits the results of learning a 20-component spectral mixture GP using GVM with different spectral metrics,
E6: assossos the ability of GVM to produce initial conditions which are then passed to an ML routine that learns spectral mixture GP with 4, 8, 12 and 16 components,
E7: presents toy example where GVM is used to fit GP with an isotropic SE kernel to multi-input (5-dimensional) data.
The code for GVM is available at /attps://github. httpps://gggaaa....o///lSSSSSSSL Cctp://HMSSSUtiiiic/eeaaaaaaaaaa minimal working example of the code presented in Appendix D. The benchmarks were implemented on MOGPTK (de Wolff et al. 2021
5.1 E1: Stability with respect to initial conditions (spectral mixture kernel, L2. temporal)
This experiment assessed the stability of GVM with respect to random initial conditions and dif- ferent realisations. We considered a GP with a 2-component spectral mixture kernel K(t) Σ of exp( exp(-7)7)) 27)7) c0T(2T)))) + o² &
8 with hyperparameters σi = 2, 71 = 10-4, μ = 10 2 We produced 4000-point rcalisations from the GP and 50
02 = 2,72 10-4 H2 10 10-2 T noise random initial conditions {Ør) j according to [#"), ~ Uniform[]e, }0., where A, is the i-th true hyperpa- rameter.
9

---

Published in Transactions on Machine Learning Research (04/2023)
fr=1 using ML, and ii) compute GVM from {0r),² and then perform ML. starting from the value found by GVM. Each procedure was implemented using a single realisation (to test stability wrt 0r) and different realisations (to test stability wrt the data). Our estimates Og were assessed in terms of the NLL and the relative mean absolute error (RMAE) of the parameters Ei 00 00-0,11 1//cll
We considered two settings: i) train from p0r
Fig. shows the NLL (left) and RMAE (right) versus computation time, for the cases of sixed (Lop) and different (bottom) observations; all times start from t = 1 and are presented in logarithmic scale. First, in all cases the GVM initialisation (in red) took about half a second and resulted in an NLL/RMAP virtually identical to those achieved by ML initialised by GVM, this means that GVM provides reliable parameters and not just. initial conditions for MT. Second, for the fixed observations (top), the GVM was stable wrt Or unlike ML which in some cases diverged. Third, for the multiple observations (bottom) GVM-initialised ML ML diverged.
diverged in two (out of 50) runs, which is far fewer than the times that
NLL: fixed data & random initial condition (50 runs)
relative abso lute error: fixed data & random initial condition (50 runs)
6x
ML from random
10
ML from random
ML. from
Inttnnisation
M. from Hntt
402
20'
3x1
10*
2'100
10a
10a
10P
16t
time [seconds]
time [seconds]
NLL: random data & random initial condition (50 runs)
relative abso lute error: random data & random initial condition (50 runs)
6'000
from rendom M. from in Intialtsation
107
HL random
init
Inkialsation
4x100
101
3'
10a
10
101 time [seconds]
10F
16
time [seconds]
Figure 1: GVM as initialisation for ML starting form 50 random initial conditions: proposed GVM (red), standard ML. (green) and ML. starting from GVM (blue). The top plots consider single dataset, while the bottom plots consider different realisations for each initial condition. The L-BFGS-B optimizer was used with the default gradient tolerance in order for the results to be comparable.
5.2
E2: Sensibility of GVM wrt the Periodogram (exact case: Wz. location-scale)
We then considered kernels with location-scale PSDs and the W2 metric over the PSDs. This case has a unique solution but requires us to compute S, in equation 8 this experiment evaluates different ways of doing so. We produced 4000 observations evenly-sampled in [0, 1000] from GPs with square exponential and rectangular PSDs (which correspond to the Exp-cos and Sinc kernels respectively see Table D each with location μ ~ U[0.025, 0.075) and scale ( U/0.01, 0.02]. We computed S via the Periodogram, Welch and Bartlett methods with different windows. Table 2 shows the percentage relative error (PRE) y averaged over 50 runs. The estimates of both parameters are fairly consistent: their magnilude. does not change radically for different windows and Periodogram methods. In particular, the estimates of the location parameter are accurate with an average error in the order of 2% for both kernels. The scale parameter was, however, more difficult to estimate, this can be attributed to the spectral energy spread in the Periodogram, which, when using the Wasserstein-2 metric, results in the scale parameter to be overestimated. For both kernels and и 0.05, 0.01, the case of the (windowless) Periodogram is shown in Fig, 2 and the remaining cases are all shown in Appendix E In the light of these results, we considered the Periodogram (no window) for the remaining experiments.
3PRE 110 (8-8) where is the true parameter and 9 its estimate.
10

---

Published in Transactions on Machine Learning Research (04/2023)
Performance of GVM learning GPs with the Exp-cos and Sine kernels under different sampling settings, Periodogram methods and windows. Each entry of the table shows the average percentage relative error for location (left) and scale (right), with their standard deviations, separated by the symbol '/", True parameters where Uf0.025, 0.075] and U(0.01, 0.02]: averages and standard deviations computed over
Table 2:
50 runs.
Kernel
Window
Periodogram
Bartlett
Welch
none
2.00=1..1/33.41113222 2.17
2.6511.4/24487144333
Exp-cos
hamn
2.98 2,0852.114/34.0-114.333
3.02:2.:3//2223311..00
2.52=1..75/31.8831133322
hamming
2.924-2.001/34.8--111..6
2.00-2..9/22:35+13322
2.40-1.78/82.22+411322
2.8.91.66//.9986..3
2.63:223/33.232226
2.31=1.77/38.57115..9
Sinc
hann
2.68::2.02/10.001-...22
2.58-1.77/58886-117..8
2.111.78/19.52211...22
hamming
2.0022.1199.22--.333
2.5621.755550.21155.77
2,13:1.76/165511100777
Learnt parameters 0.0501, scale: 0.0089
Learnt parameters 0.0506,
scale: 0.0108
150
Periodogram (w: None)
Periodogram (w: None)
200
kornel
100
Ground truth Exp-cos kemel
Ground truth Sinc keme!
100
0.00
0.01
0.02
0.03 frequencies
0.04
0.05
0.06
0.07
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
frequencies
Figure 2: GVM estimates for Exp-cos (left) and Sine (right) kernels shown in red against (blue) and true kernels (black). This case: Periodogram, no window.
5.3 E3: Linear complexity (exact case)
We then evaluated the computation time for the exact case of GVM (W2 distance and location-scale family) for increasing amount of observations. We considered unevenly-sampled observations from an single component SM kernel (μ 0.05, 0.01) in the range [0, 1000). We compared GVM against i) the ML estimate starting from the GVM value (full GP, 100 iterations), and ii) the sparse GP using 200 pscudo inputs (Snelson & Ghahramani 2006). Fig. 3 shows the computing times versus the number of observations thus validating the claimed linear cost of GVM and its superiority wrt to the rest of the methods. The (solid line) interpolation the plot is of linear order for GVM, linear for sparse GP since the number of inducing points fixed, and cubic for the full GP.
Empirical training cost
102
GVM
Seconds (log scale)
(ML) Sparse
101
10
10-1
200
400
600
800
1000 Observations
1200
1400
1600
1800
2000
Figure 3: Training times number of datapoints for the proposed GVM, full GP and sparse GP.
5.4 E4: Performance and cost of spectral and temporal metrics over the same time series
This experiment compares the performance and computational cost of the temporal and spectral implemen- tations of GVM for common time series and GP models of increasing complexity. We used real-world
11

---

Published in Transactions on Machine Learning Research (04/2023)
audio signal from the Free Spoken Digit Datase GVM was implemented with the metrics L1 and L2 both in the spectral and temporal domain to learn sample from the above dataset, which was 4300 samples long. The kernel considered was a spectral mixture (SM) (Wilson & Adams 2013), Fig. shows losses and running times as a function of the number of components the SM kernel; all runs use the Powell optimiser.
Spectral L
10a
10a
Spectral L2
Time [seconds]
Temporal L
Temporal Lz
Loss
10a
10-1
100
10
1 2 3
5 7 8 14 18 19 Number components
11111121111111117777111111
Number of components
Figure 4: Temporal and spectral implementation of GVM: Losses (Ieft) and running times (right) as a function of the order of the spectral mixture kernel.
From the left plot in Fig. (losses) let us recall that each implementation has its own metric and thus they not comparable directly, however, notice that both spectral implementations respect to the model order. Furthermore, the fact that the Temporal L2 loss increases with the model order suggests that the optimisation in the time domain is more challenging than in its spectral counterpart. In terms of computational complexily, the right plot in Fig. 4 confirms monotonicity of the computational cost with the number of kernel parameters, however, notice that the spectral implementations reached a plateau after 10 components; further suggesting the superiority of the spectral implementation in terms of ease of optimisation.
monotonic with
Lastly, Fig. 5 shows the fitting of both implementations and metrics in their respective domains. In the time domain (top plots) some peaks of the autocorrelation function (empirical kernel estimate K,), while the L2 metric (top right) aims to reach the peaks of Kn at the cost of missing central parts of it. In the spectral domain (bottom plots) we see only minor discrepancies for both metrics, with perhaps the most interesting feature being the difference the peak at around frequency 0.04, where Lı matched the peak and La provided wider fit.
see that the Li metric (top left) provides a generally acceptable fit but misses
Temporal distance:
Temporal distance: L2
1.0
1.a
Autocarrelation functtion
Autocomrlattio function
Leart 10-comp Spectral Mix.
Learnt 10-comp Spectral Mix.
Mm /
0.5
0.5
M wwwwamann
0.0
0.0
100
200
300
400
100
200
300
400
time
time
Spectral distance: L1
Spectral distance: L2
Periodogram hann)
Periodogram (w: hann)
Learnt 10-comp Spectral Mix.
Leamt 10-comp Spectral Mix.
100
100
Trequencies
frequencies
Figure 5: Example of matching for the temporal (vop) and saectral (bottom) implementation of GVM for the 10-component SM using the L1 (left) and L2 (right) metrics.
12

---

Published in Transactions on Machine Learning Research (04/2023)
With this comparison, validate the intuition that for kernels with concentrated spectral information (such as the spectral mixture) the spectral implementation of GVM is more robust to the metric and faces a less challenging optimisation task.
5.5
E5: Fitting a 20-component spectral mixture (different spectral metrics)
This experiment shows the effect of different spectral distances in the GVM estimates, also using an audio signal from the Free Spoken Digit Dataset (different from Experiment 4). Based on the promising perfor- mance of the spectral implementation of GVM, we trained a 20-component SM, kernel known to be difficult to Irain, and considered the spectral distances L1, L2, Wi and W2 (spectral); Itakura-Saito and KL were unstable and left our of the comparison. Fig. D shows the results of the GVM: observe that under almost all metrics, the 20-component spectral mixture matches the Periodogram (considered to be the ground truth PSD in this case). The exception is W2 which struggles to replicate the PSD peaks due to its objective of averaging mass horizonlally.
Spectral distance:
Spectral distance: L2
300
Periodogram (w: None)
300
Periodogram (w: None)
Leartt 20-camp Spectral Mix.
Learnt 20-comp Spectral Mix.
200
200
100
100
0.00
0.02
0.04
0.06
0.08
0.00
0.02
0.04
0.06
0.08
frequencies
frequencies
Spectral distance: W1
Spectral distance:
W2
Periodogram (w: None)
Periodogram None)
Leart 20-comp Spectral Mix.
Learnt 20-comp Spectral Mix.
200
200
100
100
0.00
0.04
0.05
0.08
0.00
0.02
0.06
0,08
frequencies
frequencies
Figure 6: GVM matching Periodogram with a 20-component SM under different spectral metrics.
5.6
E6: Learning spectral mixtures and parameter initialisation (L2, multiple orders)
4 components
8 components
12 components
16 components
not initialised BNSE-init SL-GP-init
not initialised
not initialised BNSE-init SL-GP-init
not initialised BNSE-init SL-GP-init
3K
BNSE-inkt SL-GP-init
3K
3K
2K
2K
2K
2K
-400
-400
400
-600
1K
-600
1K
600
1K
-600
1K
-000
0-
1450
1500
0-
1500
500
1000
1500
1000
1500
500
1000
1500
500
1000
1500
Figure 7: NLLs for the SE spectral mixtures (4,8,12 and 16 components) with different initialisation strate- gies.
In this experiment, GVM was implemented to find the initial conditions of GP with SM kernel (4, 12 and 16 components) to real-world 1800-point heart-rate signal from the MIT-BIH database We considered the L2 metric (spectral) minimisod with Powell and then passed the hyperparameters to an ML routine for 1500 iterations (using Adam with learning rate = 0.1). This methodology was compared against the random initialisation and provided by MOGPTK based on Bayesian nonparametric spectral
http:://ccg..ee mit. edu/timee seniees
13

---

Published in Transactions on Machine Learning Research (04/2023)
estimation (BNSE) Tobar 2018). Fig. 8 first shows the GVM approximations to the heart-rate PSD using 16-component spectral mixtures, for both both kernels.
Table 3: Computation times (secs) for fitting SMs.
4-comp
8-comp
12-comp
16-comp
GVM init
2.6
7.6
10.9
23.9
BNSE init
74.6
68.9
72.0
74.1
ML
420.3
461.6
529.2
582.6
Fig. shows NLL for the cases considered. Observe Chat: i) the non-initialised ML training becomes trapped in local minima in all four cases, ii) the initialisation provided by GVM provides dramatic reduction of the NLL, even wrt to the BNSE initialisation, iii) the "elbow" at the beginning of the GVM-initialisel case suggests that the ML training could have run for a few iterations (e.g,, 100) and still reach a sound solution. Table shows the execution times and reveals the superiority of GVM also in computational ional time.
Training time 23.18 [secs]
Training time 14.55 [secs]
Perlodo gram (W: None)
Porlodogram None)
400
400
Learnt 16 6-comp Spectral Mix.
Leart 16-comp Sinc Mix.
200
m m
0.0025
0.0075 0.0100 frequencies
0.0175 0.0200
0.0000
0,0050 0.0075 0.01.00 frequencies
Figure 8: GVM approximations of the PSD of a 1800-sample heart-rate signal using 16-components SE (top) and rectangular (bottom) mixtures. Training time shown above each plot.
5.7 E7: Learning isotropic SE kernel (standard variogram, 5-dimensional inputs, I.2)
Though our work focuses the single-inuutddmmnsioon case, the rationale behind the proposed GVM method is applicable to datasets of arbitrary input dimension, therefore, to motivate the use of GVM for multi- input GPs, we present minimal multi-input example. We considered 5-dimensional GP with SE kernel exp( = 2/
K(r) = o2 we assuned that all 5 dimensions had the same scale parameter, this is the usual setting the Variogram method Geostatistica (Cressie 1993 Chiles & Delfiner 1999). For each set of hyperparameters, we sampled 1000 points and then implemented GVM with the L2 (temporal) distance learn the GP. Table 4 the estimate error and standard deviation for four sets of hyperparameters averaged over 100 runs, from which the applicability of GVM to address the isotropic multi-input case can be confirmed. For illustration and resemblance to the standard variogram method, Fig. 9 shows the case o² = 5, 1 1= 1, and noiss = 1, where the learnt hyperparameters were 6² = 5.85, i = 1.09, and ánoise = 1.02.
1 [|+||?] + Se and considered four sets values for the hyperparameters. Notice that
Variogram for 5-dimensional isotropic SE kerne
Empirical variance
(GVM)
Lag magnitude
Figure 9: Empirical covariance and fitted covariance function via GVM. The data consisted of 1000 5- dimensional datapoints.
14

---

Published in Transactions on Machine Learning Research (04/2023)
Table 4: GVM applied multi-input data for four runs of synthetic data using different hyperparameters, 100 runs were executed to calculate the mean and
i.C., hoise- For cach set of hyperparameter, standard deviation of the estimates (i..,, 0², I, noise)
Run
a²
I
noose
0²
0 Doise
5
2
1
5.20 + 2.07
0.50 + 0.09
0.91 + 0.68
2
15
2
2
14.12
+ 13.14
2.07 + 1.09
2.04 + 1.51
3
3
3
0.5
2.68
= 3.33
3.24
+ 1.65
0.97
+ 0.71
0.5
5
0.5
0.45
= 0.62
3.96 = 2.50
0.56
0.30
6 Conclusions
By direct minimisation of the discrepancy among covariance functions, we have proposed novel methodology to pretrain stationary Gaussian processes, which avoids computation of the (cubic cost) likelihood function. The found hyperparameter values can then be passed to an MI.-based training routine as initial conditions to conclude the training of the model. Our approach, termed Generalised Variogram Method (GVM), represents a critical improvement in terms of computational complexity: we have shown, both theoretically and empirically, that for the particular case of the 2- Wassersstein spectral distance and location-scale PSDs, GVM is convex and its solution can be computed in single step. In experimental terms, we have shown the following properties of GVM in the general case: i) applicability to multi-input data, ii) stability wrt to different ways of computing the Periodogram, iii) consistency under different realisations of the GP unlike ML, iv) computational efficiency wrt ML and sparse GPss, v) a realistic alternative to compute initial conditions for ML resulting in considerable reduction of ML iterations, and lastly, vi) ability to train kernels of large number of components that are challenging to train from random initial conditions using MT.
In the general formulation of the proposed GVM (i.e., using either temporal or spectral divergence) the only requirement in our setup is stationarity, however, particular results in our proposal have specific requirements. First, when using spectral divergence (the main novel contribution of our work), it is needed that the kernels' Fourier transform (i.c., the PSD) can be computed and is Lebesgue integrable so that Dy So) in equation can be calculated; this condition is rather general and admits most stationary kernels used in the literature with the exception of the white noise variance see Sec/4.4. Second, for the exact case outlined in Thm. | we require that the PSD belongs a location-scale family, this includes standard covariance functions such as the (single-componnnt) spectral mixture, square exponential, sine, and cosine. PSDs that are not location-scale, such mixtures of the above kernels, can too be dealt. with spectral divergences, however, the solution is not exact and it has to be computed using numerical optimisation methods from equation Third, for all other covariances (including the white noise onc) we can use GVM with temporal divergences, which only requires the kernel to be stationary and evaluated pointwise, to find the solution via equation 7
Spectral kernels are a large class of covariance functions, therefore, the GVM is widely applicable in real- world scenarios in audio, seismology, astronomy, fault diagnosis, finance and any other fields where repetitive temporal patterns arise. In this sense, we hope that our work paves the way for further research in conceptual and applied fields. In theoretical terms, envision extensions towards non-stationary data using, e.g., time- frequency representations or mini-batches. In practical terms, we aim to address the general-input dimee--imension case and that our developed companion software will help others make use of GVM as and initialisation method for ML or to directly find the required hyperparameters.
15

---

Published in Transactions on Machine Learning Research (04/2023)
Acknowledgments
We would like to thank the anonymous referees for their valuable questions and comments, which allowed us to clarify the presentation of the paper and refine our results.
Part of this work was developed when Felipe Tobar was a Visiting Researcher at the Institut de Recherche en Informatique de Toulouse (IRIT).
We thank financial support from Google; ANR LabEx CIMI (grant ANR-LI-LABX-ONEE within the French State Programme "Investissements d' Avenir"; Fondecyt-Regular 1210606; ANID-Bashe Center for Mathemat- ical Modeling FB210005; ANID-Basal Advanced Center for Electrical and Electronic Enginecring FB0008: and CORFO/ CORFO/ANTD ANTD International Centers of Excellence Program 10CF1--9157 Inria Chile, Inria Challenge OcéanIA, STICAmSud EMISTRAL, CLIMATAmSud GreenAl and Inria associated team SusAIn.
References
S. Amari. Information Geometry and its Applications, volume 194. Springer, 2016.
F. Bachoc, F. Gamboa, J.M. Loubes, and N. Venet. A Gaussian process regression model for distribution inputs. IEEE Transactions on Information Theory, 64(10))6620-66377 2018.
M. Basseville. Distance measures for signal processing and pattern recognition. Signal processing, 18(4): 349 369, 1989.
E. Bernton, P. Jacob, M. Gerber, and C. Robert, On parameter estimation with the Wasserstein distance. Information and Inference: A Journal of the IMA, 8(4):657-676, 2019.
S. Bochner. Lectures on Fourier integrals, volume 42. Princeton University Press, 1959.
D. Burt, C.E. Rasmussen, and M. Van Der Wilk. Rates of convergence for sparse variational Gaussian process regression. In International Conference on Machine Learning, pp. 862-871. PMLR, 2019.
E. Cazelles, A. Robert, and F. Tobar. The Wasserstein- distance for stationary time series. 1EEE Transactions on Signal Processing, 69:709-721, 2021.
J. P. Chiles and P. Delfiner. Geostatistics, Modelling Spatial Uncertainty. Wiley, 1999.
C. Cortes, M. Mohri, and A. Rostamizadeh. Algorithms for learning kernels based on centered aligmment. Journal of Machine Learning Research, 13(28):795-82) 828, 2012.
N. Cressie. Statistics for Spotial Data. Wiley, 1993.
N. Cristianini, J. Shawe-
Shawe-Taylor, Taylor, A. Elisseeff, and J. Kandola. On kernel-target alignment. In Advances in Neural Information Processing Systems, volume 14. MIT Press, 2001.
A. Damianou and N.D. Lawrence, Deep Gaussian processes. In Proc. of the International Conference on Artificial Intelligence and Statistics, volume 31 PMLR, pp. 207-215, 2013.
T. de Wolff, A. Cuevas, and F. Tobar. MOGPTK: The multi-output Gaussian process toolkit. Neurocom- puting, 424:49 53, 2021.
T. Henderson and J. Solomon. Audio transport: A generalized portamento via optimal transport. anXiv preprint Xrx:::1906 00763, 2019.
J. Hensman, N. Fusi, and N.D. Lawrence. Gaussian processes for big data. In Proc. of the Conference on Uncertainty in Artificial Intelligence, pp. 282 290, 2013.
M. Hoffman and Y.A. Ma. Black-box variational inference as distilled Langevin dynamics. In Proc. of the International Conference on Machine Learning, pp. 9267-9277, 2020.
16

---

Published in Transactions on Machine Learning Research (04/2023)
F. Itakura. Analysis synthesis telephony based on the maximum likelihood method. In The 6th International Congress Acoustics, pp. 280-292, 1968.
V. Lalchand, W. Bruinsma, D.R. Burt, and C.E. Rasmussen. Sparse Gaussian process hyperparameters: Optimize or integrate? In Advances in Neurol Information Processing Systems, 2022.
S. Liu, X. Sun, P.J. Ramadgo, and R.P. Adams. Task-agnostic amortized inference of Gaussian process hyperparameters. In Advonces in Neural Informalion Processing Systems, volume 33, pp. 21440 21452. Curran Associates, Inc., 2020.
A. Mallasto and A. Feragen. Learning from uncertain curves: The 2-Wasserstein metric for Gaussian processes. In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
V. Masarotto, V. M Panaretos, and Y. Zemel. Procrustes metrics on covariance operators and optimal transportation of Gaussian processes. Sankhya A, 81 (1):1122-213, 2019.
G. Peyré and M. Cuturi. Computational optimal transport: With applications to data science. Foundations and Trends@ in Machine Learning, 11(5-6):35--677, 2019.
S. Popescu, D. Sharp, J. Colc, and B. Glocker. Hicrarchical Gaussian processes with Wassorstein-2 kernels, 2022. arXiv preprint. 2010.14877.
J. Candell and C.E. Rasmussen. unifying view of sparse approximate Gaussian process regres- sion. The Journal of Machine Learning Research,
2005.
C.E. Rasmussen and C.K.I. Williams. Gaussian Processes for Machine Learning. MIT Press, 2005.
G. Rios and F. Tobar. Learning non- non-Gaussian time series using the Box-Cox Gaussian process. In Interna- tional Joint Conference on Newral Networks, pp. 1-8, 2018.
Salimbeni and M. Deisenroth. Doubly stochastic variational inference for deep Gaussian processes, Aduances in Newral Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
A. Schuster. The periodogram of magnetic declination obtained from the records of the Greenwich Observatory during the years 1871-1895. Trans. Cambridge Philos. Soc, 18:107-135, 1900.
E. Snolson and Z. Ghahramani. Sparse Gaussian processes using pscudo-inputs. In Advances in Neural Information Processing Systems, volume 18. MIT Press, 2006.
P. Stoica and R. Moses. Spectral analysis of signals. Pearson Prentice Hall Upper Saddle River, NJ, 2005.
M. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Proceedings of the
Twelth International Conference on Artificial Intelligence and Statistics, volume 5, pp. 567-574, 2009.
F. Tobar. Bayesian nonparametric spectral estimation. In Advances in Neural Information Processing Systems 31, pp. 10148-10158, 2018
F. Tobar, T. Bui, and R. Turner. Learning stationary time series using Gaussian processes with nonpara- metric kernels. In Advances in Neural Informalion Processing Systems 28, pp. 3483 3491, 2015.
M. der Wilk, ST John, A. Artemev, and J. Hensman. Variational Gaussian process models without matrix inverses. In Proceedings of The 2nd Symposium on Advances in Appronimate Bayesian Inference, volume 118, pp. 9. PMLR, 08 Dec 2020.
M. Vetterli, Kovačević, and V.K. Goyal. Foundations of Signal Processing. Cambridge University Press, 2014.
C. Villani. Optimal transport: old and new, volume 338. Springer, 2009.
17

---

Published in Transactions on Machine Learning Research (04/2023)
K. Wang, G. Pleiss, J. Gardner, S. Tyree, K.Q. Weinberger, and A.G. Wilson. Exact Gaussian processes on million data points. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
A. Wilson and R. Adams. Gaussian process kernels for pattern discovery and extrapolation. In Proceedings of the 30th International Conference Machine Learning, volume 28, pp. 1067 1075. PMLR, 2013.
A Background: Fourier analysis of continuous-time stochastic processes
We present a brief description of the elements of Fourier analysis used in our proposal. For a more in-depth introduction of the subject, the reader is referred to (Stoica & Moses 2005 and (Vetterli et al. 2014).
Let us consider a stochastic process over RR denoted (yA)EER"; N. Since the usual condition of stationarity might be restrictive in some cases, we consider the following weaker version of stationarity.
Definition 3 (Wide-sense stationarity). The stochastic process y is wide-sense stationary (WSS) if its mean function is conslant and ils aulocorvelation only depends on the temporal difference. That is,
Bly+] p(t) =
(14)
E Elyr, Wta] = c(1i, (2) - -c(11- -2) 12)-
(15)
Though strict stationarity implies WSS, the implication does not hold in the opposite direction. However, for Gaussian processes (GP), whose distribution is fully determined by the first two moments, strict stationarity and WSS are equivalent conditions.
Definition 4 (Power spectral density). Let y be a WSS stochastic process with an absolutely integrable correlation funclion c(r) [E(y,y--)) -). The Fourier transform of c(.) given by
S(E) f c(T)e- J2mEE dr
(16)
is called power spectral density.
Observe that, if both and S satisfy the conditions for the inversion of the Fourier transform, then e and S are Fourier pairs, meaning that we also have
(()) dr)= f S(E)e de.
(17)
Equations 16 and 17 are a consoquence of the Wicner-Khinche Theorem (Vetterli et 2014 p. 292), which relates autocorrelation structure of WSS process with its distribution of energy across frequencies. This result is instrumental in the construction of GP: under the observation that GPs are uniquely determined by their autocorrelation (or covariance) function, the design of GP can be conveniently performed in the frequency domain by parametrising the PSD. Furthermore, recall that for zero-mean GPs the autocorrelation and autocovariance functions coincide; thus, we refer to the latter as the covariance kernel.
There are well-known pairs of kernels and PSDs that follow from the above observation. Figure 10 ratess five cases, showing the covariance kernels, their PSD and a sample of a GP with the corresponding kernel.
In practice, we regard data realisations of a stochastic process and we need to estimate the covariance or the PSDs from the available datasets. In the signal processing community, these quantities are usually estimated in a nonparametric fashion. For instance, the sample covariance in Definition | is an estimate of the covariance, which in the case evenly-spacod dataset ..a....,yNA.. takes the standard form
N k
0(kA) = N
HnAD(n+H) n=1
(18)
k
18

---

Published in Transactions on Machine Learning Research (04/2023)
Squaro exponential: Asg 0² exp (
PSDse = (M(Kse)(w)) = |ollye exp 2R
u(1) Ni00 Kee)
15000
10
-0.5
-90
Time
Frequenc
Time
Rational quadratic: Kaq= a 11+ 300
PSD mg [w) mumerical
M(1 Nio. Kao)
0.5
20
M
0
Proquency
Periodic: Koxa o² exp (R cos(unt))
PSDrnn F F[Koza] (woff cumerical
NEE) N R KPER)
-10
-5
10
0.5
0
0.5
11
-10
-5
10
15
Time
Frequency
Time
Cosine: Kes o² confort)
PSD 1 Row (w)) 2 1 3mm
((1) NO10, New
0.2
0.1
-10
10
-0.5
-15
-10
10
15
Time
Frequecc
Time
Spectral mixture: Kaar Yahe 2r1 (2m(2)))))
PSD
ΣM
exp
# (c mOP)
M(9) ('0, Kear)
Proquency
Figure 10: Illustration of the relationship between covariance and spectral representations of GPs. Left to right: kernel, PSD and a GP sample. Top to bottom: square exponential (SE), rational quadratic (RQ), periodic (Per), cosine (cos) and spectral mixture (SM) kernels. For the RQ and Per kernels, the PSD was computed numerically using the discrete time Fourier transform.
The problem of recovering the PSD from finite collection of observations ( Bty } referred as spec- tral estimation and the classic nonparametric method to perform this cstimation is called the Periodogram. This technique builds on the observation that, under mild assumptions, the definition of the PSD in equation
16 is equivalent to
27
S(w)
lim
yee just di
(19)
T
2T
Therefore, when only a finite dataset of observations is available, the limit and the expectation can be ignored in the above expression thus yielding the natural sample estimate of the PSD given by
S(w)
Σ
Vi jur
(20)
N
where N in the denominator replaced 2T as normaliser.
known as the Periodogram, a widely used method for estimating the PSD of WSS stochastic process introduced by Schuster (1900). Today, however, the concept of Periodogram refers wider class of techniques for spectral estimation that build the original formulation but at the same time address some of the known drawbacks related to the biasedness and large variance of the Periodogram.
Equation 20
Two widely used extensions of the Periodogram are the Bartlett and the Welch methods. The former operates by splitting the data into segments and computing the standard Periodogram over each one of them, to then average over all computed Periodograms with the aim reduce the noise the estimate. The latter follows the same concept but also multiplies each segment by a window so as to mitigate the effect of the border of the segments in the estimation. Usual choices for the windows are the Hann and Hamming functions.
19

---

Published in Transactions on Machine Learning Research (04/2023)
B Definition of distances and divergences considered
For two functions fi and f2, we have the following general distances:
1-Fuclidean TA(i, f2) = ((((())) - f2(5))d5
2-Fuclidean [2(S1: (Q) [2(E))
Furthermore, when /i and f2 are densities with quantile functions Q and Q2 respectively, have the additional divergences:
1- Wasserstein Villani 2009
Peyré & Cuturi
2019
/ (2(())-2((P))) (p())dpp
Wi(f1, f2)
2- Wasserstein
Villani
2009
Peyré & Cuturi
2019
Wa(fi, [ (Qi(p) d2(p))
Kullback-Leibler:
(fill)))))= log f1(5)
fi(E)
f2((E)
Itakura-Saito
Itakura 1968
Dis(A1lfa) /
(116 log fI(E) (2(e)
[i(9)
1) de [2(E)
Bregman divergences convex,
Amaai
2016
for function
IR + R that is differentiable and strictly
Dclfi, G(fi) (VG())) f2) f2)
Here, we have assumed that both fi and fz integrate unity, in the cases where this condition is not met, the densities can be normalised before computing the distance.
C Proofs
C.1 Convexity of spectral loss for W2 and location-scale family
Proof of Theorem
We recall that
W Sp.o) - d(p))dp.
(21)
From direct application of the rule of differentiation under the integral sign, we obtain the gradient for the location-scale family:
V Vn W3, 2 J' (Qno(p) (o
(22)
= / - ( + (p)
1 2 /' (n - Q(p)) (go160) dp.
20

---

Published in Transactions on Machine Learning Research (04/2023)
The Hessian for the location-scale family:
Qp,
H H..o (5,
(Q.1(H) Qua(1)3 dp.
(23)
The determinant of the Hessian (via Jensen's inequality):
H H/22- [/ Qo.(1))dp (( > Qo.1())ppppp "
do 0.
(24)
where the inequality is strict due to the strict convexity of (-)².
Therefore, the first order conditions are given by:
L (1 = == - - =
(Q(p) dQo1()) / Q(D)dp
(25)
and ((u+ (wt
(Qu, - (
(26)
I (Q(p) - (o)) So (p,
fo Q(p)Qo.
(27)
So Q%,
where the last expression we have used the fact that the location of the prototype So,1 is zero and so is its mean, meaning that if x S we can write So pQo. pEx~ [x] =
]
C.2 Learning from data (W₂ distance and location-scale family)
Proof of Proposition First, recall that in the location-scale family 0 = (u, We denote by Q and Q.. the respective quantile functions of 8 and S,. Then, following the solutions in Theorem and Jensen's inequality we can compute the following upper bound for the location parameter μ²:
(Bl 11,1)2 & BBIp" - 12112 =Fx== Qw))dp fooddp-f"ee
S. Qn Co
< E L// 1Q(p) - Qn(())) Es,
and by hypothesis S.) + 0. Similarly, now using Hölder's inequality, we obtain the following bound for the scale parameter o":
Elo" - o = [ Q()) S. (..
/ (((()) Qo,
< E
< (// (Q(p) - @n(o)))a)) $ $]
( 10. (2).()))22))
E(W2(S, Sw.)) (f. (0. (Q..)69)
which tends to 0 by again Jensen's is inequality, (.)))) < R -
21

---

Published in Transactions on Machine Learning Research (04/2023)
D
Code
We have developed short self-contained toolbox which is available The purpose of the code is to facilitate the use of the proposed method by the community and in particular to replicate all our results. For the reader's convenience, Jupyter Notebook ExpO_minimad wwm.miinmmmaa.....m. _ex___________ ex ipynb showing minimal example of the toolbox implementation presented here.
Minimal working example of provided code
[1]: #general imports import numpy as np import batplotlib.pyplo. as plt #our package from waflgp import
import utils
[2] #load data signal np loadtxt (aaattxx("Daa/h12.... DDData/h22.txt
[3]
#ins iinstantiate model, sum of 16 Gaussians
9 16
suagaaeaouut aim 'learning', kernel 'qSM') #Spectral Niz #set frequencies (optional) freqs np. gp pp.set_freee (oe_frrgssffeqq
gP
[4)
#train with pertodogram, L2 metric and components gP. wewtttrrannaeeaeee (method 'periodogram', metric 'L2,, order=q) #plot Periodogram and best PSD fét pe.poo_ppddtitle f Minimal example')
Optimization terminated successfully.
Current function value: 0. 0.0028
Iterations: 33
Function evaluations: 21261
L2-ok
Minimal example
Periodogram (w: None)
400
Learnt 16-comp Spectral Mix.
200
0
0.0000
0.0025
0.0050
0.0075
0.0100
0.0i25
0.0150
0.0175 0.0200
frequencies
22

---

Published in Transactions on Machine Learning Research (04/2023)
E
Additional figures for E2
Experiment E2 showed the scnsibility of GVM to the choice of Periodogram method and window for two presents all the figures corresponding to the estimates in Table for the Exp-cos
kernels. Figures II and 12
and Sine kernel respectively.
Learnt parameters are loc: 0.0501, scale: 0.0089
Learnt parameters are loc: 0.0515, scale: 0.0076
Perlodogram Nonc)
Porlodogram (w:
Leart kemel
Learnt kemee
100
Ground truth Exp cos kornel
100
Ground truth Exp-cos cos kornel
1/6
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.00
0.01
0.02
0.03
0.04
0.05
0.07
frequencies
frequencies
Leart parameters are loc: 0.0514, scale: 0.0077
Learnt parameters are loc: 0.0500, scale: 0.0113
Periodogram hamming)
None, Learnt kernel
Leart kemel
40
100
Ground truth Exp-cos kernel
Ground truth Exp-cos kemel
20
0.00
0.04
0.05
0.06
0.07
0.04 frequencies
0.08
0.10
frequencies
Learnt parameters are loc: 0.0495, scale: 0.0101
Learnt parameters are loc: 0.0495, scale: 0.0099
Bartlett
40
Bartlett hamming, 10)
Learnt Kernel
Learnt ker nel
Ground Exp-cos kernel
Ground truth Exp-cos kemel
frequencies
frequencies
Learnt parameters are loc: 0.0504, scale: 0.0100
Learnt parameters are loc: 0.0506, scale: 0.0093
Welch None, bins:
(w: 10)
40
Learnt kernel
Learnt kernel
Exp-cos
Exp-cos
20
frequencies
frequencies
Learnt parameters are loc: 0.0506, scale: 0.0091
Welch (w: hamming, 10)
Learnt kernel
Ground Exp-cot cos
0.03 frequencies
0.07
Figure 11: Exp-cos kernel.
23

---

Published in Transactions on Machine Learning Research (04/2023)
Learnt parameters are loc: 0.0506, scale: 0.0108
Learnt parameters are loc: 0,0505, scale: 0.0071
Perledogram Nonc)
Parlodogram (w:
Leart kemel
200
Learnt kemee
Ground truth Sinc kemel
Ground truth Sinc kerne
100
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.00
0.01
0.02
0.03
0.04
frequencies
frequencies
Learnt parameters are loc: 0.0505, scale: 0.0071
Learnt parameters are loc: 0.0506, scale: 0.0244
100
Periodogram (w: hamminng
Nono, bins:
Leart keme!
Learnt kemel
Ground truth Sinc kemei
Ground truth Sinc kernel
100
0.00
0.01
0.02
0.03
0.04
0.05
0.00
0.02
0.04
0.06
0.08
0.10
0.12
frequencies
frequencies
Learnt parameters are loc: 0.0504, scale: 0.0210
Learnt parameters are loc: 0.0504, scale: 0.0195
100
100
Bartlett (w: hann, bins: 10)
Bartlett (w: hammning bins: 10)
Leart keme!
Learnt kernel
Ground truth Sinc kemel
Ground truth Sinc kernel
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.00
0.01
0.02
0.03
0.04
0.05
0.06
frequencies
fre frequencies
Learnt parameters are loc: 0.0502, scale: 0.0168
Learnt parameters are loc: 0,0505, scale: 0.0135
100
Welch (w: None, bins: 10)
Welch (w: hann, bins: 10)
kemel
kemne
Ground truth Sinc keme!
Ground truth Sinc kernel
0.00
0.02
0.04
0.06
0.08
0.00
0.01
0.02
0.03 frequencies
0.04
0.05
0.06
frequencies
Learnt parameters are loc: 0,0505, scale: 0.0127
100
Welch (w: hamming, bins: 10)
50
Ground truth Sinc kernel
0
0.00
0.01
0.02
0.03 frequencies
0.04
0.05
0.06
Figure 12: Sine kernel.
24


Reviews:
[{"review_id": "DwcLuxDn1B", "reviewer": "Reviewer_Z3Qg", "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\nIn my opinion, the paper tackles an interesting idea for a well-known problem in the GP community. There has been always a general interest and intuition around the idea that spectral methods or signal-processing inspired tools could solve many of the issues that GP learning usually faces. In this case, I enjoyed the spirit behind the connection between the expected log-marginal likelihood and the divergences, since it opens new doors to the problem of learning hyperparameters. The general idea of bringing PSD/Kernel distances to be used as a target loss to minimise in order to reach meaning results makes sense to me as well. Some of the experimental results on toy data also seem promising.\n\n** Weaknesses **\n\nEven while I enjoyed reading the paper, there are points or ideas that are not entirely explained or extremely clarified. Some other details or experimental results seem somehow limited to me, in the sense that they do not compare with other SOTA methods who have also tried to solve this problem (see references below).\n\n- One example, for instance is the connection between the KL divergence and the expected log-marginal likelihood. The claim that both quantities are equivalent is not entirely true to me. First of all, there are terms weighted in a different manner, and the quantities are not the same at all. Also, the expected log-marginal likelihood is not the quantity of interest for us from a Bayesian perspective, right? Since we are observing a single training dataset. Later on, the use of divergences between covariance matrices seems odd to me, or at least I would expect this sort of metrics being computed wrt probability densities, not covariances instead.\n- Another weak point is Section 3 and 4, where there is a bit of lack of clarity. From the reader's perspective, one feels that the thread of the paper fuses too many concepts from different sides, but not really contributing to the central idea of the paper that is the learning of hyperparameters. For instance, the presentation of Algorithm 1 and Algorithm 2 is somehow very light, and not very reproducible from the perspective of a computer scientist. How are exactly this losses computed, how is the iteration through the data on this algorithm? We would need answers like that.\n- On the experimental side, I find that the results are promising, but somehow preliminary to me. As a first view/intuition of the performance of the method is right, but to really be a mind-changing method for the GP community, one should include rigorous comparisons with other methods in the literature. Also playing with larger data, different datasets (i.e. UCI datasets at least) and other well-known benchmarks... that would be the most convenient.\n- Finally, from the contribution aspect, one might easily feel that the proposed solution are L1,L2, and Wasserstein distances between estimated covariances and PSD functions, which itself is fine, but I feel that some of them could have been partially considered in the very long literature of GPs since its boom in the ML community after 2006. At least, a larger review would be required to be entirely trusted, in my opinion.\n\n\n** References **\n\n1] V. Lalchand et al. \"Sparse Gaussian Process Hyperparameters: Optimize or Integrate?\" NeurIPS 2022\n\n2] S. Liu et al. \"Task-Agnostic Amortized Inference of Gaussian Process Hyperparameters\" NeurIPS 2020 (see related work in this paper and references to approaches for learning GP hyperparameters).\n\n3] D. Burt et al. \"Rates of convergence for sparse variational gaussian process regression\" ICML 2019 (I think this one is related to what is proposed in the paper, even if it focuses on sparse methods. I say this bc authors here use the Nystrom approximation to the data covariance matrix and a bound on the KL divergence between the true posterior and the approximated one.)\n\n4] Rassmussen et al. \"Gaussian Processes for Machine Learning\" 2006 - Ch. 5 including the references to cross validation and other related methods.\n", "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": "Inspired in signal processing/geostatistics methods, the paper proposed a new technique for learning the kernel hyperparameters in Gaussian process (GP) models. The key idea behind is to avoid the evaluation/computation of the likelihood function. To do that, the paper first makes an observation on the similarity between the (expected) log-marginal likelihood  of the GP and the negative KL divergence between two multivariate Gaussian distributions. In particular, this leads to consider distance/divergence based methods between the kernel-derived covariance matrix and the true empirical one. The focus is then placed into both distance based i.e. L1, L2 norms and spectral metrics (between the PSD of the true and fitted covariance matrices). In practice, estimating both leads to methods that are promising when learning the hyperparameters, both in computational cost and performance (reaching the true ones).\n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "To me the paper has three main points that should be fixed or at least improved:\n\n- Establish new connections between different approaches is generally difficult, and usually makes one to incur in odd mix of notations, concepts, ideas, etc. In this case, the task would be between the general approach of learning GP hyperparameters, spectral methods in signal processing and divergences/distances. In my opinion, this could be improved significantly, mainly in Sections 3 and 4, where the goal of learning hyperparameters is somehow secondary behind the presentation of several techniques and ideas.\n\n- Related work and a deep review of the papers working on the learning of GP's hyperparameters -- this at the end, should be also accompanied by extra experiments, proving the similarities and dissimilarities between the proposed method and previous efforts in SOTA. On this, I want also to remark that the results seem promising to be, but not entirely well presented in the GP literature framework.\n\n- I would like to see the presentation of the computational cost, algorithmic solution better presented or at least clearer for the authors. If it is computationally cheap and fast, one should feel that is easy to code or implement.", "broader_impact_concerns": "I do not detect particular ethical implications of the work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "JuxuX54Ftw", "reviewer": "Reviewer_dRzy", "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 presentation is reasonably clear although it does presume quite a lot of material (e.g. several kernels that are not defined such as sinc). Given the page limit this seems difficult to address. The bigger challenge I had with the exposition was keeping track of what theoretical and experimental results apply to what approaches. I suspect there is some nice summary that could be provided: maybe a table where the rows are different approaches and the columns are different properties (e.g. convexity)? \n\nOther than the minor concerns below I'm positive about the paper: good initialization is one of the important but overlooked aspects of statistical ML, and the theory and empirical results here should both be useful to the field. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors consider how to initialize the hyperparameters of (time-series) GPs using ideas related to but extending the variogram. Both theoretical results and practical implementations are presented, first for a more restricted setting (location-scale power spectral densities, PSD [an unfortunate acronym!] with 2-Wassertein distance) and then for more general setups. A nice series of experiments of increasing complexity are conducted, from a simple squared exponential kernel initialized via the variogram to complex spectral mixtures. \n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I would have liked some more clarity around what results apply to multi-input settings (and how). The intro and math setup is very focused on time-series analysis, but the experiments section then starts out with E1 being a multi(=5) input example! My understanding is that for E1 the lengthscale is shared across the input dimensions: can this be relaxed to one per dimension? In that setting smart initialization would seem particularly important because there are potentially many hyperparameters. \n\nFor E1 and E3 the experiments should vary the true hypers (at least l and mu, maybe sigma) over some range – otherwise for all we know GVM always returns a number around 1 in the E1 setting which just happens to be the \"correct\" answer. Then the plot of true vs estimated hyperparameter could be inspected. \n\nI've annotated minor comments on the pdf here: \nhttps://www.dropbox.com/s/shoogfeux3th8gs/computationally_efficient_init.pdf?dl=0", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "i5GrldEcQ7", "reviewer": "Reviewer_VErd", "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": "### Strength: \nIt’s a novel and interesting idea. The proposed method provides an alternative for learning GP hyperparameters. The main advantage of the method is the scalability (linear time complexity) and stableness to random initialization. \n\nThe work has both theoretical and empirical contributions. The experiments are organized with questions, making it easier to consume.\n\n### Weakness: \n\nA **critical** point is that I didn’t see any results using eq(4) (please correct me if I am wrong). To show learning GP hyperparameters in the spectral domain helps, it is critical to compare the proposed method to that baseline.\n\nClarity and self contained writing. \nI find it difficult to build the context when reading the paper, partially due to the fact that I am less familiar with Fourier Transform and not using them much in my own work. But I think the authors could also make the paper more self contained for the readers in a slightly different area. Below I give some examples:\n\n- In Theorem 1 where the PSD is restricted to location scale family and using W2 distance, the optimal location and scale parameter can be found in one step. But how to get back to the original kernel hyperparameter (I assume by using inverse Fourier transform)? Stating it explicitly, even better with some examples, could help. \n\n- In the end of Proposition 2, the author mentioned “windowing and the Welch/Bartlett technique” and they are used in Table 1, but they are never explained in the text. \n- In Section 4.1, authors briefly state how to handle general cases but in my opinion, they would be better presented in pseudo code with clear input and output.\n\n- In Section 5.5, how to get the groundtruth Periodogram for the Free Spoken Digit Dataset?\n\n\nLimitation: The method only supports PSD with a specific parametric family and explicit inverse Fourier transform. This suggests, as stated also by the authors in the end of Section 3.1, a set of 6 kernels (Dirac delta, Cosine, Square Exponential (SE), Student’s t, Sinc, Rectangle) and their mixtures can be covered. Is this a complete list? The author should discuss this and its impact on real-world applications explicitly.\n\nTime series data. The experiments are mostly 1-d input time series data. The only exception is in Section 5.1, but it is a synthetic one. It would be nice to have more real-world multi-dimensional input regression tasks. If the author focuses on time series data, that should be stated explicitly in the contribution. \n\n Also, there is a whole area of kernel alignment and it can be seens as a variant of eq(4). A good pointer for kernel alignment is:\n\nCortes, Corinna, Mehryar Mohri, and Afshin Rostamizadeh. \"Algorithms for learning kernels based on centered alignment.\" The Journal of Machine Learning Research 13 (2012): 795-828.\n \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work starts with the observation that maximizing likelihood for GP is equivalent to minimizing KL-divergence between the true GP and the GP to be learnt. Then it looks for methods to find GP hyperparameters by minimizing a loss function on 1) sample covariance and 2) its Fourier transform through power spectral density (PSD). \n\nThe author proved, for location scale PSD and W2 distance, one can find the optimal location and scale parameters analytically (Theorem 1) and it converges to the true minimizer (Proposition 1) and also in more general cases under certain conditions (Proposition 2).\n\nThe empirical results show that the proposed method can recover the true hyperparameter with good accuracy (Sec. 5.1), stable with initializations (Sec. 5.2), sensibility (Sec. 5.3), scalability (Sec. 5.4) and kernel mixtures on real-world time series data (Sec. 5.5 and 5.6). \n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I would suggest the authors:\n- Adding results using eq(4) and comparing to the PSD based methods.\n- Making the paper more self-contained, see some examples in the weakness section.\n- Explicitly discussing the focus and the limitation of the methods. Whether it focuses on time series (1d) data and what kernels can not be supported currently and what is the implication of that.\n- Including kernel alignment as variants for eq(4) and even better, providing empirical results.\n", "broader_impact_concerns": "No.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All reviewers agreed on the technical novelty of the paper. The reviewers were also happy with the updates made by the authors and the additional supportive evidence they provided. This work is of interest to the community so I am supporting acceptance as is -- this is the third revision provided by the authors, so the paper is pretty polished.
Decision: Accept (poster)