Integrated Variational Fourier Features for Fast Spatial Modelling with Gaussian Processes
Authors: Talay M Cheema
Carl Edward Rasmussen
Venue: Accepted by TMLR ()

Abstract: Sparse variational approximations are popular methods for scaling up inference and learning in Gaussian processes to larger datasets. For $N$ training points, exact inference has $O(N^3)$ cost; with $M \ll N$ features, state of the art sparse variational methods have $O(NM^2)$ cost. Recently, methods have been proposed using more sophisticated features; these promise $O(M^3)$ cost, with good performance in low dimensional tasks such as spatial modelling, but they only work with a very limited class of kernels, excluding some of the most commonly used. In this work, we propose integrated Fourier features, which extends these performance benefits to a very broad class of stationary covariance functions. We motivate the method and choice of parameters from a convergence analysis and empirical exploration, and show practical speedup in synthetic and real world spatial regression tasks.

Published in Transactions on Machine Learning Research (04/2024)
Integrated Variational Modelling with Gaussian Processes
Fourier Features for Fast Spatial
Talay M Cheema
tmc/90cam.ac.
Deportment of Engineering University of Cambridge
Carl Edward Rasmussen
Department of Engineering University of Cambridge
Reviewed on OpenReview: https: // openreview. net/ net/forum? formm? id= PdPPtWaaa
Abstract
Sparse variational approximations are popular methods for scaling up inference and learning in Gaussian processes to larger datasets. For N training points, exact inference has O(N3) cost; with M & N features, state of the art sparse variational methods have O(N cost. Recently, methods have been proposed using more sophisticated features; these promise O(M)) cost, with good performance low dimensional tasks such as spatial modelling, but they only work with a very limited class of kernels, excluding some of the most commonly used. In this work, we propose integrated Fourier features, which extends these performance benefits to very broad class of stationary covariance functions. We motivate the method and choice of parameters from a convergence analysis and empirical exploration, and show pract ical speedup in synthetic and real world spatial regression tasks.
1 Introduction
Gaussian processes (GPs) are probabilistic models for functions widely used in machine learning applications where predictive uncertainties are important for example, in active learning, Bayesian optimisation, for risk-aware forecasts [Rasmussen & Williams 2006 Hennig et al. 2022 Garnett 2023). The hyperparameters of those models are often learnt by maximising the marginal likelihood, so is important that this quantity can be evaluated fairly cheaply, especially to facilitate comparison multiple models, or multiple random restarts for robustness. Yet, for N datapoints, the time cost is O(N3), associated with calculating the precision matrix of the data and its determinant. This prohibitively large for many datasets of practical interest. A particularly important class of problems is spatial modelling, where often Gaussian process regression is applied in low (2-4) dimensions to very large datasets, and ideally the choice of prior process, encapsulated by the prior covariance function, is guided by domain-specific knowledge.
One popular approach for improving scalability to use a sparse variational approximation, wherein M < inducing features are used as compact representation, and a lower bound of the log marginal likelihood is maximised. In the conjugate setting, where the measurement model affine with additive, white, and Gaussian noise, the variationally optimal distribution of the inducing features available in closed form (SGPR: Titsias 2009). This reduces the size of the precision matrix from N N to M x M, but in practice multiplications by the cross covariance matrix between the data and the features dominates the computation cost at M²)
SGPR with inducing points where the features are point evaluations of the latent function - can interpreted as replacing data measured with independent and identical noise with a reduced pscudo-datasct measured at different locations and with correlated and variable noise. This pseudo-dataset selected to minimise distortion in the posterior distribution over functions (or, equivalently, minimum discrepancy

---

Published in Transactions on Machine Learning Research (04/2024)
between the lower bound and the log marginal likelihood). This is particularly advantageous in the common case that the data is oversampled, where it is possible to set M < N with asymptotically vanishing distortion (Burt et al. 2019 2020a). Inducing points are the state of the art solution, but the scaling with N is problematic. One popular way to avoid this is to use batches of data (SVGP; Hensman et al., 2015). But this necessitates the use of stochastic, lypically first order, optimisers; in the conjugate setting, this leads to iterative learning of the variational distribution which is otherwise available in closed form.
Ideally, would like to find approximate inference method which avoids the O(N) scaling for any dataset reasonably
and any prior by careful design of the inducing features. But this is more generality than expect, and methods are generally restricted in the prior covariance functions they support, which in turn restricts the freedom of modellers. Existing work on zonal kernels on spherical domains and tensor products of Matérn kernels on rectangular subsets of RD give a recipe for taking the O(N) part of the computation outside of the optimisation loop for low dimensional datascts (Dutordoir et al. 2020 Hensman et al. 2017
Cunningham et al., 2023).
In this work, we propose integrated variational Fourier features (IFF), which provide the same computational benefits, but for much broader class of sufficiently regular stationary kernels on RD | We achieve this by allowing for modest numerical approxinations in the evaluation of the learning objective and posterior predictives, rather than searching for mathematically exact methods. Yet, in contrast to those previous approaches, we also provide convergence guarantees. This both provides reassurance that the number of features required scales well with the size the training data, and shows that the numerical approximations do not significantly compromise performance.
In Scction
review variational GP regression in the conjugate sctting, and we review related work in Section 3 In Section and guidance for tunable parameter selection follows in Section 4.1 Finally in Section S we evaluate our method experimentally, showing significant speodup relative to SGPR in low dimensions, and competitive performance compared to other fast methods, with broader applicability. summary of our contributions is as follows.
f we present our IFF method, and the complexity analysis; the main convergence results
We present a new set of variational features for Gaussian process regression, whose O(M) mem- ory cost and O(M)) per optimisation step computational cost greatly increases scalability for low dimensional problems compared to standard approaches demonstrated on large scale regression datasets and can be applied using broad class of stationary covariance functions on RD.
We provide converge results demonstrating the number of features required for an arbit good approximation to the log marginal likelihood grows sublinearly for a broad class of covariance func- tions.
We provide reasonable default choices of parameters in our algorithm, including the number of inducing features M, based on an empirical study and motivated by our theoretical results.
2 Background
In the conjugate setting, the probabilistic model for Gaussian process regression is
f~ f~ Un = f(xn) Pr Pn ~((,00
(1)
for n E N1;N), with each I, E RD and Per Un E R, with the covariance function or kernel k RD x RD + R symmetric and positive definite. Let f = [..., f=,,),...] ...]T, and Kạt be the covariance matrix of finite dimensional random variables a, b. For example, Kÿt is the N x N matrix with [Ki]]an' k(2,, 2.). The postcrior predictive at some collection of inputs x. and marginal likelihood are as follows, where L1.N.y = VI,N, and K, is defined analagously to KBp. Then, the posterior predictive distribution and
We assume the kernel's spectral density has bounded second derivative.
2

---

Published in Transactions on Machine Learning Research (04/2024)
marginal likelihood are follows
Rasmussen & Williams
2006
Chapter 2).
p(f(z.))z, v) = (f((z)) + o K., K,j(Ki + -011)-1Kp p(ylz) (V(00,, Ky I) L -=logplyla
K.,
(2)
We optimise the latter with respect the covariance function's parameters. The data precision matrix A (Kş + o which depends the value of the hyperparamoters, dominates the cost, as for cach evaluation of the log marginal likelihood L, need to compute its log determinant and the quadratic form y Ay, both of which incur O(N²) computational cost in general. Note that the posterior predictive is the prior process conditioned on = f p.
For the variational approximation, we construct an approximate posteriof ((f) = f p(fl((((()((())) = p(fly) where u #1:M is a collection of inducing features with prior distribution p(u). That is, we condition on u instead of y and average over an optimised distribution on u. The classic choice inducing points, where Um (=m) for some Zm IR We maximise lower bound the log marginal likelihood (DKL is the KT divergence).
F f= a( ply.f.u) {()
- df C-Dr----(l)))
(4)
More generally, um chosen to be a linear functional om of [Lázaro-Gredilla & Figueiras- 2009). denoted (dm.f) C with associated parameter 2m; in order that u is Gaussian priori. For features other than inducing points, these are termed inter-domain features. Let on be such that (0m: f) (dm: f)* (the complex conjugate) and let I. be the covariance operator corresponding to k. That is, (x.) (o.,,k
k(z,, ))> Then (Bogachev
1998. Chapter 2; Lifshits 2012
(dm: f ~ N(0, (pm, KΦm))
E Mo , f3 = (dm: Kdm!)
(x.) f`] = (d'm k(Xos
and for convenience define
(m(xx)) = of (2 xx) def (o b(Ixs -)) c² (Ix, 2m) = (6mm k6,,,1,)) km.m. = k(zms 2mm) def (0m, KDm)) ) = (@m c(Zmm )) = (5m c(,, zm))
=
(5)
(6)
which give the entries of the covariance matrices Kaf and Kan- With indueing points, c k k. But, more generally, p(u) NO, Kau), and the optimal q(u) available in closed form
q(u) N(uu E.) with E, K (Kuu - Ku
Pu =
" K Knçy
(7)
with corresponding training objective
Titsias 2009
1 - -
Fipm, En) = log Ny/0, K Kut K Kuf +211)
Kug K- Kuf)
(8)
2
wherein the structured approximation to the data precision A' kKi + 021)-1 However, by exploiting standard linear algebra results (Appendix A), the inverse and log determinant be isolated to B = (Km-to "K, (which is the precision of the appropriately noise-corrupted features u+ Kusp) and K-1 both of which are only M x M. However, in practice, the dominant cost is O(NM²) to form Kup K since generally M < N and the cross-covariance matrix depends nonlincarly must be recalculated each lime Burt et al.) (2020b). Put differently, the features (.) are dependent on the hyperparameters. By choosing the linear functionals carefully, we aspire to find features which do not depend on the hyperparameters, so that K can be precomputed and stored, reducing the cost to O(M3), without compromising on feature efficiency.
the hyperparametors,
In standard minor abuse notation, we write the distributions over densities, though none exist.
3

---

Published in Transactions on Machine Learning Research (04/2024)
The posterior predictive at new points n, is calculated
(1, / p(((s.))a,,,, du
( Ki Ku= K K Kus + Ki.KK Ky En K..).
(9)
Moreover, high quality posterior samples can be efficiently generated (for example, when the number of inputs in I, is very large) by updating a random projection approximation of the prior (for example, using random Fourier features) using samples of the inducing variables (Wilson et al. 2020).
We note that the lower bound property of this training objective makes meaningful: increases in L involve either increasing the marginal likelihood with respect to the hyperparameter, or reducing the KT. divergence from the approximate posterior to the true posterior. This KL divergence is between the approximate and true posterior processes, giving reassurance on the quality of posterior predictive distributions (Matthews et al. 2016). Moreover, the inducing values act as meaningful summary of the training data which can be used in downstream tasks, for example in order to make fast predictions (Wilson et al. 2020
3 Related Work
There are two other main, broadly applicable, approaches to reducing the cost of learning, which are com- plementary:
using iterative methods based on fast matrix-vector multiplications (MVMs) to approximate the linear solve and log determinant, and
2. directly forming a low-rank approximation to the kernel matrix Kjt-
In the former case, the cost is reduced to O(N²) in exchange for modest error, since only limited number of steps of the iterative methods are needed to get close to convergence in practice. This is particularly advantageous when performing operations on GPU (Gardner et al.] 2018b Pleiss et al. 2018). and when A has some special structure that permits further reductions due either to structure in the data or k (Saalçi 2011 Cunningham et al. 2008). Direct approximations of Ktt include by projections onto the Fourier basis for stationary kernels (Random Fourier Features, RFF (Rahimi & Recht 2007) and variants [Lázaro-Gredilla et al. 2010: Gal & Turner 2015), interpolating from regular grid points (stochastic kernel interpolation (Wilson & Nickisch 2015] Gardner et al.] 2018a, SKI)), or projecting onto the highest variance harmonics on compact sets (Solin & Särkkä 2020). Notably, SKI makes use of fast MVMs with structured matrices to obtain costs which are linear in N for low D. In contrast to the variational approach, these methods tend lo approximate the posterior indirectly and the approximations may be qualitatively different
to the exact posterior (see, for example, Hensman et al. (2017)). ). Variational methods can also be viewed making the low rank approximation Ky K Kuf to the kernel matrix, but note that the training objective differs from simply plugging in this approximation to the marginal likelihood, it has an additional trace term (Equation (8)). Recently, authors have attempted to incorporate nearest neighbour approximations into a variational framework (Tran et al. 2021 Wu et al. 2022).
One notable approach which does not into these categories is using Kalman filtering: Gaussian process regression can be viewed as solving a linear stochastic differential equation, which has O(N) cost given the linear transition parameters (Särkkä et al. 2013).. In practice, if the data does not have additional structure such regularly spaced inputs, computing these transition parameters will dominate the cost.
Finally, by careful design of the prior, we can create classes of covariance function for which inference and learning are computationally cheaper (Cohen et al.) 2022 Jø rgensen & Osborne 2022). However, these are not broadly applicable in the sense that the classes of covariance function (and hence the prior assumptions) are limited, and only suitable to certain applications.
Fourier features If we restrict the prior be stationary, that is k(æ,x') = k(x x') K(T), then k has unique spectral measure. We assume throughout that the spectral measure has a proper density
4

---

Published in Transactions on Machine Learning Research (04/2024)
sE) = I ( & dr. Note that according to the convention we use here, s de k(0)
features is Fourier features, appropriately normalised by the spectral density: (01.;.) = f()) 12 /s(E) dr. These are independent with unbounded variance
A first attempt at byperparameterinddpeedde
Lázaro-Gredilla & Figueiras- Vidal. 2009 Lilshits 2012 Chapter 3), which can be shown as follows.
(1z1,))= = / Mxx',x) -2m ==t"== /s(() dr = e -22755
(10)
(5.55) = f cla' (g)e=2ne 2 of 1s(5))
= f f (=2xx"T(-e)
- (()/s(()
(11)
Here, 5 is the Dirac delta. These features are unsuitable for constructing the conditional prior p(flu) informally, the prior feature precision K vanishes but Kaj is finite, so the conditional prior mean Ki, Kuu u is zero (Figures la and lb
Yet the general form is promising since (i) K fu depends only on the chosen features and the training inputs F, so indeed if we can fix the frequencies to good values, then can precompute this term outside of the optimisation loop, and reduce the cost of computing A' to ((M3) per step (Appendix A and also (ii) the features are independent, so Kun. would be diagonal.
Modifications to Fourier features include applying a Gaussian window (rrzilo-Gredilla & Figuciras 2009) which gives finite variance but highly co-dependent features, and Variational Orthogonal Features, where (dm, f) = [ cf2 "bm(()/ (s)/v/(()) dff((z) dæ for pairwise orthogonal gim. This approach yields inde- pendent features, so Kun is diagonal, but it is challenging lo find suitable sets of ort orthogonal functions. In both of these cases, the cross- coos-covv Cm. still depends upon the hyperparameters, and so there is usually little to no computational advantage (Burt et al. 2020a)..
Variational Fourier Features (VFF) (Hensman el. al. 2017) set (0m; f) to a reproducing kernel Hilbert space (RKHS) inner product between the harmonics on (a, b). < b E R and f in ID. Due to limiting the domain to compact subset, Fourier transforms become discrete that is, they become generalised Fourier series. Consequently, conditioning on only finite subset of the frequencies works, and this gives diagonal + low- rank structure in Kun for lower order one-dimensional Matérn kernels. However, the covariance functions are defined on IR x R rather than [a, b] x [a, b], so it is not straightforward to evaluate their spectra. Replacing the Euclidean inner product with an RKHS inner product permits to do this for lower order Matérn kernels, but this is not easily extended to other covariance functions, such as the spectral mixture kernel, or products of kernels. Moreover, in higher dimensions it is necessary to use tensor product of 1D Matérn kernels, which is limiting, [ensman et al. (2017) and Dutordoir et al. (2020) note that using regular grid of frequencies as in VFF significantly increases cost for D >>> 1 unecessarily, since features which are high frequency in every dimension are usually very unimportant. However, we demonstrate that it is possible to filter out these features (Section 4p.
This approach could be generalised by replacing the Fourier basis with some other basis. Then Om is calcu- lated using RKHS inner products with other basis functions, always yielding a hyperparameter independent Cm. with sparse matrices if the basis functions have compact support with little overlap. However, the need to calculate the RKHS norm of the basis functions for the elements of Kan limits these methods to kernels whose RKHS has a convenient explicit characterisation. In practice, this means using tensor products of ID Matérn kernels. The recent work of Cunningham et al. (2023) is a specific example of this which uses B-splines.
Dutordoir et al (2020) used spherical harmonic features for zonal kernels on the sphere, and this can be applied to RRD by mapping the data onto sphere. In this case the inducing features are well defined and independent, and this can be generalised to other compact homogeneous Riemannian manifolds. However, the harmonic expansion of l on the domain must be known; for isotropic kernels on DED restricted to the manifold, these can be computed from the spectral density (Solin & Särkkä 2020). Yet isotropy is too limiting an assumption; one can effectively incorporate different lengthscales in each dimension by learning the mapping onto the sphere, but Kuf also depends on this mapping, and so the cost returns to O(N
We seek a method which can be used with a broader class of covariance functions, but retains the key computational benefits.
5

---

Published in Transactions on Machine Learning Research (04/2024)
4
2
f(g)
0
-2
-2
-3
-2
2
-0.4
-0.2
0.0
0.2
0.4
(a) Prior
4
2
f(g))
-2
-2
-4
-3
2
-0.4
-0.2
0.0
0.2
0.4
2
(b) Conditioning on Fourier features
0.4
0.5
0.2
f(x)
f(E)
0.2
-05
*
*
-0.4
-1.0
-2
0
2
-0.4
-0.2
0.0
0.2
0.4
(c) Conditioning Integrated Fourier Features
Figure 1: Illusrationn of the Integrated Fourier Feature construction. We plot the mean function (dashed), between one and three standard deviations (shaded) and sample functions in both the data and frequency domains for a squared exponential kernel with unit lengthscale. The sample functions in the data and frequency domains correspond to one another.
(a) The prior's Fourier transform white Gaussian noise whose variance is given by the spectral density. (b) We cannot condition meaningfully on come finite collection of frequencies (red stars), this gives information about the other frequencies the conditional prior fl)) in the data domain is unchanged. (c) We show only the inducing values in the frequency domain, which are averages of the surroundig region. The conditional prior is now mcaningful, and the residual uncertainty is due to high frequency content not included in the features.

---

Published in Transactions on Machine Learning Research (04/2024)
Choosing It is well known that optimising inducing imputs is usually not worth the extra computational cost compared to good initialisation. We briefly review the initialisation methods for different features described above.
For inducing points, Burt et al. (0200 show that sampling from a k-DPP (determinental point process; with k M, and the kernel used in the DPP the same as the GP prior's) performs well, both theory they alternate between optimising
and in practice. Since that initialisation is dyperparaneeeerdepeneeeee the hyperparameters and sampling in variational expectation maximisation (EM) approach.
For VFF as described by Hensman et al. (017)) a rectangular grid of regularly spaced frequencies must be used, which they select (optimally in ID) to be centred around the origin. In higher dimensions, the regular grid leads to including suboptimal frequencies in the corners of the grid. A construction which leads to a more feature efficient set of frequencies is the following. Create a rectangular grid, and then discard the features which are not within a given cllipsoid, where the cllipsoid's axes should be chosen to the proportional the bandwidth of the spectral density in that dimension (which is inversely proportional to the lengthscale). This corresponds to discarding the corresponding rows in Kuf: and the corresponding row and column in Kun-
For spherical harmonics, the optimal choice is to use the frequencies which have the highest variance. For many covariance functions (for example, those constructed from monotonically decreasing stationary kernels on D using the method of Solin & Särkkä (020)) this corresponds to choosing the first M frequencies.
For B-spline features, a grid of regularly spaced basis functions which covers a rectangular domain containing the data is used, and the sparsity of the matrices is used to make the method efficient features which are not strongly correlated with the data also contribute less to the computational cost.
4 Integrated Fourier Features
We are not able, in general, to integrate out the Dirac delta in Equation [1 and retain the desirable com- putational properties without introducing further approximations. We propose to average Fourier features over disjoint intervals of width & (Figure [e]. and approximate the spectral density as constant over the integration width. We focus D to lighten the presentation here; to enforce that the intervals are disjoint we require Izm- 2ml for any + m'.
/m/te/
/ f(z)e 12TE= dE = =1 franter
te/2
(02m. f) def
_ e-
2m/2 (01., de
-//2
Am t2/2
c2(2', 2m) =
C1 (2
.) ze 12
(12)
6/2
/ +e/2
3m a +e/2
(2 Zmn') -2
E')
de de 2 5m-m (s((m)
(13)
2m/2 J 2m/-2/22
Now, 5 is the Kronecker delta. Note that if the intervals were not chosen to be disjoint then only the last line would change. The advantage of choosing disjoint intervals is to make Kan diagonal, which will simplify later analysis, as well as modestly reducing the computational cost. But recall that the inversion and log determinant of B still required, and this matrix remains dense.
4.1 Convergence
Now, if we calculate covariance matrices using the numerical approximation detailed above, and use this to evaluate the collapsed variational objective of Equation { we no longer have a proper variational objective in the sense of Equation (8). In order to distinguish the proper objective from our approximation, we introduce the notation 3 for the approximate objective.
In this section we show that the approximation & converges to C a reasonable rate as M + oo, for a broad class of covariance functions, showing that these features are efficient and produce good approximations for
7

---

Published in Transactions on Machine Learning Research (04/2024)
the purposes of hyperparameter optimisation. Since longer have the interpretation of reducing the KT, between posterior processes, provide additional results to give reassurance about the quality of the approximate posterior predictive distribution.
Firstly, we subtly transform the features to simplify the analysis. Note that applying an invertible linear transformation T to the features has no impact on inference or learning. That is, if we transform u' with mean and covariance Mu,, En K...- TK. Then from Equation 7 it follows that if we optimise after transforming, the optimal De = TT and the optimal jul Tpu', would be expected from optimising before transform- ing. Furthermore the collapsed objective of Equation ( and posterior predictive mcan and covariance of Equation are left unchanged.
Tw', then Kuf [l)
[[T'''((x']] TKws, and similarly
For the analysis, instead of normalising by the spectral density, we normalise by its square root. This has the advantage that we do not need any approximation for k, only for c. This simplifics later calculations.
(03m f) def
Fm t=/2
/
[(x)e
((o))--2nço 12n50 dx/ v a(E) = = f
Im Jzm-e/2
+e/2
(b) (V(())
Zm 2/2
Cm (a') = cg(x', =mm)= =
G (1 , ()()) de volzm)e 12nama = (x1)
2/2
Rm Am/2
Fm f
km.m² kimm Fm/) = =
(E(EEE)) de' de 6-
2m e2/2
2m/2
We proceed without the subscript 3 hereafter for brevity. The new approximate features are a straightforward invertible linear transformation of the previous features (in particular, 6m 03.m V (2m We analyse how well é approximates u, but in practice features ú2 in order to have è independent of the hyperparameters (Appendix B We defer the details of the
real-valued version of the equivalent approximate
proofs, particularly for higher dimensional inputs, to Appendix C
For convergence of the objectives, we use the following result of Burt et al.
2019
Ey, Eg(Dkr(o())
= -EylC - F
(14)
g² K KüjKau Kul)
t= t=t((Kii
(15)
Qu
where y is distributed according to Equation The first equality follows for any well-defined inducing features (Matthews et al. 6016)-that is, those where the inducing features can be priori described as a linear functional of the prior process. We now detail the additional assumptions.
A1 Let $ = s/o2²) be the normalised spectral density, which is assumed to exist, with e = b(z,z)/00 We assume that the normalised spectral density has a tail bound
B(E) dE S 3p
(16)
for any p >>0 and some $.q > 0.
A2 The second derivative of is bounded.
A3 The first derivative has relative bound
ds(E)
d s(E) < s({)
< 22s(EE)
(17)
dg
dg
for some L O, where the implication follows wherever s(E) 0. For example, this is satisfied if and its first derivative bounded everywhere, which is the case for widely used covariance functions.
8

---

Published in Transactions on Machine Learning Research (04/2024)
A4 The frequencies are chosen according to the regular grid 2m (m 1/2)e - M/2 for M even. This requirement can be relaxed in practice. The key requirement. is that Uml e/2, Zm e/2] - IR; could also be varied as a function of m, with the convergence rate dominated by the largest.
The higher dimensional generalisation of these are the standard assumptions. We use the following simple result, proved the supplement.
Lemma 4.1. Under assumptiona A3 and A4,
Cn (m( = 2m (x)
Theorem 4.2. For sampled according to Equation there ecists Mo. on > O such that
under assumplions A1 to A4, for any A, 6 > 0
C
IP
< 6
for all M 2 Mo, and with
00
M
S
Moreover, there erists M1,aj 01 > such that for all M > M1.
C F z A
P
for all M 2 M1, and with
2+3
M s
N)
24
A8
Proof. We sketch the 1D case here. Let f = tr(Kşt - Rū K Kut) = tr(Kí Qt)) First we show that
i/Na² E O0MM24//443))
1
Σ
k(In; Fnl) - Σ Cm
(18)
No²
No²
1122
1
3(2m)
(19)
(M+1))//2
5(E) d())
de
EEE
(20)
(M+1)e/2
2u
S(E) de v1i
(21)
Me/2
where E, = S (M+1)//2
J-(M+1)e/2 s(E) - 3(zm) +M The integral term in the last line is in O((M)) by assumption A1, and EΙ O from standard bounds on the error of the midpoint approximation (using A2).
We must have that - 0 as M - to make the midpoint approximation exact, yet we must have Me - to ensure the features cover all frequencies. By optimising the trade-off, we get the stated bound.
To complete the proof, for we could immediately apply Equation (14). But for f, we adapt the result of Equation (14) to show EylcC 3//N| s i/No² for sufficiently large M, using Lemma 4.1 (assumptiona A3, A4). By applying Markov's inequality, we complete the proof.
Remark 4.3. The case of subgaussian spectral densitics (such as for the squared exponential covariance function) is q + 0o, which yields M O(VN). Note that this is due the ((M=2) terms which arise due to approximating the spectral density constant. Intuitively, it appears that if there were numerical
9

---

Published in Transactions on Machine Learning Research (04/2024)
approximation, the cost would be dominated by the amount of spectrum in the tails, such that convergence for subgaussian tails would be possible with M E Oloo N) for sufficiently small with inducing points or eigenfunctions with the squared exponential covariance function Burt et al. (2019).
Though this demonstrates that the features suitable for learning, we may wish to the same features in making predictions. We no longer use the bound in the process KL, but we show that the posterior predictive marginals converge with comparable or better rate than the objective for any choice of variational distribution.
Theorem 4.4. For the optimised Pur En (to mazimise 3), let the posterior predictive at any test point 2+ using the exact features u have mean and variance μr, E, and with the approzimate features à have mean and variance ú, E. Then, under assumptions A3 and A4,
I A E O(ME2) Σ - E
In porticular, allowing 8 ~ M (913) as in the proof of Theorem
s+1
4.2
(see Appendir C we have
11 1p - Al M
9+3
(22)
3
Σ E
M
(23)
3+3)
Proof.
Use the definitions in Equation 7 and apply Lemma 4.1
and the triangle equality.
4.2 Higher dimensions
In higher dimensions, we modify the assumptions as follows.
A1 Assume that K's spectral measure admits density, and denote this by and assume that the density admits a tail bound
B(E) dei ... dED S 3p 9D
(24)
for any p and some 8,4 > 0.
A2
The spectral density's second derivative is bounded.
A3 The spectral density's first derivatives are bounded as
as(E) < 2Ls(E)
Ov8(E)
< LV(((E)
(25)
δEa
dEa
for some L (where the second expression follows wherever s(E) > 0). For example, this would be satisfied the spectral density and its first derivative are both bounded everywhere, which includes widely used covariance functions.
A4 Let the inducing frequencies be analogous regular grid in higher dimensions, with MI/D Z even. That is, for multi-index mI.D.
Iemran = (-(M+/D + 1)/2 ma)e for d {1: D}.
(26)
Then the main result (Theorem 4.2 and the predictive bounds in terms of q in Theorem 4.4 are unchanged. See Appendix C for details.
Assumption A4 requires full regular grid of features, which means the number of features increases expo- nentially in D. In common with previous work (Hensman et al.] 2017 Cunningham al. 2023) this
10

---

Published in Transactions on Machine Learning Research (04/2024)
Gaussian kernel
Matérn-3/2 kernel
1.4
1.2
Normalised E (eW
1.0
0.36
0.8
.50
170
0.6
0.4
0.2
0.4
0.6
0.8
1.0
1.2
1.4
18100
2.0
Normalised bandwidth (AME)
Figure 2: Gap between the log marginal likelihood and the training objective (C - 3) for different settings of M, e for data sampled from a GP with a Gaussian (left) or Matérn-3/2 (left) kernel. In each case the hyperparameters are set to their groundtruth values, where the lengthscale is Å. The inputs are samples from a uniform distribution centred on and with width Wa. The horizontal line is at 0.95.
be avoided the covariance function can be assumed be additive over dimensions. Otherwise, we can improve the computational cost by using a subset of the grid points where the spectral density is highest. Although this improves the cost in higher dimensions, the scaling is in general still exponential in D. This limits IFF's applicability to lower dimensional spatial spatiotemporal modelling applications.
4.3 Choosing the approximation parameters
Choosing .M These need to be selected whenever using SGPR. In IFF, the faster decays, the lower M we should need for a good approximation to the log marginal likelihood (Theorem 4.2) and Figure / shows that M need not be very large (in each dimension) to get good performance; we need Me at around the approximate bandwidth of the covariance function, which increases as the lengthscale À reduces. If we know priori how small the lengthscale might become - for example by examining the Fourier transform of the data, from prior knowledge, or by training a model on small random subsets of the data then we can use this to select M. Note that if the lengthscale is shorter, we would expect M to need to be larger for SGPR also, as we would need inducing points to be placed closer logether. In practice, as with all SGPR methods, M controls the trade-off between computational resources and performance, and can be set as large as needed to get satisfactory performance for the applicaiton, or as large the available resources allow.
Choosing 2 Given a fixed E, M, the optimal choice to choose frequencies spaced by E in the regions of highest spectral density, comparable to spherical harmonies. For monotonically decreasing spectral densities maximised at the origin, such the Gaussian Matérn kernels, this corresponds exactly to our refined construction for VFF in Section g which involved choosing grid points which were contained within an ellipsoid whose axes were inversely proportional to the lengthscale in each dimension. To avoid dependence on the hyperparameters (which would undo the benefits being able to precompute A'), we opt in practice for spherical threshold.
Choosing s When adding more features, we can either cover a higher proportion the prior spectral density or reduce e. In the full proof of Thcorem 4.2 (Appendix C. as the spectral density gets heavier tailed, approaches O(M-1). Then ME is almost constant, which suggests this part tends to dominate. That is, once Me is sufficiently large, should add more features by reducing rather than increasing
11

---

Published in Transactions on Machine Learning Research (04/2024)
M. However, we explore the trade-off between increasing Me and decreasing numerically (Figure by plotting the gap between the IFF bound and the log marginal likelihood, varying the bandwidth covered and the size of e relative to the inverse of the data width (W₂ R max, 2n min, Tn). The model and data generating process use a kernel with lengthscale X. We see that in practice, as long as E is below the inverse of the data width, the gap is not very sensitive to its value. Thus for the experiments, we conservatively set. (000//mmxp... (maxn,d Fnd and))
ed = 0.95/
This result appears to be in tension with the theory. However, we note that in the proof we are effectively interested in producing good approximation of the function across the whole domain (we make no as- sumptions about the input locations). But in Figure g we explicitly take into account where the data is. Intuitively, our construction involves approximating the function with regularly sampled frequencies. Then it should be possible to construct a good approximation to a function on an interval of width Wa as long as the sample spacing is no more than 1/W,, as Fourier dual to the classical Nyquist-Shannon sampling theorem (soe, for example, Chapter 5, Vetterli et al. 2012
Other covariance functions We have so far assumed that the spectral density is available in closed form. However, we only need regularly spaced point evaluations of the spectral density, for which it suffices to evaluate the discrete Fourier transform of regularly spaced evaluations of the covariance function. This adds, at worst, O(M?) computation to each step.
4.4 Limitations
IFF can be used for faster learning for large datasets in low dimensions, which matches our target applications. Typically, it will perform poorly for D z 4, and both in this case and for low N, we expect SGPR to outperform all alternatives, including IFF, and analysis and evaluation limited to the conjugate setting.
IFF is limited to stationary priors; while these are the most commonly used, they are not appropriate for many spatial regression tasks of interest, and a fast method for meaningful non-stationary priors would be a beneficial extension. Amongst those stationary priors, we require that the spectral density is sufficiently regular; this is satisfied for many commonly used covariance functions, but not for periodic covariance functions, where the spectral measure is discrete. While IFF can be used with popular covariance functions for modelling quasi-periodis functions (such functions, or the spectral mixture kernel) if the data has strong periodic components, the maximium marginal likelihood parameters will cause the spectral density to collapse towards a discrete measure (for example, learning very large longthsealcs).
product of squared exponential and periodic covariance
Finally, we have left. some small gaps between theory and practice, both in how to select the tunable parameter s, and in characterising the quality of the posterior predictive distribution.
5 Experiments
We seek to show that IFF gives a significant speedup for large datasets in low dimensions, with particular focus on spatial modelling. Amongst other fast sparse methods, compare against VFF and B-Spline features. For spherical harmonics, learning independent lengthscales for each dimension is incompatible with procomputation. In any case, we found that we were unable to successfully learn reasonable hyperparameters with that method in our selting, except if the nunber of feature was very small. For a conventional (no pre-compute) sparse baseline, we use inducing points sampled according to the scheme of Burt et all (2020a). For our synthetic experiments, we also used inducing points initialised using k-means and kept fixed. For the real-world spatial datasets, we also tested SKI, due to its reputation for fast performance, and its fairly robust implementation.
12

---

Published in Transactions on Machine Learning Research (04/2024)
0.030
0.025
Gap to L
0.020
0.015
0.010
0.005
0.000
128
64
32
16
128
64
16
Execution time (normalised)
Execution time (normalised)
0.04
Gap to groundruth C
0.03
0.02
0.01
0.00
Execution time (normalised)
Execution time (normalised)
IFF (ours)
IFF maive (ours)
0.04
VFF redaced
VFF (Hensman al, 2017)
Inducing points k-DPP
Gap to C
0.03
(Bart et al, 2019)
Indocing points k-means
0.02
0.01
0.00
0
200
400
600
800
1000 0
1000
2000
3000
4000
M
M
Figure 3: Comparing standard sparse Gaussian process regression (black) to IFF (red) for data generated from prior with Gaussian covariance function in 1D (left) and 2D (right). Lower and the to the left is better. The groundtruth L is evaluated at the groundtruth hyperparameters, whereas in the other rows, L is evaluated at the learnt hyperparameters. The gaps are normalised by N, and execution time is normalised by the longest. The bottom row shows feature efficiency, whereas the upper rows show computational efficiency.
13

---

Published in Transactions on Machine Learning Research (04/2024)
0.04
0.03
Gap to C
0.02
0.01
0.00
64
Execution time (normalised)
Exccution time (normalised)
0.04
Gap to groundtruth L
0.03
0.02
0.01
0.00
256
128
64
256
128
64
32
16
Execution time (normalised)
Execution time (normalised)
0.04
IFF (ouns)
IFP naive (ours)
0.03
VFF (Hensman al, 2017)
Gap to L
[Burt dl 2019)
Inducing points 1-DPP
Inducing points k-means
0.02
0.01
0.00
200
400
600
800
1000
1000
2000
3000
4000
M
M
but with the data sampled from a Matérn-5/2 GP. The picture broadly comparable, but VFF now more closely matches the prior, so the drop feature efficiency is far less higher dimensions.
Figure 4: As Figure 3
14

---

Published in Transactions on Machine Learning Research (04/2024)
5.1 Synthetic datasets
First we consider synthetic setting where the assumptions of Theorem 4.2 hold. We sample from a GP with a Gaussian covariance function, and compare the speed of variational methods in 1 and 2 dimensions. We use a small (N = 10000) dataset in order that we can easily evaluate the log marginal likelihood at the learnt hyperparameters. Where possible, we use the same (squared exponential) model for learning; for VFF, we use a Matérn-5/2 kernel in 1D, and a tensor product of Matérn-5/2 covariance functions, since this is the best approximation to Gaussian kernel which is supported. Further details are in Appendix D
Additionally, in the 2D setting, we use both the naive set of features (a regular, rectangular grid), and the refined set of features described in Section 3
IFF generally has slightly lower gap to the marginal likelihood at the learnt optimum for any M than other fast variational methods (Figure 3 bottom row), but because the O(NM²) work is done only once, and the other fast. sparse methods are much faster to run than inducing points (Figure 3 top two rows). Note the logarithmic time scale on the plots: for a specified threshold on the gap to the marginal likelihood, IFF is often around 30 times faster than using inducing points.
The experiments demonstrate the issues with the limited choice of prior with methods such as VFF. In ID, the Matérn-5/2 kernel is a good approximation to Gaussian, so the performance is similar to IFF. But for 2D, the product model is a much worse approximation of the groundtruth model. We expect to see a similar pattern for B-spline features, but we were unable to run the method in our syuthetic setting due to unresolved issues in the implementation of Cunningham et al.] (2023) which did not arise in the real-world experiments. When we reproduce the same experiment, but with data sampled from a Matérn-5/2 GP, the results are similar, but with much smaller gap between VFF and the other methods in the 2D case (Figure 4/.
We note that in the ID setting, when the data is sampled from a prior with Gaussian covariance function, the k-DPP method gets stuck at a local optimum, which is a known risk with that method. We did not find this to be an issue in any of the other experiments.
The refined feature set for higher dimensional IFF and VFF (solid lines) is indeed slightly more feature efficient than the naive approach (dashed lines; Figure 3 bottom right panel). But in practice, we find that the computational overhead to selecting better features outweighs the savings from using fewer features for VFF, thought not for IFF (Figure 3 upper and middle right pancls).
5.2 Real World Datasets
We now compare training objective and test performance on three real-world spatial modelling datasets of increasing size, and of practical interest. We plot. the root mean squared error (RMSE) and negative log predictive density (NLPD) on the test set along with the training objective and run time in Figures 5 and 6 using five uniformly random 80/20 train/test splits. For inducing points, we always use the method of Burt et al. (020b). The time plotted is normalised per split against inducing points. Further training and dalaset details are in Appendix D.
We add SKI to the comparison here, since it is the most widely used alternative to variational methods in this setting. For the variational methods, both the time and performance are implicitly controlled by the number of features M; as M is increased the performance improves and the time increasos, that is we move along the curve to the right. This is very useful property, since we can select M according to our available computational budget and be fairly confident of maximising performance. With SKI, the equivalent parameter is the grid size, and similarly increasing the grid size generally improves performance. However, when the grid size is low, optimisation can take longer, so for example for the temperature dataset, we move to the left. as the grid size is increased (Figure G top row). For the other datasets, we only plot the SKI with the grid size automatically selected by the reference implementation.
SKI generally has very good predictive means, leading to low RMSE, and is very fast, but the predictive variances are poor, generally being too low and producing some negative values. Ignoring the negative
15

---

Published in Transactions on Machine Learning Research (04/2024)
Test RMSE
Text NLPD
Training loss
0.45
1.0
0.2
0.40
0.0
Temperature
0.35
0.6
0.2
0.30
-0.4
0.2
0.6
0.0
1151
M72 250
1584
012
270
1054
M12
350
0.6
4.6
24
0.5
Precipitation
4.4
20
0.4
4.3
0.3
16
4.1
0.2
2008 1004 502
2002
2004
312
350
126
3048
1004
512
300
0.78
1.22
0.51
0.76
1.20
House price
0.50
0.74
1.18
0.49
0.72
1.16
0.70
0.48
1.14
0.68
0.66
1.12
2048 1024
032
2048
1024
512
250
125
16
2548 1024
312
Time [pormalisex]]
Induring
IFF (ouns)
VFF (Hemenan al. 2017)
Bxpline (Cumingham et 2223
Figure 5: Performance curves for real world datasets of increasing size (the top row is the smallest). Lower and the left better.
16

---

Published in Transactions on Machine Learning Research (04/2024)
Test RMSE
Test NLPD
Training loss
0.45
0.2
1.50
0.40
1.25
0.0
Temperature
1.00
0.35
-0.2
0.75
0.30
0.4
0.50
0.25
0.25
0.66
0.20
0.00
1003
612
256
0.6
0.55
Procipitation
4.6
20
0.4
4.4
18
0.3
16
14
0.2
1026072
20
31
1021077
200
22
1924 002 270
0.32
1.225
1.1
1.200
0.50
1.0
House price
1.175
0.48
1.150
1.125
0.46
1.100
0.7
0.44
226
312
(normalisec)
Incuiing points IFF (ours)
VFF (Himetmin et al, 2017) appinsss (Cunninglan al, 2023)
SKI (Wilson & Nickedh 2015)
Figure 6: As Figure 5 but with SKI included. SKT exhibits very favourable and fast predictive mean performance, but its predictive variances are poor, leading to very large NLPD.
17

---

Published in Transactions on Machine Learning Research (04/2024)
values, the NLPD is much worse for SKI than for the variational methods, which is highly undesirable in probabilistic method.
We exclude SKI in Figure 5 in order to zoom in on the curves for the variational methods. We are interested in the regime where M < N; as we move to the right and M is similar to N, inducing points will become competitive with the faster methods, since the (MS) cost dominates. Comparing TFF VFF, we see that always performs at least as well, and produces substantially better performance for a given time on the temperature and precipitation datasets, due to more flexible choice of covariance function in particular, note that the training objective of IFF is substantially lower than VFF on these datasets, but comparable to that of inducing points, which uses the same covariance function.
The B-spline features are also limited in choice of covariance function, but the sparse structure of the covariance matrix leads to a much better performance. Nonetheless, despite involving a dense matrix inverse, IFF has a significant advantage on the precipitation dataset (around 25-30% faster on average), and is comparable on the houseprice dataset.
Compared to the idealised, synthetie, setting, inducing points become very competetive in the higher resource setting towards the right hand side of cach plot, particularly for the smallest dataset (temperature). But for more performance thresholds, we find that fast variational methods offer a substantial improvement, with typically more than a factor of two speedup.
6 Conclusions
Integrated Fourier features offer a promising method for fast Gaussian process regression for large datasets. There are significant cost savings since the O(N) part of the computation can be done outside of the loop, yet they support a broad class of stationary priors. Crucially, they are also much easier to analyse than previous work, allowing for convergence guarantees and clear insight into how to choose parameters.
They are immediately applicable to challenging spatial regression tasks, but significant limitation is the need to increase M exponentially in D. Further methods to exploit structure in data for spatiotemporal modelling tasks (D (D=3 3 or 4) is an important line of further work. More broadly, an interesting direction is to consider alternatives to the Fourier basis which can achieve similar results for non-stationary covariance functions, which are crucial for achieving stale of the art performance in many applications. Finally, a worthwhile direction for increased practical use would be to develop suitable features to quasi-periodic priors such those arising from the spectral mixture kernel, or product of periodic and Gaussian kernels
References
Vladimir I Bogachev. Gaussian Measures. American Mathematical Society, 1998. ISBN 978-0-8218-105445
David Burt, Carl Edward Rasmussen, and Mark van der Wilk. Rates of convergence for sparse variational Gaussian process regression. In 36th International Conference on Machine Learning (ICML), 2019.
David R Burt, Carl Edward Rasmussen, and Mark van der Wilk. Variational orthogonal features, 2020a.
URL bttpp::/aaxxv..rg/abb///////
13170
David R Burt, Carl Edward Rasmussen, and Mark van der Wilk. Convergence of sparse variational inference in Gaussian processes regression. Journal of Machine Learning Research (JMLR), 21(131):1 63, 2020b.
URL http://jmlr: http://rw..rr//ppppp///////11000 html
Michael K Cohen, Samuel Daulton, and Michael A Osborne. Log-linear-time gaussian processes us- ing binary tree kernels. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and
A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 8118 8129. Cur- ran Associates, Inc., 2022. URL https::/prooeeddings neurips. http::/s....ppppe....ssa/////p/p /iles//aper/002///iiel BBBcdbbacccsseccccceeeeeee-----e
Harry Jake Cunningham, Daniel Augusto de Souza, So Takao, Mark van der Wilk, and Mare Peter Deisen- roth. Actually sparse variational gaussian processes. In Francisco Ruiz, Jennifer Dy, and Jan- Willem de
18

---

Published in Transactions on Machine Learning Research (04/2024)
Meent (eds.), Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, volume 206 of Proceedings of Machine Learning Research, pp- 10395-10408. PMLR, 25-27 Apr 2023. URL https://proceedings mlr. vrres/vv006/ 2tt220///uunnngan22a2a... html
John P Cunningham, Krishna Shenoy, and Maneesh Sahani. Fast Gaussian process methods for point process intensity ostimation. In 25th International Conference on Machine Learning (ICML), pp. 192-199, 2008.
Vincent Dutordoir, Nicolas Durrande, and James Hensman. Sparse Gaussian processes with spherical har- monic features. In 37th International Conference on Machine Learning (TCML., 2020.
Yarin Gal and Richard Turner. Improving the Gaussian process sparse spectrum approximation by repre- senting uncertainty in frequency inputs. In 3g International Conference on Machine Learning (ICML), 2015.
Jacob Gardner, Geoff Pleiss, Ruihan Wu, Kilian Weinberger, and Andrew Gordon Wilson. Product kernel interpolation for scalable Gaussian processes. In 21st International Conference Artificial Intelligence and Statistics (AISTATS), 2018a.
Jacob R. Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, and Andrew Gordon Wilson GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration, 2018b. URL https://arre org/ a////80//809 11165
Roman Garnett. Bayesian Optimization. Cambridge University Press, 2023.
Philipp Hennig, Michael A Osborne, and Hans P Kersting. Probabilistic Numerics: Computation as Machine Learning. Cambridge University Press, 2022. doi: 10. 10.1017/9781316611111
James Hensman, Alexander G de G Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification, In 18 International Conference on Artificial Intelligence and Statistics (AISTATS),
2015.
James Hensman, Nicolas Durrande, and Arno Solin. Variational Fourier features for Gaussian processes. Journal of Machine Learning Research, 2017.
Martin Jø rgensen and Michael A Osborne. Bezier gaussian processes for tall and wide data. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Ad- vances in Neural Information Processing Systems, volume 35, pp. 14854-24360 Curran As- sociatos, Inc., 2022. URL https://prroceedlngg newriseee.ppp... wewrips.cc/papeeee cc/paper filess/paperr2222////
Miguel Lázaro-Gredilla and Aníbal Viddal. Inter-domain Gaussian processes for sparse inference using indueing features. In 26th Conference on Neural Information Processing Systems (NeurIPS), 2009.
Miguel Lázaro-Gredilla, Joaquin Quinonero-Candela, Carl Edward Rasmussen, and Aníbal R Figueiras- Vidal. Sparse spectrum gaussian process regression. The Journal of Machine Learning Research (JMLR), 11:1865-1881, 2010.
Mikhail A Lifshits.
Lectures on Gaussian Processes. Springer, 2012. ISBN 978-3-642-24988-9)
Alexander G de Matthews, James Hensman, Richard Turner, and Zoubin Ghahramani. On sparse variational methods and the kullback- leibler divergence between stochastic processes. In Arthur Gretton and Christian C. Robert (eds.), Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, volume 51 of Proceedings of Machine Learning Research, pp. 231 129-299,422 Cadiz, Spain, 09-11 -1 May 2016. PMLR. URL https::/ mlr. www.mmmr..errrre
html
Alexander G. de G. Matthews, Mark van der Wilk, Tom Nickson, Keisuke. Fujii, Alexis Boukouvalas, Pablo Leôn- Villagrá, Zoubin Ghahramani, and James Hensman. GPflow: A Gaussian process library using Tensor Journal of Machine Learning Research (JMLR), 2017.
19

---

Published in Transactions on Machine Learning Research (04/2024)
Geoff Pleiss, Jacob Gardner, Kilian Weinberger, and Andrew Gordon Wilson. Constant-tinns predictive distributions for Gaussian processes. In 35th International Conference on Machine Learning (TCML), 2018.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In 24th Neural Information Processing Systems (NeurIPS), 2007.
Carl Edward Rasmussen and Christopher K 1 Williams. Gaussian Processes for Machine Learning. MIT Press, 2006. ISBN 978-0-262-18258-9,
Yunus Saatçi. Scalable Inference for Structured Gaussian Process Models. PhD thesis, University of Cam bridge, 2011.
Simo Särkkä, Arno Solin, and Jouni Hartikainen. Spatiotemporal learning via infniite-dimensional Bayesian filtering and smoothing: A look at Gaussian process regression through Kalman filtering. IEEE Signal Processing Magazine, 2013.
Arno Solin and Simo Särkkä. Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Computing, 2020.
Michalis Titsias. Variational learning of inducing variables sparse Gaussian processes. In 12th International Conference on Artificial Intelligence and Stalistics (AISTATS), 2009.
Gia-Lac Tran, Dimitrios Milios, Pictro Michiardi, and Maurizio Filippone. Sparse within sparse Gaussian processes using neighbor information. In Marina Meila and Tong Zhang (eds.), Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 10369 10378. PMLR, 18 24 Jul 2021. URL htttps://proceeddingg mlr. wwt.e/ereese//tteaeeeeeee html
Martin Vetterli, Jelena Kovačević, and Vivek K Goyal. Foundations of signal processing. 2012.
Andrew Gordon Wilson and Hannes Nickisch. Kernel interpolation for scalable structured Gaussian processes (KISS-GP). In 32nd International Conference Machine Learning (ICML), 2015.
James T Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Deisenroth. Efficiently sampling functions from Gaussian process posteriors. In 37th International Conference on Machine Learning (ICML), 2020.
Luhuan Wu, Geoff Pleiss, and John P Cunningham. Variational nearest neighbor Gaussian process. In Kamalika Chaudhuri, Stefanie Jegelka, Te Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (eds.), Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 24114-24130. PMLR, 17-23 Jul 2022. URL https::/ppoocecingss... mlr. preas/v162/wu22. htmm
20

---

Published in Transactions on Machine Learning Research (04/2024)
A
Computation
Recall the collapsed objective.
1
F(uu; En) (ylo, K Kus + 0²1) tr((K Kuj kut)
(27)
1 =- log Kip Kud Kuf 0² 2"
(Kaf KKu KunKu Kuj + o² 2/ tr((g Kus K Kay)
(28)
With y= Kufy, En U2 v², we apply the Woodbury identity to the inverse in the quadratic form.
K Kuf +
- o= KK + t Kut Kuf 0²1)=1 +
y K Kus + o yy= o -22 o (Kuu + a Kuf -
For the log determinant, we can use the matrix determinant lemma.
1 Kij Kuil= = (Ku + to "K
-2 KuiKKullKKK 102 Inl
Finally, we write down the trace directly using the fact that Kua is diagonal.
tr(Kp Kij K uu Kus)
Σ k(Inn In)
Σ ICm (xm) P//mram
We combine the above to get an easy to evaluate expression for F, and replacing the inter-domain cross covariance matrices with their numerical approximations,
evaluate the IFF objective 3.
Notably, x² E Rz0 depends only on y, can be precomputed and stored with only O(N) cost, and ÿ E R- depends only on æ, and z, so can also be precomputed and stored with O(NM) cost. Similarly, the matrix KufKū, can precomputed with O(NM²) cost. For large N, we split the data into chunks of 10 000 to save memory in this precompute stage. During optimisation, each calculation of this objective is reduced to the ((M)) cost associated with the inverse and log determinant caleulations, which perform cheaply after first performing the Cholesky decomposition.
The prediction equation is
((f(x=)(x,, I.)) = / p(f(o.)l,, 2, (N(((x))) Kid pus K.. Kū Km Kas K Kr1 1-1E E.g Kyd K
(29)
where & stands for I. in the subscripts. This is just the sparse, inter-domain version of Equation Pas En are given in Equation
and
B Real-valued features
As noted in Section applying an invertible linear transformation to the features does not change inference
or learning.
To simplify the presentation, and generalise the result, we change notation slightly from the main text. Let
the complex valued featured be " which is the feature corresponding to the frequency zmp...,.ppp
and let be antisymmetric in every axis (that is, for each d.
-m1 mm etc) and require ma +
For example, if we use a regular grid, with M1/D an integer,
E/2 + mie M1/D/2
Zem DDD
E/2+ + MDD - M1/DD/2
21

---

Published in Transactions on Machine Learning Research (04/2024)
which indeed satisfies this property. Now,
12m
cos(22zx' T {) + isin(2 {)
D
Σ
Σ
II sin 2redad cos 2m
1=0
dgg
so let mDD, for positive my only, with S $ {1: [}. be defined as
/ f(z) II
Jn-6/2
sin a
II
cos a drde: ... dED
SD-0/2
dES
d#S
that
( l. II sgu(ma).
E(-i)(s1
dES
Lemma B.1. The real representation is equivalent to the complex representation used the main text.
Proof. It suffices to show that this transformation be expressed matrix with lincarly independent rows. The row corresponding lo each v' mm....mp has only non-zero entries for Upa I.../mmls for any S. Hence, if the absolute values of mı,...; mp differ, then the rows are linearly independent. Now, suppose the absolute values are fixed and consider an arbitrary collection of indices S C {1: D} to have negative sign. This corresponds to a particular u', hence a particular row the matrix. Then the non-zero entrics in the row corresponding to S' each correspond a choice of S, and the sign flipped (relative to the case where each ma positive) if JS | odd.
The rows only linearly dependent if they have all their signs flipped all their signs not flipped. That is, if there exists S' + S" such that |S S"I and [S n niS' have the same parity for all S S {1 : D}. But since they are distinct, there must be at least one E {1 D} which is in S' but not S" (or vice versa) and so for S = {d], the parity differs.
C Convergence
We follow the steps of Section 4.1 generalising to higher dimensions and filling in the details.
Recall that inference and learning with with the modified features ug and their approximation úg is equivalent to inference and lcarning with the features use in practice. In this section, for brevity, use à ûg. and u ug-
Following from Burt et al. 2019 the gap between the log marginal likelihood and the training objective is bounded as
(xx EE,IC -F
(30)
6² Kúp K Kut)
t t= tr(KK
(31)
Qu
when are valid inducing features, and the data are generates according to the Equation 1 We defer the effect of approximating with â until later in the proof. First, we set out the technical assumptions.
Let = v where k(xxx) (that is, define v as the signal to noise ratio).
A1 Assume that i's spectral measure admits density, and denote this by and assume that the density admits tail bound
B(E) dei den S 8p RD
(32)
for any p and some 8,4 0.
22

---

Published in Transactions on Machine Learning Research (04/2024)
A2 The spectral density's second derivative is bounded.
A3 The spectral density's 's first derivatives are bounded as
Ōs(E) dg
(v)(5) dg
S (L((E)
S L.vs(E))
(33)
for some L > 0 (where the second expression follows wherever s(g) > 0). For example, this would be satisfied if the spectral density and its first derivative are both bounded everywhere, which includes widely used covariance functions.
A4 Finally, let the inducing frequencies Zm (-(M + 11/2 + m)e in one dimensions, and an analogous regular grid in higher dimensions, with M E Σ even. That is, for multi-index m1.D.
lama (-(M-/D + 1)/2 ma)s for d {1: D3.
(34)
Notationally, we usually use a single index m and, for brevity, let
II [am]a 8/2, [zmld + E/2),
0= -= -=UDm
(35)
We might wish to deviate from A4 in practice for example, to prioritise higher importance frequencies when using a finite number, or to vary in each dimension or as function of location. We note that the proofs which follow can be generalised to these cases, with the rate of convergence controlled by the largest width used
In the rest of this section, we refer to these as the standard assumptions, and we now reiterate the definitions of and à.
-- S(// [(s)e 12n50 (enee dx/ s(E) de = = S (Vv.s)))ss())
(ds,, 1) =
(xm) cg(x', 2m)) = =
( ~ v Vs(ēmje Je 12xmmm0' = cm
n.
km.mm = ka(āms Zml) = = [ / (((E) d== = 6m-m'
D...
Then à is defined implicitly through the definition of Ems kmmm above.
Lemma C.1 (Lemma Under assumptions A3 and A4,
(m)))
(36)
Proof. We first consider the case. Let (m(E) = Vo(E) (8(2m). Then by Taylor's theorem, |Em (E)
S
LiS zm (onlvs(() for some e' E Dm- In particular, let g' = arg maxge Vs(E). But we also have VE" <
V LV 5(E')E'' 2ml V + s/s()))// Then for < 2/L L it follows that
V /s(E") S V ((m))
1 LE Vo(2m) + O(s))
(37)
and so
(E) T.gg 2ml ( (1 0()).
(38)
23

---

Published in Transactions on Machine Learning Research (04/2024)
Then, considering at first the upper bound
/ Vs(E)e 1255 de
(m)x)= === Dm
(39)
=E- / (Vo(zm) + Em(E))e
E 2m dE
(40)
12 -1Rssma Vo(mm)e
S
1-i2m(8-2m1
de + t / Em(E)e -2205= -12m5 de
(41)
D...
f [Em(E) de
S (2T + +
(42)
LE
S (x( (inn((2TEx) +
(43)
2
Here sinc(a) = sin(a)/a. The sine term of constant order. The lower bound is found by subtracting the magnitude of the error term instead of adding. The result then follows.
In higher dimensions, we follow the samne argument, and the new upper bound is
LE
Cm (x) S (x)
ssnnc (2TEE) +
(44)
from which the result. follows.
Remark C.2. The error bound in Equation 44) generally tightens as falls, but loosen as increases, with
the approximation value vanishing when 1/2.
Theorem C.3 (Theorem 4.2 of the main text). Under the assumptions A1-A4, for any A 0, there
epists Mo, 0o for all M Mo,
C 8
P
8
N
N
and with
M <
a
N)
Ad
Proof. Let = tr(K5t - Rickod K-d Rut) = tr(Kn - Q11) follows.
Then we can show that i/No² E O(M 20/(9+31) as
Σ
k(Zng -5 Σ 3m (2n (2m)
(45)
N
No²
n02
Σ 3(2m)
(46)
/ ((E) de
(47)
2v
3(E) de VEi
(48)
RRL
where E = (c) - 8(zm) + 3(M33D). The integral term in the last line is in ((A by assumption A1, and E, E O(Me from standard bounds on the error of the midpoint approximation (since the second derivative of the integrand is bounded by assumption A2).
We must have that & - as M + so to make the midpoint approximation asymptotically exact, yet we must have Me + to ensure the features cover all frequencies. We optimise the trade-off between these two.
24

---

Published in Transactions on Machine Learning Research (04/2024)
In particular, let EgM- p/P for some Eo and p (1/3, 1). Then,
i
No O(M 4(1 1(1-p) + M -0 (8)-
(49)
The overall rate asymptotically dominated by the worse these two rates, so we optimise p as which is the which sets both rates equal at
= arg max min - p), - 1} q +
(50)
+ 3
2q q +
Altogether have
No O(M 2+3)
(51)
For bounding E[(- - in terms of i, we could immediately apply the result Equation [ B)/N] we cannot, since à are not exact features. But following the proof of Lemma of Burt ct al. 2019 we have
30 But for
c 3 S
i
By
No² +0(M33D)
(52)
No²
provided Q11 2 and (log Q11 + 0 - log Kjj + o2II//No2 < 3(M For the first condition, apply Bochner's theorem. Consider the covariance function
(0,, x') Σ
(53)
m
which is used to form the elements of 1n. Its Fourier transform is
[F(((E) = Σ (2m))8(E 2m
(54)
which is indeed a positive measure, so q is a positive definite covariance function, so Q11 2 0.
For the log determininant term, we have from Lemma 4.1 (which holds due A3, A4) that the relative error of Qfí from Qi is symmetric and in O(M)*D), hence the error in the log determinant is in O(NMe the scaled identity shift does not change the order of this relative error. Then,
(log 11 + 0²2 log K + o211)/NO2 < (log 14 + 0²1 - log K + o2Il)/No² + 3D)
(55)
E (OME))
(56)
where the last step follows from (log 1Q + 0²T - log Kft + o²Tl) & 0 (see proof of Lemma 2 of Burt et al. (2019h). Thus, we have
L
3 E O(M
Ey
(57)
N
Then the first part of the results follow by straightforward application of Markov's incquality. For the second part, we replace the big O notation with an explicit constant ag > 0.
L $
85 E[[C /N
(58)
A/N
N
M
(59)
M S
N
(60)
6A
7
25

---

Published in Transactions on Machine Learning Research (04/2024)
This result demonstrates that the objective use for hyperparameter optimisation, even with the numerical approximation, converges at reasonable rate to the log marginal likelihood, at least. if the spectral density has sufficiently light tails. This means it is a good surrogate for learning when we can set M large enough.
With the proper inducing features u, we can be reassured that the posterior predictives will not be too bad, since the whole process KT, from the approximating to exact posterior bounded as - Matthews et al. (2016). With û, we require some additional reassurance, which the is the subject of the next result.
Theorem C.4 (Theorem 4.4 of the main text). For the optimised Pa, Eu (to mazimise 8); let the poste- rior predictive at any test point I, using the exact features u have mean and variance E, and with the approzimale Jealures â have mean and variance it, E. Then, under assumplions A3 and A4,
1p - A O(M)2D) E - O(M2
(61)
(62)
412
and in particular, allowing M Dig+3) as in the proof of Theorem C.3
we have
| - M E M- 9+3
(63)
3
Σ - Σ E 0 M +
(64)
This result shows that the predictive marginals using the approximation converge a reasonable rate to those without the approximation, for any fixed Mhay Σn as long as the spectral density is not too heavy tailed (g 3). We note that we do not comment the rate at which the optimal variational means and covariances converge to each other, nor the consequent rate of convergence for the posterior predictives according to each objective's optimal variational distribution.
Proof. The predictive distributions are
(f((.)) = -N(((x.))
K K May k(x,, xxn) - Ky Kml Kut + K Kyr EuKu Kut
(65)
5.
a(((.)) = -N((xx)) Ris Katthu k(oa) ) - Ru Kmd Rus + Rū Kud 1r-1 En Kad Ru.
(66)
Recall that K E and that [[K.. Kurlml (cm(xx)) cm (2.) O(=P) by Lemma 4.1 for the means follows
The result
l A = I(K,, - Rin)kad E 2(M
(67)
which follows since K Pa. is an M dimensional vector with each element in O(eP). For the covariance, we use the triangle inequality.
Σ S KK 1a Kun Kink Karl E
(68)
+ K K K Kur - R2 K K. K- Rael
(69)
Now, each of the terms on the right hand side is a one dimensional marginal variance, so we rewrite them as follows.
K Ki! Kus
R K,11 Kasl =
"
-1mm -1
(70)
E 2
(71)
Kn Kur R Kid En K- -
2D Σ
[ulmm (Cm (om(xx)) (za) cm (x.)) (Ix))
(72)
m.m!
(33)
The first term is En (x,))2 Iêm (x=) (x+)² ) E 3(M)32) following from the proof of Theorem other two terms of the same order 2D) using Lemma 4.1
C.3 The
26

---

Published in Transactions on Machine Learning Research (04/2024)
D Experimental Details
We include the code for the experiments and figures which can be referred to for full details.
For Figure 2 we used N = 1 000 data points, whose input locations were sampled from a zero mean Gaussian distribution with standard deviation Wz/6. We randomly sampled a signal standard deviation of ~ 0.5422 and a signal to noise ratio (SNR) & Rs 1.255, and used a unit lengthscale (A 1). To verify that the pattern is broadly consistent for different parameter choices, we also reproduced the figure with lengthscale A 0.5, and, with both lengthscales, we used both half and four times the SNR. In each case, we set W. 6Λ / Flnally, we consider the case with uniformly sampled data. The results are plotted in Figures 7 to 12
For the synt hetic experiment, we generated N = 10000 data points in 1 and 2 dimensions by samping from a GP with a Gaussian or Matérn-5/2 covariance function, with unit (or identity) lengthscale, unit variance, and set the SNR to 0.774 (arbitrarily chosen, poor signal to noise ratio for a challenging datasct). In 1D we sample the training inputs uniformly on width G/N/2 centred interval. In 2D we do the same in each dimension, but with width 5. We then fit each model plotted, training using LBFGS and using the same initialisation in each case, other than necessary restrictions on the choice of covariance functions described in the main text. The initial values were lengthscales of 0.2, and unit signal and noise variances. We did multiple random trials, and plot the 2 standard deviation error bars; this is mainly for uncertainty in timing, but for inducing points we also have uncertainty due to randomness in the inducing point (re)initialisation method. Inducing points with inducing inputs optimisod takes far longer than the other methods, so was excluded.
For the real-world experiments, we used a similar sctup the synthetic experiments in terms of initialisations. Guided by the synthetic results, we use the full rectangular grid of frequencies for VFF, but use spherical mask for IFF. We set as described in the main text. Comparably, for B-splines and VFF, we set the interval (a, bj as 0.1 wider than the data (that is, Og = minn,d Fnnd - 0.1, by = maX,,d Fn,d + 0.1); we use fourth order B-splines. The experiments were generally run on CPU to avoid memory-related distortion of the results, with the exception of SKI, which was run on GPU since it depends on GPU execution for faster MVMs.
We implemented IFF, VFF and spherical harmonics using gpflow (Matthews et al., 2017 which we also used for inducing points, reusing some of the code of Burt et al. (020bb) for the k-DPP (re)initialisation method. The spherical harmonics implementation depends on the backend-agnostic implementation of the basis functions though as noted in the main text, we were unable to produce comparable results. We used the gpytorch implementation of SKI Garddne et al.] 2018b We use the publicly available Tensorflow 2 implementation for B-splines"]
D.1 Real-world dataset information
In all cases, normalise both the inputs and targets to unit mean and standard deviation in each dimension. However, we report the test metrics (RMSE and NLPD) averaged over test points but on the unnormalised scale. The number of training and test points, and the standard deviation of the outputs, for each dataset is reported in Table [[
The precipitation dataset is a regularly gridded (in latitude and longitude) modelled precipitation normals in mm in the contiguous United States for Jamuary 2021 (publicly available with further documentation at https://water weather wwweeer.gvv/preee gow/pprecip/tommlad php) note the data at the source is in inches). We downsample the data by 4 (by 2 in each dimension). The data is highly nonstationary, and so challenging target for GP regression with typically used, usually stationary, covariance functions. In particular, the lengthscales are fairly large across the plains and in the southeast in general, but quite small near the Pacific coast, especially in the Northwest. For a stationary model, the high frequency content in that region leads to globally low lengthscale.
27

---

Published in Transactions on Machine Learning Research (04/2024)
Gaussian kernel
Matérn-3/ 3/2 kernel
1.4
1.2
Normalisod E (&W,
1.0
0.8
0.6
0.4
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1810
1.5
2.0
2.5
Normalised bandwidth (AMe)
Figure 7: As Figure but with half the SNR.
Gaussian kernel
Matérn-3/2 3/2 kernel
1.4
1.2
Normalisof E (
1.0
10
01
1220
990
5555
0.8
0000
0.54
0.6
0.4
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.80
1.5
2.0
2.5
Normalised bandwidth (AMs)
Figure 8: As Figure
2
but with four times the SNR.
28

---

Published in Transactions on Machine Learning Research (04/2024)
Gaussian kernel
Matérn-3
1/rnnn33222 3/2 kernel
1.4
1.2
10
Normalisod E (&W
1.0
$520
0.8
1210
0.6
0.4
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1810
1.5
2.0
2.5
Normalised bandwidth (AMe)
Figure 9: As Figure 2
but with half the lengthscale (twice the bandwidth).
Gaussian kernel
Matérn-3/ 3/ kernel
1.4
1.2
Normalisod E (
1.0
0.8
0.6
0.4
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1810
1.5
2.0
2.5
Normalised bandwidth (AMe)
Figure 10: As Figure
2
but with half the lengthscale (twice the bandwidth) and half the SNR.
29

---

Published in Transactions on Machine Learning Research (04/2024)
Gaussian kernel
Matérn-3 3u2érnn3/2aa kernel
1.4
1.2
Normalisod E (
1.0
0.8
0.6
0.4
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1810
1.5
2.0
2.5
Normalisod bandwidth (AMe)
Figure 11: As Figure
but with half the lengthscale (twice the bandwidth) and four times the SNR.
Gaussian kernel
Matérn-3/2 // kernel
1.4
1.2
Normalisod E (
1.0
0.15
0.8
11
9010
0.6
0.4
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1810
1.5
2.0
2.5
Normalised bandwidth (AMe)
Figure 12: As Figure but with data sampled uniformly on [- W, W./2], and with higher SNR (se 2.067).
30

---

Published in Transactions on Machine Learning Research (04/2024)
Table Summary of real world datasets.
Datasct
Test points
Output standard deviation
Temperature
12 947
3236
2.766
Precipitation
23 144
5785
58.855
House price
106 875
718
0.642
The temperature dataset the change mean land surface temperature (*C). over the year ending February 2021 relative to the base year ending February 1961 (publicly available from https://dataa http://aaaa.gg.s...... nasaa nasa. www..gassoo..oo. gov/ gistemp/mapa) It. is also regularly gridded, over more of the globe.
The house price dataset is snapshot of house prices in England and Wales, which is not regularly gridded. We use a random 20% of the full dataset, and target the log price to compress the dynamic range. It is based on the publicly available UK house price index Inttps://landreiisstyy wwwpiitttrydaaa..... gov. wwt.upkkkppppppiip and we enclose the exact dataset we use.
31


Reviews:
[{"review_id": "S6ys7h6nSW", "reviewer": "Reviewer_bFVe", "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": "I think the overall goal of this paper (achieving $O(M^3)$ runtime for general kernels) seems useful and the proposed approach makes sense at a high level. However, I found the theoretical portions of this paper too unclear to verify the main results. While the topic of the paper is far from my expertise and I wasn't previously familiar with the methods (e.g., SGPR) this work is building upon, I think a substantial amount of the difficulty is due to the writing of the paper itself. Below, I ask a few questions and suggest several ways to improve the presentation.", "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": "Existing methods are either limited to kernels with specific structure (e.g., isotropic or 1D Matern kernels) or have $O(NM^2)$ computational cost, where $N$ is the sample size and $M$ is the number of features used. This work proposes a method that takes $O(M^3)$ (with $M \\ll N$) per optimization step (though it was hard for me to understand the exact assumptions under which this occurs; see below).\n\nIn experiments on synthetic data ($N=10,000$), IFF consistently provides a speedup roughly of $2^5$ over standard sparse GPs (inducing points), in both 1D and 2D. In 1D, performance is similar to VFF in 1D, but significantly better in 2D. In experiments on three somewhat larger ($N=12,947$, $N=23,144$, and $N=106,875$) real datasets, IFF fairly consistently provides among the best trade-off between compute time and test performance, typically beating VFF and B-splines by a small amount.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Major**\n\n1. The assumptions (e.g., disjoint intervals, tail bound, first and second derivative conditions) made in Section 4 are scattered throughout the presentation, making it difficult to understand how the assumptions relate to the results and their proofs. I suggest labeling each assumption when it is first stated (e.g., A1., A2., etc.) and then explicitly referencing each assumption in the result where it is needed (e.g., instead of \"Under the standard assumptions\", say \"Under assumptions A1 and A2\").\n\n2. As far as I can tell, the proofs of the theoretical results (Appendix C) assume the inducing features are a regular grid (Eq. (29)), but this is not clear in the main paper Section 4. In fact, the disjoint interval assumption in the first paragraph of Section 4 seems to suggest that a regular grid is not necessarily assumed. I am quite confused by this.\n\n3. At the beginning of Section 4, it is assumed for simplicity that $D=1$, and the assumptions are stated for this case, but then the results are stated for general $D>=1$. Thus, the reader cannot understand the stated results without referring to the more general assumptions in the Appendix. I suggest sticking to one version ($D=1$ or general $D$) throughout this section.\n\n4. Although the tail bound and second derivative assumptions on the normalized spectral density are fairly standard, I don't think the first derivative assumption (Eq. (17)) is as standard. I therefore suggest (a) adding a sentence discussing this assumption (e.g., whether is is satisfied by commonly used kernels) and (b) explicitly stating the multidimensional generalization of this assumption.\n\n5. I couldn't quite follow how Eq. (33) follows from Eq. (32). Is there a missing \"$1 +$\" or something in Eq. (33)?\n\n6. The bound on $|\\Sigma - \\hat{\\Sigma}|$ in Theorem C.4 needs more details:\n - Why are the elements of $K_{uu}^{-1} \\Sigma_u K_{uu}^{-1}$ of order O(\\epsilon^D)?\n - Some more details about how Theorems C.3 and 4.1 imply the reported convergence rates of the kernel matrices $\\hat K_{u*}$ and $\\hat K_{u*}^*$?\n - \"The other two terms are of the same order $O(M \\epsilon^{2D})$\": Why doesn't this $O(M \\epsilon^{2D})$ doesn't dominate the overall convergence rate of $O(M \\epsilon^{3D})$ in Theorem C.4? Or am I misunderstanding what \"The other two terms\" refers to?\n\n7. At some points, the discussion is too vague to follow. Some examples:\n - After Eq. (15), \"The first equality follows for any variationally correct features (Matthews et al., 2016).\": What does \"variationally correct\" mean?\n - Before Eq. (14), \"we use the result of Burt et al. (2019) to bound\": Which result (they have several)?\nAs far as I could tell, these points aren't clarified in the Appendix.\n\n**Minor**\n\n8. Page 2, First Sentence, \"iterative learning of the variational distribution which is otherwise available in closed form, which results in slower learning overall\": I didn't quite understand why an iterative approach is necessarily slower than a closed form solution. For example, I think it is frequently the case that using SGD to optimize the ridge regression objective is often faster than calculating it in closed form.\n\n9. Page 2, Last Sentence, \"the hyperparameters\", and elsewhere: The paper frequently refers to \"the hyperparameters\", but it's not clear to me what this refers to. Are there hyperparameters other than $\\sigma$? Relatedly, when $\\sigma$ is introduced in Eq. (1), the paper should explicitly point out if this is a hyperparameter.\n\n10. Page 3, Sentence 1, \"the quadratic form $y^⊤Ay$, both of which incur $O(N^3)$ computational cost\": To clarify, the $O(N^3)$ cost is to compute $A$, right? The quadratic form itself should only take $O(N^2)$ times, right?\n\n11. Page 9, Theorem 4.4: I suggest more explicitly writing $\\mu(x_*)$, $\\hat{\\mu}(x_*)$, $\\Sigma(x_*)$, $\\hat{\\Sigma}(x_*)$, to make it clearer that the quantities being bounded are scalars.\n\n12. Figure 1: This should probably appear one page earlier, closer to where it is first referenced (near the top of page 5).\n\n**Typos**\n\nThis paper has *many* minor typos, and needs extensive proofreading. A few examples are listed below:\n\n13. Page 2, Paragraph 2, Last Sentence, \"Existing work... Cunningham et al., 2023)\": This sentence seem to parse grammatically, and I couldn't understand what it was trying to say. Is there an error here?\n\n14. Page 2, Paragraph 3, \"for for\"\n\n15. Section 3, \"to approximate the solution to the linear solve\"\n\n16. Page 5, Paragraph 3, \"$dxdx$\"\n\n17. Page 9, Remark 4.3, \"this is due to the $O(M \\epsilon^3)$ which arise due to...\": It's not clear what quantity is in $O(M \\epsilon^3)$. I guess this should be the error term $E_1 \\in O(M \\epsilon^3)$?\n\n18. Page 11, Last Paragraph, \"For out synthetic experiments\"\n\n19. Page 14, Paragraph 1, \"$N = 10/, 000$\"\n\n20. Page 22, Paragraph 3, \"$M^{1/D} \\in $\" and \"That isve\"\n\n21. Page 22, Eq. (30): I think the definition of $\\square_m$ is intended to be a Cartesian product, but I don't think I have ever seen $\\oplus$ used for the Cartesian product. Should this be, e.g., $\\prod$ or $\\otimes$?\n\n22. Page 22, defs. of $u$, $c_m$, and $\\bar{k}_{m,m'}$: In the multidimensional case discussed here, I think these integrals should be over $\\square_m$ rather than the 1D intervals $[z_m-\\epsilon/2, z_m+\\epsilon/2]$.\n\n23. Page 22, Lemma C.1 (Theorem 4.1): I think \"Theorem 4.1\" should be \"Lemma 4.1\", consistent with Page 8.\n\n24. Page 24, after Eq. (44), \"asumptotically\"\n\n25. Page 24, after Eq. (46), \"the result ??\" (missing reference)\n\n26. Page 25, \"if the spectral density is sufficiently heavy tailed\": I think this should be \"light tailed\" or \"not too heavy tailed\"?\n\n27. Page 25, Theorem C.4: Should \"$O(M^{−2qq+3})$\" be \"$O(M^{−\\frac{2q}{q+3}})$\"?\n\n28. Page 25, after Eq. (63), \"The first term is...\": There is a missing absolute value (\"$|$\") sign in this expression.", "broader_impact_concerns": "I have no concerns regarding ethical implications of this work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "v8Xuca3WU8", "reviewer": "Reviewer_y6ZF", "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\nOverall, I think the paper introduces novel ideas and relevant advances for obtaining scalable approximations for GP regression. I particularly liked how the authors took and presented ideas from previous works in the context that they were interested in before jumping into their main technical propositions. Despite the derivations with Fourier features and inducing methods, I did not find any theoretical/technical typo at first sight, so the proposed equations are likely correct. \n\nThe empirical results look really good and seem to outperform previous approaches (Burt's, VFFs and standard inducing points) in precision/time, which I think is very important to remark.\n\nI also liked that the authors clearly stated the limitations of the method and the gaps that they left for further development or just bc they were out of the scope of the current manuscript. In that direction, I think the work is extremely transparent with the ideas and the justification of the decisions taken, which is important for its future impact and understanding of readers.\n\n**Weaknesses.**\n\nHaving stated the main strengths, and once again remarking that I think is a paper of high-quality ideas --- I think is worth saying that the manuscript should be updated to make such ideas and contributions \"shine\" in the text. In that regard, the paper loses clarity and some details are omitted without any reason at some points of the story. It is very clear while reading that after developing the project, some ideas are extremely clear for the writer, and then they are partially omitted or not clearly explained. These sorts of issues are of course possible to be fixed in this reviewing round, so I encourage the authors to improve the manuscript to make the paper stronger and of better clarity for its future impact in the community. For helping in this process, I will add some notes and points that I think should be improved in the next section.", "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 work proposes integrated Fourier features to perform fast-and-scalable GP regression in O(M^3) with stationary kernels and of regular spectral density (i.e. non-periodic kernels). The proposed ideas build on top of inter-domain GPs and recent methods in the context of frequency-based representations of kernels to obtain scalable sparse approximations for GPs. While following this line of thinking, the propositions show a good performance on the experimental results on both synthetic and real-world datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Some changes I would like to see updated to vote for a strong acceptance of the paper and likely for improving its clarity and impact in the future.\n\n**Point 1.** Page (1) --- \"the scaling with N is problematic ....\" I see the point of saying this and reasoning that scaling with N and using stochastic optimization makes the learning process slower in general. However, if the paper needs this sort of thinking for avoiding the use/introduction of stochastic gradient descent methods, I think more attention should be put on the computational advantage of the method, like really showing the time-difference between using Fourier features with sparse approximations and other stochastic gradient-based methods. I basically say this bc many readers are now familiar and native in the use of SGD, even with GPs, and stating that the proposed ideas are much more faster, even showing clear empirical results would convince many of them to believe this is a good direction. Indeed, this last sort of things would be better than the last paragraph.\n\n**Point 2.** Fixing writing typos. I caught quite a few writing typos like \"pesudo\" instead of pseudo or \"for for\" just in the first pages. A deep look into the text and some rewording/rewriting of such typos would help.\n\n**Point 3.** The introduction is great in general, but when one jumps in the background things get a bit messy. So it took me quite a few reads to remember/realize the key point of everything, which is that Kuu is diagonal because the Fourier features are independent. Of course, this in the end is what makes the computation to be reduced to O(M^3) bc the trace can be reduced down to the equation written in Appendix A. This point is kind of hidden in the text, as one realizes it when looking into Appendix A. Just notice that this Kuu is usually the GP prior on standard inducing points, so any practitioner might accidentally think that Kuu is still fully correlated and then have issues to understand what is going on. So I think is important to remark and remind all around the paper that this Kuu is diagonal and features independent in all cases.\n\n**Point 4.** Could the notation be matched with D. Burt's papers or VFF? I think is not a lot of work, (i.e. using bold letters for inducing features and matrices for instance and maybe using more standard acronyms for frequency domains), and it would save a lot of thinking-time for the reader.\n\n**Point 5.** Some lingering questions in my mind that maybe are obvious but could be perhaps clarified. The issues with the Dirac delta: a) unsuitable for constructing the conditional prior (this needs more explanations and at least to be smoother for the reader), and b) difficulties to integrated out the Dirac delta in Eq. (11) --- is this bc it cannot be done, not numerical methods available, intractability, difficult functional evaluation, etc... the reasons are not clear to me here and i think is super important.\n\n**Point 6.** Spectral density, c and k. These three elements are the ones presented and indexed at every step of the paper. Indeed, the different methods have different indices 1,2,3.. However, this is a bit difficult to follow and not very clear. Perhaps, it would help to cleary state the role of every one of them, and also use other way to say that they are different approximations or methods.\n\n**Point 7.** I see the utility of averaging the Fourier features over disjoint intervals of fixed width. But as said before, this is very related to the problems on integrating out the Dirac delta and so on. Could it be explained a bit more why the averaging wants disjoint intervals, if this ones are exactly near each other, the size of \\epsilon chosen, etc...?\n\n**Point 8.** The transformation of the features by the invertible linear transformation T seems useful, but I do not understand what is going on and how it helps to have only one approximation... This part is super difficult to follow for me, and also do not understand why the normalization by spectral density turns into the normalization by square root. I guess bc T is used, but not clear how to me.\n\n**Point 9.** Why only for D<4? I see that it performs poorly for larger D, but where is the bottleneck or issue? I would also add a bit more of analysis/detailed context in the results, and connecting things a bit better with the Appendix. For instance, explaining a bit the nature of the real-world datasets and the characteristics of the data to easily understand the impact of the performance.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "DJo6gEXf3W", "reviewer": "Reviewer_ZPJr", "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:\n1. The proposed features are motivated by a convergence analysis and empirical exploration.\n2. The paper demonstrates practical speedup in synthetic and real-world spatial regression tasks.\n\nWeakness:\n1. The paper's experiments involve synthetic datasets with arbitrarily chosen parameters, such as a signal-to-noise ratio (SNR) of 0.774, which may not be representative of real-world scenarios.\n2. The paper mentions inducing points with inducing inputs optimized taking longer than other methods, but it does not delve deeper into the implications of this observation.", "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 introduces integrated Fourier features for Gaussian processes (GPs) to enhance their scalability for larger datasets. The method extends performance benefits to a broad class of stationary covariance functions. Overall speaking, this paper proposed a speedup method for Gaussian process applied on large dataset. Experiments show that the proposed IFF method can be applied on a various source of dataset and show great scalability, reducing the  computational costs especially when M << N.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I think there is nothing much to change. Maybe it will be good to have a table to compare the computational costs and other properties of current methods for solving Gaussian process problem.", "broader_impact_concerns": "There doesn't seem to be any direct ethical concerns related to the paper's content.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper presents a method for speeding up Gaussian process modeling of spatial data, using what they call integrated variational features.  Previously there were methods that took advantage of the spatial structure of the data to reduce the computational complexity of GP inference from N^3 to M^3, where N and M are the number of examples and number of feature dimensions respectively.  However, those methods were restricted to a small restricted set of covariance functions, limiting the modeling flexibility.  This work extends the M^3 complexity to many more covariance functions, using the variational Fourier features.  This would allow for much more flexible modeling of large spatial datasets.

All three reviewers voted for accept (two leaning, one accept).  The reviewers found that the paper was well motivated, addressed an important and relevant problem, the methods were well justified and empirical results convincing.  One reviewer found that the theoretical results were somewhat confusing and unverifiable.  However, after reading over the author response, that reviewer now seems satisfied.  The reviewers all found that the claims made in the paper were well supported by the empirical and theoretical evidence.  Overall, this is a strong paper that advances the state-of-the-art in efficient spatial modeling using Gaussian processes and seems like a good candidate for acceptance to TMLR.  Therefore I recommend acceptance.
Decision: Accept (poster)