Gauss-Legendre Features for Gaussian Process Regression

Gaussian processes provide a powerful probabilistic kernel learning framework, which allows learning high quality nonparametric regression models via methods such as Gaussian process regression. Nevertheless, the learning phase of Gaussian process regression requires massive computations which are not realistic for large datasets. In this paper, we present a Gauss-Legendre quadrature based approach for scaling up Gaussian process regression via a low rank approximation of the kernel matrix. We utilize the structure of the low rank approximation to achieve effective hyperparameter learning, training and prediction. Our method is very much inspired by the well-known random Fourier features approach, which also builds low-rank approximations via numerical integration. However, our method is capable of generating high quality approximation to the kernel using an amount of features which is poly-logarithmic in the number of training points, while similar guarantees will require an amount that is at the very least linear in the number of training points when random Fourier features. Furthermore, the structure of the low-rank approximation that our method builds is subtly different from the one generated by random Fourier features, and this enables much more efficient hyperparameter learning. The utility of our method for learning with low-dimensional datasets is demonstrated using numerical experiments.


Introduction
Gaussian processes (GPs) [41] provide a powerful probabilistic kernel learning framework, which allows learning high quality nonparametric regression models via methods such as Gaussian process regression (GPR). Indeed, GP based methods are widely used in machine learning and statistics. They have been applied to a wide variety of problems, such as data visualization, Bayesian optimization [28], modeling dynamics and spatial data analysis [30]. One of the key advantages of the GP formulation of kernel regression is that the marginal likelihood is a function of the kernel hyperparameters, and that it can be computed via a closed-form formula. By maximizing the marginal likelihood, one can learn the hyperparameters from the data, thereby tuning the method in a principled manner.
However, learning GPs comes with an hefty computational price-tag. Given a training set of n points of dimension d, exact GPR requires solving a (usually dense) linear equation, and thus requires O(n 3 ) FLOPs. Prediction costs O(nd) FLOPs per test point. Such costs are problematic for datasets with more than a few thousand points. The situation is even more severe if we consider the hyperparameter learning phase: here the cost is O(n 3 ) FLOPs per hyperparameter in a learning iteration (assuming we use a first-order optimization method). Hyperparameter learning of exact kernel models on large-scale data is even more unrealistic than training such models.
Given the ubiquity of GPs, it is unsurprising that there is a rich literature on scaling GP-based method, e.g. [23] and [42]. One attractive approach is to approximate the kernel matrix (also known as covariance matrix) K θ as a sum of a diagonal matrix (often a multiple of the identity) and a low rank matrix [31]: In the above, K θ ∈ R n×n denotes the kernel matrix, where the subscript θ denotes the dependence of the kernel matrix on the hyperparameters θ (discussion of our notation appears in Section 2.1), Z has s n columns, and W(θ), D(θ) are diagonal matrices. The various steps of GPR can be much more efficiently conducted on a kernel whose kernel matrix has the structure of the righthand side of Eq. (1), e.g. training takes O(ns 2 ) (see Section (4.2) for details on efficient GPR with low-rank approximations with an even more restricted structure in which D(θ) is a multiple of the identity).
In the kernel learning literature, methods for forming a low rank approximation of kernels can be roughly split into two approaches: methods that use data-dependent basis functions, and methods that use independent basis functions. An example for the first kind is the Nyström method [40]. Such methods utilize the given training data, and thus may outperform methods that use independent basis functions, especially when there is a large gap in the eigenspectrum. However, data dependence can incur additional costs. For example, the Nyström method requires keeping some of the data as part of the model.
Another class of methods for building low rank approximations of kernel matrices are methods that use independent basis functions, and thus approximate the kernel function directly. One such important and highly influential method is the random Fourier features approach suggested by Rahimi and Recht in 2007 [24]. Following the publication of [24], there has been extensive research on random features, including works that attempt to improve the approximation quality of the method (e.g., [33,9]), works that focused on using in random features to learn huge datasets (e.g. [16,2]), and works that focused on theoretical analysis of random features (e.g. [44,29,1]). The previous list is far from exhaustive. In the context of our work, worth mentioning is [1] which showed that if the kernel matrix of the approximate kernel spectrally approximates the kernel matrix of the true kernel then the excess risk when using kernel ridge regression with the approximate kernel is not much larger than the excess risk when using the true kernel.
Random Fourier features, and random features methods in general, are based on writing the kernel function as an integral and then using numerical integration schemes in order to construct a low rank approximation of that function 1 . In random Fourier features, a shift-invariant kernel is rewritten as an integral via an application of Bochner's theorem, and Monte-Carlo integration is used to build the low rank approximation. The use of Quasi Monte-Carlo in lieu of Monte-Carlo integration was explored in [3]. Bach explored the connection between random Fourier features and kernel quadrature rules in [4], however without providing any practically useful explicit mappings for kernels. Monte-Carlo and Quasi-Monte Carlo integration admit only slow convergence rate. As a consequence, the number of features required for spectral approximation when using Monte-Carlo or Quasi Monte-Carlo integration must be polynomial in quality parameter of the spectral approximation. In this paper we argue that in the context of Gaussian process regression a stronger notion of spectral equivalence is required. The slow convergence rate of Monte-Carlo or Quasi Monte-Carlo integration implies that at best the number of features required for spectral equivalence is linear in the number of training points, which is obviously undesirable.
Random features approaches based on Monte-Carlo and Quasi Monte-Carlo suffer from another serious defect when it comes to GPR: the low rank approximation they build has the form While for training and prediction, this structure works equally as well as the structure in Eq. (1), when it comes to hyperparameter learning this is no longer the case; see Section 4.3. One can construct faster converging low-rank approximations using numerical quadrature rules such as Gaussian quadrature. Dao et al. considered the use of Gaussian quadrature in the context of kernel learning [12]. Gaussian quadrature rule is a method for numerically approximating weighted integrals (i.e., integrals of the form´∞ −∞ f (x)w(x)dx where w(x) ≥ 0 is a weight function) that is optimal in some formal sense. In the context of approximating kernel functions, the weight function w(·) is determined by the kernel function and the value of the hyperparameters. Once the weight function has been determined, in 1 A rank k bivariate function f (x, y) is a function that can be written as f (x, y) = k j=1 σ j φ j (x)ψ j (y) for some σ 1 , . . . , σ k , φ 1 , . . . , φ k , ψ 1 , . . . , ψ k [35]. order to use a Gaussian quadrature the nodes and weights corresponding to that particular weight function must be computed. Efficient algorithms exist, but these algorithms require the computation of integrals as well. For a single kernel, that is when using a fixed value of the hyperparameters, and when using a fixed number of quadrature features, computing the nodes and weights is a one-time offline task. However, if the hyperparameters are not fixed, e.g. when they are set using hyperparameter learning, Gaussian quadrature becomes unrealistic. Furthermore, the fact that the nodes and weights change with the hyperparameters implies that we must use an approximation of the form of Eq. (2) and not of Eq. (1), which is less desirable. The connection between random Fourier features and quadrature rules was also explored in [21].
Low rank approximations for kernels matrices have also been widely used in the statistics literature, and in particular the spatial statistics literature [11,14,5,15,17]. Possible limitations of the low rank approximation approach in the context of spatial statistics have been noted in [23,5,32,26], and analyzed mathematically in [31]. The use of random features in the context of spatial statistics was explored in [34].
In this paper we propose a quadrature based low-rank approximation approach for efficient GPR involving a wide class of kernels which includes shift-invariant kernels (i.e., stationary covariance functions). Unlike previous literature which uses quadrature features in the context of GPR, our method forms an approximation of the form of Eq. (1), and so is able to efficiently perform hyperparameter learning in addition to training and prediction. Our method achieves this by using a fixed set of quadrature nodes and weights, and designing the approximation so that varying the hyperparameters corresponds to only changing the integrand.
Specifically, our method uses a Gauss-Legendre quadrature. Gauss-Legendre quadrature is a Gaussian quadrature for the uniform weight function on a finite interval. Thus, the weight function does not change with the hyperparameters, and with it the quadrature nodes and weight stay fixed, whereas only the integrand varies. Changing only the integrand translates to a simplified parametric form for the approximate kernel matrix (Eq. (1)) which is more amenable to efficient computations. Our proposed method, which we call Gauss-Legendre Features, is described in Section 4.
The Gauss-Legendre quadrature is designed to approximate integrals with an integration area which is a finite interval. However, for most widely-used kernels the integrand has infinite support. We address this issue by utilizing the fact that for such kernels the integrand decays quickly, so we can approximate the integral by truncating the integration area. The truncation cutoff is determined by a parameter of our method. Another parameter is the number of features (i.e., quadrature nodes) used in the approximation. In order to set these two parameters correctly, we need a method for assessing the quality of one kernel function approximation by another. To that end, we introduce the notion of spectral equivalence, and argue that if one parameterized family of kernels is spectrally equivalent to another one, then that first family is a good surrogate for the second family in the context of GPR. These results are summarized in Section 3.
We rigorously analyze how to set the truncation cutoff and the number of features to achieve spectral equivalence (these results are reported in Section 5). Here another advantage of using Gauss-Legendre quadrature becomes evident: the Gauss-Legendre quadrature converges much faster the Monte-Carlo or Quasi Monte-Carlo integration, so typically the number of features is sublinear in the training size. Indeed, for widely used kernels like the Gaussian kernel and the Matèrn kernel, the number of features required when using Gauss-Legendre features is poly-logarithmic in the training size (see Section 6). Sublinear number of features is also likely achievable using Gaussian quadrature (kernel learning using Gaussian quadrature is suggested in [12], however without proving spectral equivalence). Yet, as explained this is rather problematic for hyperparameter learning, and in general requires a large overhead for computing the quadrature nodes and weights.
Finally, empirical results (Section 7) clearly demonstrate the superiority of our proposed method over classical random Fourier features when conducting Gaussian process regression on low-dimensional datasets.

Notation and Basic Definitions
We consider all vectors as column vectors, unless otherwise stated. For a vector x or a matrix A, the notation x * or A * denotes the Hermitian transpose. The n × n identity matrix is denoted by I n . A Hermitian matrix A is positive semidefinite (PSD) if x * Ax ≥ 0 for every vector x. Also, for any Hermitian matrices A, B of the same size, the notation A B means that B − A is PSD.
We consider n pairs of training data (x 1 , y 1 ), . . . , (x n , y n ) ∈ X × Y ⊂ R d × R, where x denotes the input vector of dimension d and y denotes a scalar response. A kernel function (aka covariance function) is a function k : X × X → R which is positive definite, i.e. for every m ∈ N and x 1 , . . . , x m ∈ R d , the matrix K ∈ R m×m defined by K ij = k(x i , x j ) is PSD. The matrix K is known by various names: kernel matrix, Gram matrix, covariance matrix. Given a dataset x 1 , . . . , x n ∈ R d , we will conveniently use X to denote the n-by-d matrix whose rows are x 1 , . . . , x n , and use K(X, X) to denote the kernel matrix corresponding to the kernel k with data X. For another kernelk we will useK(X, X) to denote the kernel matrix.
In many cases we will deal with parameterized families of kernels {k θ } θ∈Θ , where θ represents the hyperparameters vector, and Θ is a set of possible parameters values. The kernel matrix corresponding to k θ is denoted by K θ (X, X). We also group the responses y 1 , . . . , y n into a single vector y ∈ R n .
The Kullback-Leibler divergence (abbreviated KL-divergence henceforth) is a well established metric for how much one distribution is different from a reference distribution. We denote the KL-divergence between two probability distributions P on Q by D KL (P, Q), and recall the following is a well established result 2 : if

Gaussian Process Regression
Gaussian Process Regression (GPR) is a Bayesian nonparametric approach for regression. First, the following regression model is assumed: n is a (hyper)parameter). Additionally, it is assumed that f is a Gaussian Process, f (x) ∼ GP (µ (x) , k(x, x )), where k is the kernel function. This means that for any set of data points x 1 , . . . , x m ∈ R d the vector f ∈ R m defined by f j = f (x j ) (j = 1, . . . , m) is a Gaussian random vector with mean defined by µ(X) = [µ(x 1 ) · · · µ(x m )] T and covariance matrix K(Z, Z). Throughout the paper we assume, for the sake of simplicity, that the mean function µ(x) is 0. This simplifies the formulas while not really restricting generality (a nonzero mean can be easily handled). Under these assumptions y ∼ N (0, K(X, X) + σ 2 n I). Under these priors, the expected predictive value for f (x) at a test x is Consequently, training is conducted by computing the vector α := K(X, X) + σ 2 n I n −1 y From these formulas we see that assuming that evaluating the kernel function takes O(d) operations and that α is computed using direct factorization, training takes O(n 3 ) operations and prediction takes O(nd) operations. The previous description is for a fixed kernel k. Typically, the kernel k θ has hyperparameters which we represent throughout the paper by the vector θ. The hyperparameters are usually constrainted to some possible set of hyperparameters values Θ, and thus defines a parameterized family of kernels {k θ } θ∈Θ . Hyperparameter learning refers to the process of determining the value of the hyperparameters directly from the training data, and is considered one of the important advantages of the GP framework. This is typically conducted by maximizing the log marginal likelihood: In order to maximize L(θ) using a first-order optimization method it is required to compute its gradients. Using direct methods, computing the gradient takes O(n 3 |θ|) where |θ| represents the number of hyperparameters in θ.

Random Fourier Features
Random Fourier Features (RFF) [24], is one of the most popular methods for constructing a low rank approximation of kernels and scaling up kernel methods. The method targets shift-invariant kernels, i.e. kernels of the form k(x, x ) = k 0 (x − x ) for a positive definite function k 0 (·). RFF is motivated by a simple consequence of Bochner's Theorem: for every shift-invariant kernel for which k 0 (0) = σ 2 f there is a probability measure µ and possibly a corresponding probability density function p(·), both on R d , such that Let us assume that the density p(·) exist. If one chooses η 1 , . . . , η s randomly according to p(·), and defines So, an approximated kernel can be defined: The kernel matrix corresponding to the approximate kernel is where Z ∈ C n×s to be the matrix whose m th row is ϕ(x m ) * . The low rank structure ofK (RFF) (X, X) allows more efficient training (O(ns 2 )) and predictions (O(sd)), which are attractive if s n. The previous description is for a fixed kernel (and fixed hyperparameters). When using GPR with hyperparameter learning we are dealing parameterized family of kernels {k θ } θ∈Θ . This case has not been considered in Rahimi and Recht original work [24]. We discuss it in Section 4.3.

Spectrally Equivalent Kernel Approximations
Our strategy for scaling up GPR is based on approximating the kernel k θ by an approximate kernelk θ that is low-rank in some sense which will become apparent in the next section. This raises the question: how can we determine whetherk θ indeed approximates k θ well? In [1], the authors suggested that in the context of kernel ridge regression, spectral approximations of the kernel matrices allows us to reason about how well one kernel is approximated by another. The argument in [1] is based on risk bounds for given fixed hyperparameters, and so is less appropriate for GPR where hyperparameter learning is common practice. In this section, we introduce the notion spectral equivalence, a stronger form of spectral approximation, and connect it to hyperparameter learning in GPR.
Assume a bounding set X ⊆ R d for the data, then given a dataset (x 1 , y 1 ), . . . , (x n , y n ) ∈ X × R, the general assumption when using a kernel k is that If, however, we would have used the kernelk, then the assumption would have been Thus, a measure on how wellk approximates k might be devised by measuring how much N (µ,K(X, X)) is different from N (µ, K(X, X)). The KL-divergence is a well-established measure on how different one probability distribution is from a reference distribution, so arguably,k approximates k well if the KL-divergence D KL N (µ, K(X, X)), N (µ,K(X, X)) is small. Indeed, the use D KL N (µ, K(X, X)), N (µ,K(X, X)) as such a measure was suggested in the literature on spatial data analysis [5,25,31].
The notion of spectral equivalence, which we develop below, is a measure on how two matrices are close to one another. To connect it to the KL-divergence, which we use to measure how wellk approximates k, we have the following lemma, which implies that if the covariance matrices of two multivariate distributions are close, then the KL-divergence is small. Lemma 1. Suppose that µ ∈ R n and Σ 0, , Σ 1 ∈ R n×n are two symmetric positive definite matrices. Suppose that We first need the following Lemma.
Proof of Lemma 1. Since for two symmetric positive definite matrices A and B, A B implies B −1 A −1 , Eq. (5) implies that Multiplying by Σ 1/2 0 on the right and left sides gives . Thus, Also, from Lemma 2, there exist γ 1 , . . . , γ n ∈ [−n −1 , n −1 ] such that Using Eq. (3), we obtain Lemma 1 motivates the following definitions: Definition 3. We say that a n-by-n symmetric matrix A is spectrally equivalent to another n-by-n symmetric Definition 4. Let n ≥ 1 be an integer, and X be a data domain. Two positive definite kernels k andk are n-spectrally equivalent on domain X if for every X with n rows in X , the kernel matrixK(X, X) is spectrally equivalent to the kernel matrix K(X, X).
The last definition uses two specific kernels, k andk. In GPR, we usually use a parameterized family of kernels {k θ } θ∈Θ , where θ represents the hyperparameters, and Θ is a set of possible parameter values. We generally assume that Θ is bounded. Boundedness of Θ is necessary, since without it, it is possible to drive the kernel matrix to identity, thereby making it impossible to approximate it using a low rank matrix. We then approximate each kernel k θ byk θ , that is we use the parameterized family {k θ } θ∈Θ . We say that the parameterized family {k θ } θ∈Θ approximates the parameterized family {k θ } θ∈Θ well if for every θ ∈ Θ the kernelk θ approximates k θ well, as is captured by the following definition.
Definition 5. Two parameterized families of positive definite kernels {k θ } θ∈Θ and {k θ } θ∈Θ are said to be n-spectrally equivalent on domain X if for every θ ∈ Θ, k θ andk θ are n-spectrally equivalent over X .
In light of Lemma 1, if two parameterized families {k θ } θ∈Θ and {k θ } θ∈Θ are n-spectrally equivalent, then for any parameters θ and any dataset consisting of n data points, the distributions on the response assumed by the two GP models induced by these families are close in the sense that the KL-divergence is close to 1.

Gauss-Legendre Features
In this section, we present our proposed method (Gauss-Legendre Features), and show how it can be used to perform efficient Gaussian process regression. Our method includes two important parameter vectors: U and s. In the next section we show how these parameters can be set in order to obtain an approximation that is spectrally equivalent to the true kernel.

Feature Map
We begin by describing the Gauss-Legendre feature map. The proposed method builds feature maps for kernel families that can be written in the following form: In the above, the function ϕ : X ×R d → C is such that for every x ∈ X the function ϕ(x, ·) is even-symmetric (i.e., for every , and for every θ 0 the function p(·; θ 0 ) is an even probability density on R d . Note that θ 0 can be a vector.
Note that in Eq. (8) we included a ridge term σ 2 n γ(x − x ). Typically, the ridge term is omitted from the kernel but appears in various equations involving the kernel matrix due to Gaussian noise assumption in the GPR model. While we could state our theory in the more traditional way of having the noise term outside of the kernel, the definitions and theorems statements will be somewhat more cumbersome. We found it more convenient to include σ 2 n as part of the vector of the parameter set θ, and include the ridge term in the kernel definition. The resulting equations are the same.
There are quite a few kernel families that adhere to this structure. For example, due to Bochener's theorem, a shift-invariant kernel with an additional noise level term can be written in the form So, we can use ϕ(x, η) = e −ix T η to cast shift invariants kernels in the form of Eq. (8).
The underlying idea of Gauss-Legendre features is to first truncate the integral Eq. (9) to the domain , for some U = (U 1 , . . . , U d ) T , and then approximate the truncated integral using a tensorized Gauss-Legendre quadrature. The domain Q U might depend on X and Θ, but not on the concrete dataset x 1 , . . . , x n . Let (χ whereη . . . In Eq. (10), we assumed we have a bijective mapping a 1 , . . . , a s between {1, . . . , s} and {1, . . . , s 1 } × · · · × {1, . . . , s d } and then defined: Finally, the parameterized family of approximate kernels is Note that the conditions that p(·; θ 0 ) is even and ϕ(x, ·) is even-symmetric, coupled with the fact that the Gauss-Legendre quadrature is symmetric, ensures thatk θ (x, x ) is always real. Consider the first term in the right hand side of Eq. (11). It is a bivariate function which can be written as a sum of s separable bivariate functions. Thus, we can informally view s as the rank of the decomposition, and if s is small, then this is a low-rank approximation. The parameterized family of approximate kernels {k θ } θ∈Θ is composed of kernels that can be written as a low-rank bivariate function plus a ridge term. In the next subsection we show how to utilize this low-rank structure in order efficiently perform Gaussian process regression.
Of course, the crucial question is how do we choose U = (U 1 , . . . , U d ) and s = (s 1 , . . . , s d ). We want to choose these parameters such that {k θ } θ∈Θ and {k θ } θ∈Θ are n-spectrally equivalent over the domain X , where n is the target dataset size (since GPR is nonparametric, the effective rank of the kernel matrix goes to infinity when n goes to infinity, so it is impossible to approximate the kernel matrix well with a matrix of fixed rank, i.e., with s fixed, as n goes to infinity). We discuss this question in the next section. In the reminder of this section, we discuss how to efficiently perform GPR using the approximate kernel family {k θ } θ∈Θ .

Efficient Gaussian Process Regression
As a first and crucial step, we show how to write the kernel matrix ofk θ as a low-rank matrix plus a ridge term. Given a dataset x 1 , . . . , x n ∈ X , let X, as usual, denote the matrix whose row j is x T j . Define the matrix Z ∈ C n×s , Z lj := ϕ(x l , η j ) .
Note that Z depends on U and s, but not on the hyperparameters θ. Next, define We now haveK Notice that dependence on θ is confined to the diagonal matrix W(θ). This will be very helpful in deriving efficient formulas for GPR.
We now discuss each of the various stages of GPR separately. For simplicity, we assume that the GP prior has zero mean (µ = 0).
Training. Given y 1 , . . . , y n , training usually amounts to computing the vector where y = [y 1 , . . . , y n ] T . However, in our case, in order to utilize the structure ofK θ (X, X), we instead compute: In the above, the second equality is a simple consequence of the Woodbury matrix identity. Since W(θ) has positive diagonal, w can be computed using O(ns 2 ) operations, discounting the cost of computing Z.

Prediction. Given a test set x
which are distinct from the training set, the predicted vector Hence, once we have w (computed during training), we can compute y (t) using O(ts) operations, discounting the cost of compute Z (t) .
Hyperparameter Learning. Hyperparameter learning amounts to finding the hyperparameters θ which maximize the log marginal likelihood. To do so, we need to be to able to efficiently compute the log marginal likelihood, and its gradient. It is well known that the likelihood is given by and the derivatives are given by where θ i represents an hyperparameter in θ.

Proposition 6.
After an O(ns 2 ) preprocessing step of computing Z * Z, and discounting the cost of computing the partial derivatives of p with respect to the hyperparamters, the log marginal likelihood L(θ) and the gradient ∇L(θ) can be computed in O(ns + s 3 + s|θ|) arithmetic operations, where |θ| represents the number of hyperparameters in θ. Furthermore, the amount of memory storage required is O(s 2 ).
Proof. First, let us consider the computation of the likelihood. For the first term in Eq. (12), note that In the previous equations, we made the dependence of α and w on θ explicit. Obviously, once Z * Z has been computed (an O(ns 2 ) preprocessing step), we can compute both w(θ) and α(θ) in O(ns + s 3 ). The first term in Eq. (12) is now equal to −y T α(θ)/2.
For the second term in Eq. (12), using the matrix determinant lemma, we have Since W(θ) is diagonal, once Z * Z has been computed, we can compute this term in O(s 3 ) operations.
Next, let us consider the computation of each derivative of the likelihood according to Eq. (13). The crucial observations are: (using θ i to denote an hyperparameter in θ 0 ) where using last equality amounts to computing partial derivative ∂p(η j ; θ 0 )/∂θ i for j = 1, . . . , s. For the second term in Eq. (13) we have, so this term can be computed in O(s) operations once we compute w(θ) and α(θ) (which are computed during the computation of the likelihood). For the first term in Eq. (13), let Again, once Z * Z has been computed, F(θ) can be computed in O(s 3 ) operations. Now, using the Woodbury formula and cyclicality of the trace, we have (the calculations leading to these formulas appear in Appendix A). Thus, once Z * Z, F(θ) and the diagonal of Z * ZF(θ) have been computed, all these computations can be done in O(s|θ|). Note that the diagonal of Z * ZF(θ) can be computed using O(s 2 ) operations once we have Z * Z and F(θ).
In terms of memory storage, notice that once Z * Z and Z * y have been computed there is no longer any need for Z (which requires O(ns) words to store). However, in order to compute Z * Z and Z * y we do not need to form all of Z in memory, but rather can stream over the training set transforming every training point using ϕ and accumulating its contribution to Z * Z and Z * y. Thus, the dominant storage cost is for The computations can be performed more stably by utilizing various matrix identities. We delegate the details to Appendix A.
To have a computational advantage in the training and prediction steps we need s = o(n). However, since for most kernels computing the gradient of the likelihood requires O(n 3 |θ|), our method has a computational advantage in the hyperparameter learning phase even if s = Θ(n).

Comparison to Other Methods
Our proposed method is very much inspired by the Random Fourier Features (RFF) method [24]. Although originally defined only for shift-invariant kernels, the method can be easily generalized for kernels of the form of Eq. (8). We refer to the generalization as Monte-Carlo Features (MCF). RFF is a special case of MCF. In particular, given some fixed parameters θ, a MCF approximate kernel is where η 1 , . . . , η s are sampled from the density function p(η; θ 0 ). As before, we include the ridge term σ 2 n γ(x, x ) in the kernel definition. Thus, the kernel matrix approximation is Obviously, training and prediction can be efficiently executed by utilizing the identity plus low-rank structure ofK (MCF) θ much in the same way as we have done for Gauss-Legendre features, and indeed this is the reason the method was developed [24].
The above developments were for a fixed θ. However, it is less clear how to define a family of approximations for various θ, and perform hyperparameter learning. A key issue is that the kernel approximation should vary smoothly with θ, so obviously fresh η 1 , . . . , η s cannot be sampled differently for every θ. It is outside the scope of this paper to consider how to use MCF to define parameterized families of kernel approximation suitable for hyperparameter learning. Nevertheless, since we wish to use RFF as a baseline for complexity comparisons and numerical experiments, we show how it is possible to use the specific case of RFF to form parameterized families of kernel approximation and perform hyperparameter learning for a restricted family of kernels that includes the Gaussian and Matèrn kernels.
Specifically, we will consider a restricted class of shift-invariant whose kernel has the following specific form: where k 0 (·) is a positive definite function, and L is diagonal with positive entries L 1 , . . . , L d which are part of parameter vector θ (i.e., θ = [L 1 , . . . L d , σ 2 f , σ 2 n ]). Note that the number of variables in θ 0 is equal to the dimension d. Gaussian and Matèrn kernels are examples of such kernels. In this case we can write k θ in the form of Eq. (8), where ϕ(x, η) = e −ix T η . We further assume that p(η; L) is such that sampling a random vector η is the same as sampling from the distribution defined by p(·; I d ) and scaling the vector by L −1 . Thus, for such kernels we can sample η 1 , . . . , η s once from p(·; I d ) and view any change in θ 0 = L as corresponding change in η 1 , . . . , η s . Concretely, letting However, the crucial point is that now Z(L) depends smoothly on L, so we can compute gradients. Formulas quite similar to the ones derived in the previous section can be derived (we omit most details), with the main difference being in taking the derivative of the kernel matrix with respect to the parameters in L. Here we have ∂K In terms of complexity, the main difference between Gauss-Legendre Features and RFF is that for the former the matrix Z * Z stays constant when θ varies, and so the product can be computed once, while for the latter Z(L) * Z(L) varies, and so changes every iteration. This adds an additional cost of O(ns 2 ) operations for every gradient computation. Furthermore, we need to compute ∂Z(L)/∂L k for k = 1, . . . ., d in each iteration, each costing O(nsd) operations, for a total of O(nsd 2 ) operations. Furthermore, since Z(L) changes in each iteration, and Z(L) features in many of the equations, we cannot compute the matrix Z(L) * Z(L) once and reduce storage costs to O(s 2 ), and using Z(L) implicitly many times will incur a large overhead. Thus, for RFF the storage cost is O(ns), Table 1: Computational complexities (arithmetic operations) comparison between Gauss-Legendre Features and Random Fourier Features for kernels of the form of Eq. (15) (e.g., non-isotropic Gaussian an Matèrn kernels). In the table, n is the size of the training set, t is the size of the test set, s is the approximation rank, and I is the number of gradient computations for hyperparameter learning (e.g., number of gradient descent iterations).

Gauss-Legendre Features
Random Fourier Features Hyperparameter learning: A similar issue will likely arise when using features based on Gaussian quadrature, like was suggested in [12] (that paper does not discusses GPR hyperparameter learning). When θ changes, the distribution that defines Gaussian quadrature changes. Unlike Guass-Legendre features which uses fixed nodes, for Gaussian quadrature features the quadrature nodes change with θ, which prevents the use of a fixed feature matrix Z(θ). Furthermore, when performing hyperparameter learning with Gaussian quadrature features, we need to not only compute the quadrature weights but also compute their derivatives.

Parameter Computation
In order to complete the description of our method, we need to specify how to choose U and s. First, we show how to set U and s for a fixed θ ∈ Θ and n such that k θ andk θ are n-spectrally equivalent. We then consider how to set a fixed U and s such that the two families {k θ } θ∈Θ and {k θ } θ∈Θ are n-spectrally equivalent.

From Matrix Approximation to Integral Approximation
For now (and until subsection 5.4) let us assume that θ is fixed. Our goal is to set U and s such that k θ and k θ are n-spectrally equivalent, i.e. for every dataset x 1 , . . . , x n ∈ X In other words, we want to set U and s such that for every v ∈ R n , Henceforth, for conciseness, we drop X from the expressions, although the various expressions implicitly depend on X. Let Since rescaling v rescales all the terms in the previous inequality, we can assume without loss of generality that v T K θ v = 1. In that case, Eq. (16) is equivalent to Thus, the nodes η 1 , . . . , η s and weights h 1 (θ 0 ) . . . , h s (θ 0 ) function as a quadrature approximation.

Truncating the Integral
As alluded earlier, we approach the quadrature approximation Eq. (17) by first truncating the integral and then using a Gauss-Legendre quadrature for the truncated integral. In other words, we write We set U such that the first term is smaller than σ −2 f n −1 /2, and set each of the components in s to be large enough so that the second term is also smaller than σ −2 f n −1 /2. Obviously, we want to set the components in U to be as small as possible, to limit the integration area. Having a smaller integration area allows us to use smaller values in s. The minimal values in U such that the first term is bounded by σ −2 f n −1 /2 depends on how quickly p(η; θ 0 ) decays as η ∞ → ∞: the faster the density decays, the smaller is the region where the function value has significant contribution. Therefore, in our analysis, we distinguish between four classes of decay of p(·; θ 0 ), and analyze each on its own.
For a diagonal L, with positive entries L 1 , . . . , L d on the diagonal (i.e., L = diag (L 1 , . . . , L d )), let us define: The families P C,L and P (r) C,L include densities that decay at a polynomial rate or faster, while E C,L includes densities that decay at exponential rate or faster, and E (2) C,L includes densities that decay at a square exponential rate or faster. Table 2 shows four well known kernels, their corresponding densities, and their decay class.
The following proposition specifics how to set U based on the decay class and the maximum value of ϕ.
Proposition 7. Suppose that, M R is such that for every η ∈ R d and every x ∈ X we have |ϕ(x, η)| ≤ M R . Then, the following establishes a U (min) = (U for every v such that v T K θ v = 1.
C,L and r > d/2: ( 2 F 1 is the hypergeometric function). Then set U (min) k .
3. If p(·; θ 0 ) ∈ E (1) Remark 8. We recommend to set U to U (min) . For P C,L , E C,L and E (2) C,L we give explicit formulas for U (min) . For P (r) C,L , it is defined implicitly as the solution to a nonlinear equation. We recommend finding the solution numerically using root-finding methods. We also give an explicit upper bound for the value of U (min) , which can be used if one wishes to avoid solving a non-linear equation, however those upper bounds tend to be loose. Nevertheless, the upper bound is used later to derive asymptotic bounds on s.
where the first inequality is due to the Cauchy-Schwartz inequality, the second inequality follows from observing that the smallest eigenvalue of K θ is bigger than or equal to σ 2 n , and the last inequality is due to the fact that every entry in z(η) has absolute value that is smaller or equal to M R .
The case of p ∈ P C,L : In this case, So, in order for Eq. (19) to hold, we set The case of p ∈ P (r) C,L : In this case, where m = min 1≤k≤d {L k U k }, and in the second inequality we use d-dimensional (d ≥ 2) spherical coordinates [7]: sin φ j where η = t, 0 ≤ φ j ≤ π/2, 0 ≤ θ < π/2, m ≤ t < ∞, and the Jacobian is Also, in the fourth equality we use the following property of the beta function: with y = 1/2 and x = (j + 1)/2, for any j = 1, . . . , d − 2, that iŝ Now, we write the last integral of (22) in terms of incomplete beta function, as follows: The expression in the last equality is the analytic continuation of the beta function B −m −2 (r − d/2, Therefore, we obtain ˆR In order for Eq. (19) to hold, let x > 0 be the solution of the equation 1 2σ 2 f n and then set U (min) k = x/L k . Note that if we replace the expression after the second equality in (22) with the upper bound then, in order for Eq. (19)to hold, we can set .
Note also that for d = 2, we have the simplest case of spherical coordinates, which yields Also, for d = 1 we do not need spherical coordinates. In that case, we have ˆR Now, U (min) that equates the last expression with 1/2σ 2 f n is obtained by solving the same equation obtained for d ≥ 3, only with d = 1.
Note that if we bound the first integral similarly to the bound in the case d ≥ 3, we obtain the same formula for U (min) only with d = 1.The case of p ∈ E (1) C,L : In this case, So, in order for Eq. (19) to hold we set The case of p ∈ E (2) C,L : In this case, where we used the bound erfc (x) ≤ e −x 2 for non-negative x, which follows from [10]. In this paper it was shown that for any non-negative x erfc (x) = 2 πˆπ 2 0 e − x 2 sin 2 θ dθ and therefore erfc (x) = 2

Approximating the Truncated Integral
The nodes η 1 , . . . , η s and their weights h 1 (θ 0 ), . . . , h s (θ 0 ) are simply rescaled multivariate Gauss-Legendre quadrature nodes and weights 3 , and so s j=1 h j (θ 0 )|z(η j ) * v| 2 is a quadrature approximation of´Q U |z(η) * v| 2 p(η; θ 0 )dη. In this section we derive a lower bound on s that guarantees that The bound on s depends on U. Together with Proposition 7, we completely specify how to build the quadrature approximation so that Eq. (17) holds.

Decay of Chebyshev Coefficients for Multivariate Functions
Our analysis relies on generalizations of existing decay bounds for Chebyshev expansions of analytic functions in one dimension to multivariate functions. In this subsection we introduce these results. Classical decay bounds for Chebyshev expansions of analytic functions in one dimension are based on bounding the function values on the Bernstein ellipse (see [20,Section 1.4] for further details). For multivariate functions, a polyellipse is used instead.
Note that E 1,β U , where β U denotes entrywise division between β and U, is a Bernstein polyellipse, and that E U,β = U E 1,β U . So, E U,β is a polyellipse with foci at ±U k . For a multivariate analytic function f on [−1, 1] d , the multivariate tensorised Chebyshev expansion is given by where the coefficients are given by Although Theorem 10 has essentially been proven in [39], for completeness we include in Appendix B our proof of Theorem 10 which is based on a different technique.

Then for
Note that in this case Proof. For conciseness, we drop θ 0 from p throughout the proof. For convenience, we use the following form of the quadrature rule The sum in the right-hand side is a quadrature approximation off v (χ)p(χ), which we analyze.
To that end, we first bound the analytic continuation off v (χ)p(χ) on E 1,β U (where β U denotes entrywise division). For every η ∈ E U,β we have (similar to the derivation of Eq. (21)): We can now apply quadrature approximation bounds onf v (χ)p(χ) to bound the error be the multivariate Chebyshev expansion off v (χ)p(χ), and let P 2s−1 be the truncated expansion: Similarly to the strategy employed in [38] and [37,Theorem 19.3], we have In the first equality, we use the following equalitŷ which follows from the exactness of the Gauss-Legendre quadrature in one dimension: Hence, .

Finally, bounding
for each k = 1, . . . , d , gives the bound from the theorem and the statement now follows immediately.
The last theorem allows us to compute the required s based on U and the singularities in p(·; θ 0 ). In Section 6, we show concrete examples for using Theorem 11 to bound s and U. The main step is finding the polyellipse parameters. If the analytic extension of p(·; θ 0 ) has its first singularity at the pure imaginary value x 0 = ±iβ for some β = (β 1 , . . . , β d ) T such that 0 < β k < U k , thenp(·; θ 0 ) has its first singularity at x 0 = ±iβ U. Thus, we can choose the ellipses parameters to be

Handling a Parameter Domain Θ
Given a hyperparameter domain Θ, we want to set the parameters U = (U 1 , . . . , U d ) and s = (s 1 , . . . , s d ) such that Eq. (17) hold for every θ ∈ Θ. This way, the parameterized family of positive definite kernel approximations {k θ } θ∈Θ given by Eq. (11) with these U and s is n-spectrally equivalent to {k θ } θ∈Θ on the data domain X . To do that, we need to find the worst-case (over θ ∈ Θ) parameters U and s. The following gives a general end-to-end statement.
Let M U,β := sup Then for the parameterized family of kernel approximations given by Eq. (11) is n-spectrally equivalent to {k θ } θ∈Θ on the data domain X . Furthermore, if we set U = U (min) and s according the last lower bound we have To use this theorem, one needs to bound the decay of the density functions p(·; θ 0 ) over Θ and calculate an upper bound on C U,β . In general this might be hard, but luckily in most kernels display a monotonicity in their hyperparameters that helps identify the worst case for θ over Θ. For example, for the one dimensional Gaussian kernel exp(− −1 (x − x ) 2 2 //2), the various parameters in the theorems monotonically increase as → 0. In the next section we given concrete examples for using Theorem 12 for the kernels listed in Table 2.

Examples of Feature Maps for Kernels
In this section we show how to apply the theory presented in the previous section to design n-spectrally equivalent kernel approximations for a few widely used kernel functions. Throughout this section, we assume that n is fixed, the data domain is X ⊆ R d , and the hyperparameter domain is Θ. Furthermore, we assume that we have a bounding box on the domain, i.e. X ⊆ d k=1 [−R i /2, R i /2] (obviously, such a bounding box can be easily computed from the input data). Let R = [R 1 , . . . , R d ].

Gaussian Kernel
Recall that the Gaussian kernel is where we added a scaling factor σ 2 f and included the ridge term σ 2 n γ(x − x ) in the kernel definition. Note that for conciseness, we consider the isotropic version; the formulas can be modified for the anisotropic case. As discussed in Section 4, by setting ϕ(x, η) = e −ix T η , this kernel matches the form of Eq. (8). We assume that the hyperparameters are bounded as follows: we have a ridge term; our method is not able to approximate the Gaussian kernel in the absence of a ridge term).
The density is given by Therefore, for θ ∈ Θ , p(·; ) ∈ E (2) We also have |ϕ(x, η)| ≤ 1 for all x and η, so we set M R = 1. So, based on Proposition 7 we set In addition, p(·; ) is analytic on R d , and in particular it is analytically continuable to the polyellipse E U,β with β k = 2U k , ρ k = 1 + √ 2, as described in Theorem 11. Now we bound p(·; 0 ) on the polyellipse as follows, for Hence, C U,β = d 0 (2π) −d/2 exp 2 0 U 2 2 /2 . Recall that z : R d → C n is given by z(η) j = e −iη T xj , it is easy to verify that for every v ∈ R n the function |z(·) * v| 2 is an analytic function on R d . In particular, it is analytically continuable to the polyellipse E U,β with β k = 2U k , ρ k = 1+ √ 2 as required by Theorem 11. Now we bound |z(η) j | for each j on its corresponding ellipse. For any x, y ∈ R d : |z(x + iy) j | = |z(iy) j |, so we need to bound |z(iy) j | for y ∈ d k=1 [−U k , U k ]: |z(iy) j | = e −i·iy T xj = e y T xj ≤ e y 2 xj 2 ≤ e U 2 R 2/2 .
Hence, |z(η) j | ≤ e U 2 R 2/2 =: M U,β . Now, one can apply Theorem 12 with the these parameters and obtain that for 2 ln(1 + √ 2) + 1 we have the desired bound. Since U 2 = O( √ ln n) (in particular, U 2 2 = O(ln n)) and assuming the bounding box R is fixed then s k = O(ln n) suffice and s = d k=1 s k = O((ln n) d ) suffices.

Matèrn Kernel
Recall that the Matèrn kernel is where we added a scaling factor σ 2 f and an included the ridge term σ 2 n γ(x − x ) in the kernel definition. Note that for conciseness, we consider the isotropic version for a fixed ν; the formulas can be modified for the anisotropic case. As in the case of the Gaussian kernel, by setting ϕ(x, η) = e −ix T η , this kernel matches the form of Eq. (8). We assume that the hyperparameters are bounded as follows: we have a ridge term; our method is not able to approximate the Matèrn kernel in the absence of a ridge term).
The density is given by . Therefore p ∈ P , as required by Theorem 11. Now we bound p(·; 0 ) on the polyellipse as follows, for x k ∈ − where the maximum value is obtained at the nearest points to the poles: y k = ± Recall that z : R d → C n is given by z(η) j = e −iη T xj . Now we bound |z(η) j | for each j on its corresponding ellipse. For any x, y ∈ R d : |z(x + iy) j | = |z(iy) j |, so we need to bound |z(iy) j | for y ∈ d k=1 [−β k /2, β k /2]: |z(iy) j | = e −i·iy T xj = e y T xj ≤ e y 2 xj 2 ≤ e β 2 R 2 /4 Hence, |z(η) j | ≤ e β 2 R 2/4 =: M U,β . Now, for the asymptotics, consider the upper bound for U in Theorem 12. If one denotes γ = ln 2 2d+2 π −d/2 Γ(ν+d/2) Γ(ν) n0 σ 2 f 0 n 2 and δ = ln , then Theorem 12 can be applied with the these parameters and obtain that (for d ≥ 3) we have the desired bound.
Assuming the bounding box R is fixed, s k = O(ln n) suffice and s = d k=1 s k = O((ln n) d ) suffices.

Semigroup Kernels
The previous two examples were of shift-invariant kernels, and the feature mapping ϕ was based on Bochner's theorem. In this section, we demonstrate the application of our theory to a different type of kernels: semigroup kernels [43]. These type of kernels require a slight modification of our setup, which we briefly describe below, but adjusting theory itself is technical and we omit it. Semigroup kernels are well-suited for non-negative data, i.e. X ⊆ R d + , and require that the kernel value at x and x depends only on the sum x + x : k(x, x ) = k 0 (x + x ). One example of such kernel is the reciprocal semigroup kernel: where we add a scaling factor σ 2 f and an included the ridge term σ 2 n γ(x − x ) in the kernel definition. It can be shown that every semigroup kernel can be written in the following integral form [6], which is analogous to Eq. (9): where p(·; λ) is a probability density function which is supported only on R d + . For the reciprocal semigroup kernel we have p(η; λ) = λe −λ η 1 = λe −λ d k=1 η k , so p(·; λ) ∈ E (1) C,L with C = λ, L = λI d , which is analytic on R d + . Thus, we see that semigroup kernels can be represented as which is almost the same as Eq. (8), except the integration area is R d The construction of Gauss-Legendre features is quite similar to the integration area is R d , except that we replace the assumption that X ⊆ for k = 1, . . . , d. We omit the details of the construction, since they mostly repeat the construction described in Section 5. Now consider the reciprocal semigroup kernel. We assume that the hyperparameters are bounded as follows: n0 > 0 (i.e., we have a ridge term). We set In addition, since p(·; λ) is analytic on R d + , and in particular it is analytically continuable to the polyellipse E U,β with β k = 2U k , ρ k = 1 + √ 2. Now we bound the analytic continuation of p(·; λ 0 ) (which we also denote by p(·; λ 0 )) on the polyellipse as follows. For any Hence, C U,β := λ 0 e λ0 √ 2−1 2 U 1 . For semigroup kernels we use z : R d + → R n defined by z(η) j = e −η T xj as the feature map. It can be seen that for every v ∈ R n the function |z(·) * v| 2 is an analytic function on R d + , and we also have |z(η) j | ≤ 1 = M R for η ≥ 0 and for all 1 ≤ j ≤ n. In particular, it is analytically continuable to the pollyellipse E U,β with β k = 2U k , ρ k = 1 + √ 2 . Now we bound |z(η) j | 2 , where here z denotes the analytic continuation. Notice that for any Since U 1 = O(ln n) and assuming the bounding box R is fixed, then s k = O(ln n) suffice and s = d k=1 s k = O((ln n) d ) suffices.

Numerical Experiments
In this section we report experiments evaluating the performance of our proposed quadrature based approach. Our goal is to show that indeed if U and s are set to be large enough, our method yields results that are essentially indistinguishable from using the exact kernel, while offering faster hyperparameter learning, training and prediction. Clearly, from the theoretical results, our method predominately applies to lowdimensional datasets (for example, such datasets are prevalent in spatial statistics), so we experiment with one dimensional and two dimensional datasets. We experiment both with the Gaussian kernel or the Matèrn kernel.
In the graphs, we label our method as GLF-GPR (standing for Gauss-Legendre Features Gaussian Process Regression). We use the following methods as benchmark: exact GPR (labeled in the graphs as Exact-GPR) and GPR based on random Fourier features (labeled RFF-GPR). As performance metric we use the MSE error on a test set (as a function of number of features) and the time to learn the hyperparameters. Training and prediction time of both GLF-GPR and RFF-GPR are essentially the same for the same number of features, and both are faster than Exact-GPR if the number of features is smaller then the training set size. Thus, when it comes to training and prediction time, it is sufficient to explore the test error as function of the number of features. However, hyperparameter learning time can vary considerably between GLF-GPR and RFF-GPR, so we compare this quantity directly.
The various methods were implemented in MATLAB. Optimizing the hyperparameters was conducted using the MATLAB function fmincon after transforming the hyperparameters to a logarithmic scale. For each problem we defined a hyperparameter domain, e.g.
f 0 } and we take the initial hyperparameters for the optimization to be [ 0 , σ 2 f 0 , σ 2 n0 ]. Running times were measured on a machine with two 3.2GHz Intel(R) Xeon(R) Gold 6134 CPUs, each having 8 cores, and 256GB RAM.

Synthetic Data
In this subsection, we report experiments on synthetically generated data. The data is generated by noisily sampling a predetermined function, i.e. samples are generated from the formula where f is the true function and {τ i } are i.i.d noise terms, distributed as normal variables with variance σ 2 τ = 0.5 2 (for 1D) or σ 2 τ = 0.3 2 I 2 (for 2D). In these experiments we use the isotropic Gaussian kernel. First, we consider a one dimensional function: The function was sampled equidistantly on [−1, 1] with n = 800 samples. The results are reported in Figure 1, where we show how GLF-GPR with the number of quadrature points s compared to RFF-GPR.  In the top-right graph, we see that the log-likelihood of GLF-GPR merges with the log-likelihood of Exact-GPR for each of the hyperparameters, where the optimal hyperparameters are dashed in blue. RFF-GPR deviates considerably. This graph exemplifies that GLF-GPR can yield good approximation to the exact log-likelihoods, while RFF-GPR yields a poor approximation. We also see that GLF-GPR optimizes hyperparameters that are much closer to the exact values than RFF-GPR. The bottom-left plot shows the MSE error on the same test points. We see that the GLF-GPR error stabilizes on error of Exact-GPR even before the theoretical value of s. The bottom-right graph shows the runtime of the hyperparameter learning phase for different values of quadrature points s. GLF-GPR is clearly more efficient than Exact-GPR and mostly more efficient than RFF-GPR. As expected, as s becomes larger, GLF-GPR learn the hyperparameters much faster than Exact-GPR and RFF-GPR. Furthermore, GLF-GPR achieves a low error rate with less features than RFF-GPR, and thus is able to do training, prediction and hyperparameter learning much faster than RFF-GPR.
Next, we consider a two dimensional function: The function was sampled on an uniform grid on [−1, 1] × [−1, 1] with n = 4096 samples. We consider 1 = 2 so U 1 = U 2 and s 1 = s 2 , i.e., s = s 2 1 . The results are reported in Figure 1. Similar to the synthetic 1D experiment, in the top-right plot we see that GLF-GPR yields a good approximation to the exact log-likelihood. In the bottom-left plot we see that shows the GLF-GPR error stabilizes on error of Exact-GPR at a much smaller number of quadrature points than the theoretical s.

Natural Sound Modeling
Next we consider the natural sound benchmark used in [42] (without hyperparameter learning) and [13] (with hyperparameter learning). The data is shown in the top-left graph of Figure 3. The goal is to recover contiguous missing regions in a waveform with n = 59309 training points. The test consists of 691 samples. The Gaussian kernel is used for learning. Intensity

Test Data and Predictions
Test Data RFF-GPR GLF-GPR Exact-GPR Results are reported in Figure 3. In the bottom-left graph we plot the test error as a function of the number of features. Initially GLF-GPR produces poor results, but when s is large enough, the results are similar to the Exact-GPR (see also the top-right plot). We see that even when the number of quadrature points is smaller than the theoretical value required for spectral equivalence, GLF-GPR's error is stabilizes on the Exact-GPR error. In contrast, RFF-GPR's error oscillates above Exact-GPR's error. In the bottom-right graph we see that the runtime of the hyperparameters learning phase is significantly smaller for GLF-GPR.

Google Daily High Stock Price
We consider a time series data of the daily high stock price of Google spanning 3797 days from 19th August 2004 to 19th September 2019. We set the data as x ∈ {1, . . . , 3797} and y = log(Stock high ). The test is of size of 12% of the data, i.e., consists of 502 days. We use the Matèrn kernel with ν = 5/2.
We note that the theoretical number of quadrature features s required for spectral equivalence is bigger than the number of training points. Possible reasons are: the hyperparameter σ 2 n in these dataset is very small and that increases our bound in Eq. (21). This increases U which increases s. In addition, the weight function of the Matèrn kernel has singularity point which leads to the ellipse parameter to be pretty small. However in practice we see that the approximation convergences around the value s = 1550.   Results are reported in Figure 4. From the top-right and top-left plots we see that RFF-GPR is producing a poor approximation while GLF-GPR approximation merges with Exact-GPR approximation. Also, from the bottom-left plot we see that initially GLF-GPR produces poor errors, but when s is large enough, the error stabilizes on the Exact-GPR error (see also the top-right plot).

Spatial Temperature Anomaly for East Africa in 2016
Similar to [34], we consider MOD11A2 Land Surface Temperature (LST) 8-day composite 2D data of synoptic yearly mean for 2016 in the East Africa region. For the training set, we randomly sample 77404 LST locations and set x ∈ {(Longitude, Latitude)} and y = {temperature}.We examine the MSE errors on the remaining 6005 locations, but use all 83409 data points to draw maps. We also use the anisotropic Matèrn kernel with ν = 1. Again, theoretical number of quadrature features for spectral equivalence is bigger than then number of training points. However again in practice we see that the approximation convergences with less features. Due to memory and time constraints, we were unable to use Exact-GPR, and RFF-GPR results are presented up to the computer's memory capacity. Results are reported in Figure 5. In the top plot we see that GLF-GPR approximates the true function well, unlike RFF-GPR. Also, from the bottom-left plot we see that around s = 21025, the GLF-GPR error stabilizes while RFF-GPR error is still suboptimal.

Conclusions and Future Work
In this paper, we proposed the use of Gauss-Legendre feature for large-scale Gaussian process regression. Our method is very much inspired by Random Fourier Features [24]. However, our method replaces Monte-Carlo integration in RFF with a Gauss-Legendre quadrature of a truncated integral representation of the kernel function. With Gauss-Legendre quadrature our method is able to build spectrally equivalent kernel approximation with an amount of features which is asymptotically poly-logarithmic in the training size. In contrast, with RFF the number of features for spectral equivalence must be at least linear. Sublinear amount of features can also be obtained using a Gaussian quadrature (suggested in the context of kernel learning in [12]). However, this is problematic in the context of hyperparameter learning (see Section 4.3). RFF has a similar issue. In contrast, the use of Gauss-Legendre quadrature allows our method to keep the quadrature nodes and weights fixed, leading to simplified structural dependence of the kernel matrix on the hyperparameters which is more amenable to hyperparameter learning. Finally, we demonstrate the utility of our method on several real-world low-dimensional datasets.
We mention a few possible directions for future research: • Asymptotically, our method requires a number of features that is poly-logarithmic in the training size. Yet, for some moderately sized datasets our theoretical results required a number of features larger than the number of training points. However, in practice the number of features required for high quality results was much smaller than the bound. Closing this gap is an open problem.
• Our method is able to handle a rich family of kernels (see Eq. (8)) which includes stationary kernels and some non-stationary kernels (see Section 6.3). Extending our method to arbitrary non stationary kernels is an open problem.
• The number of features needed by our method is exponential in the dimension, i.e. we have not escaped from the curse of dimensionality. A future research direction is to replace the tensorized multivariate quadrature with sparse grids, and in doing so avoid the exponential dependence on d.

A.2 Efficient Gaussian Process Regression using QR Decomposition
Here we present alternative formulas to the ones presented in Section 4.2 and based on QR decomposition instead of the normal equations. Such formulas are likely to be more numerically robust. Let Z = Q Z R Z be a thin QR decomposition of Z, i.e. Q Z ∈ C n×s is such that Q * Z Q Z = I n and R Z ∈ C s×s is an upper triangular matrix. We suggest to compute the QR decomposition of Z in lieu of computing Z * Z, and keeping only R Z and Q * Z y so still only O(s 2 ) is needed. There is no asymptotic penalty in terms of arithmetic operation count since the decomposition can be computed in O(ns 2 ) operations. However, Z * Z tends to be ill-conditioned due to the squaring of the condition number of Z, so it is best to avoid computing it. Note that the QR decomposition is computed only once, and not per iteration.
Let us consider a specific iteration, and for conciseness we omit θ for the following formulas. Let We compute a thin QR decomposition A = Q A R A of A (O(s 3 ) operations) and write where Q 1 , Q 2 ∈ C s×s , i.e., Hence, which is a QR decomposition of The crux is that given the QR decomposition of Z, we can compute the QR decomposition of B in O(s 3 ) arithmetic operations instead of O(ns 2 ). We now compute which implies that w can be computed in O(ns) operations (since W is diagonal). Similarly, we also have i.e., F(θ) can be computed in O(s 3 ) operations. This allows us to compute the first two formulas in Eq. (14) in O(s) time. As for the third formula, we have (25). Since ∂W(θ)/∂θ i is diagonal, we need to compute only the diagonal of Z * Z and Z * ZF(θ). The diagonal of Z * Z is just the square norms of the columns of Z, and can be precomputed in O(ns). Furthermore, in some cases we know analytically the values of this norm. For example, for shift-invariant kernels we use ϕ(x, η) = e −ix T η so the squared norms of the columns of Z is equal to n. As for Z * ZF(θ), we have: Note that Q * 1 R Z has already been computed for F(θ). Since we only need the diagonal of Z * ZF(θ), and this is an Hermitian matrix, the diagonal is just the squared norms of the columns of Q * 1 R Z . Thus, after O(s 3 ) preprocessing, for every θ i the first term in Eq. (13) can be computed in O(s) operations .
As for the first identity in Eq. (13), we still have α = σ −2 n y − σ 2 f ZWw and Z * α = w. So, the second identity can be computed as previously described using w, which is computed according to Eq. (26).

B Analysis of Multivariate Chebyshev Approximation in High Dimensions
The following theorems are generalizations of similar one-dimensional theorems. All the proofs rely on ideas similar to the ones presented in [18,19,20] and [37, Theorem 3.1, Theorem 8.1]. Note that a similar generalization can found in [36,39]. For completeness, we present our own proof, which is based on a different technique. For convenience, denote: as the polyellipse, where each E ρ k is Bernstein ellipse with foci at ±1 and the sum of of major and minor semiaxis lengths of the ellipse is ρ k . Also denote C r := z ∈ C d : |z k | = r k ∀k = 1, . . . , d as the polycircle centered at the origin, and simply denote C 1 in the case where r 1 = . . . = r d = 1. Finally, denote A r,R := z ∈ C d : 0 < r k < |z k | < R k ∀k = 1, . . . , d as the polyannulus centered at the origin. The following proposition appears in [27] as Theorem 1.5.26. See also [8,Pages 32,[90][91] for further details.
Proposition 13. Let A r,R be the polyannulus centered at the origin, and let f be an analytic complex function on A r,R . Also, let r k < s k < R k , k = 1, . . . , d. Then f has a multivariate Laurent expansion