extfieldsofstudy:
Mathematics
Computer Science

text:
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


Distilling Heterogeneity: From Explanations of Heterogeneous Treatment Effect Models to Interpretable Policies

Internet companies are increasingly using machine learning models to create personalized policies which assign, for each individual, the best predicted treatment for that individual. They are frequently derived from black-box heterogeneous treatment effect (HTE) models that predict individual-level treatment effects. In this paper, we focus on (1) learning explanations for HTE models; (2) learning interpretable policies that prescribe treatment assignments. We also propose guidance trees, an approach to ensemble multiple interpretable policies without the loss of interpretability. These rule-based interpretable policies are easy to deploy and avoid the need to maintain a HTE model in a production environment.

Introduction
Unlike conventional A/B testing where individuals are randomly assigned to a treatment condition, internet companies are increasingly using heterogeneous treatment effects models (HTE) and personalized policies to decide, for each individual, the best predicted treatment (Garcin et al., 2013;Toguç and Kuzu, 2020). However, to flexibly model heterogeneity in treatment and outcome surfaces, state-of-the-art HTE models used for personalization Shi et al., 2019;Kennedy, 2020) frequently leverage black-box, uninterpretable base learners such as gradient boosted trees and neural networks. To generate a HTE prediction, some of them combine outputs from multiple models via stacking, weighting, etc. (Künzel et al., 2019;Nie and Wager, 2021;Montoya et al., 2021). Moreover, their treatment effect predictions are sometimes used as inputs to additional black-box policy learning models (Imai and Strauss, 2011;McFowland et al., 2020).
There are several challenges with these approaches. Firstly, when there are multiple (not just binary) outcomes, the creation of a policy is no longer as trivial as assigning the treatment group that maximizes a single treatment effect. Secondly, as dataset size and the number of features used for HTE modeling or policy learning increases, the cost of maintaining such models in deployment increases (Sculley et al., 2015a). Finally, the black-box nature of HTE models and resulting policies makes them difficult to interpret, which can deter their uptake for critical applications. This paper describes methods for: (1) learning explanations of HTE models, to better understand such models; (2) generating interpretable policies that can perform on par with or even better than black-box policies. Our explanations and interpretable policies are constructed from if-else rules, using individual features as inputs (Figure 1). We demonstrate a variety of different ways to learn interpretable policies, from methods that act on HTE model predictions to methods that work on conditional average treatment effects. We also propose ways to ensemble multiple interpretable policies while still remaining interpretable.
Crucially, all of our proposed methods work in the general setting of more than two treatment groups and multiple outcomes, which is needed in production environments at internet companies. For example, a new website feature may increase user engagement at the expense of revenue. We do so by leveraging techniques such as model distillation (Bucilua et al., 2006;Hinton et al., 2014) with multitask (Caruana, 1997) and multiclass student models, and linear scalarization of multiple outcomes from the multi-objective optimization literature (Gunantara (2018)).

Related Work
We do not attempt to comprehensively survey the state of HTE estimation since our contribution is agnostic to the choice of HTE model. Most modern HTE estimation methods are compatible with our proposed methods.
Subgroup finding. Our work on interpreting HTE models by surfacing segments with heterogeneity produces a representation analogous to subgroup finding methods. These methods are divided into: (1) statistical tests to confirm or deny a set of pre-defined hypotheses about subgroups (Assmann et al., 2000;Song and Chi, 2007); (2) methods that aim to discover subgroups directly from data (Imai and Ratkovic, 2013;Chen et al., 2015;Wang and Rudin, 2021;Nagpal et al., 2020); (3) methods that act on inferred individual-level treatment effects to discover subgroups Foster et al. (2011); Lee et al. (2020). However, the comparison of subgroup finding algorithms is still an open question. Loh et al. empirically compared 13 different subgroup-finding algorithms (Loh et al., 2019), finding that no one algorithm satisfied all the desired properties of such algorithms.
Policy learning. The problem of policy learning from observational/experimental data has received significant attention due to its wide applicability (Beygelzimer and Langford, 2009;Qian and Murphy, 2011;Zhao et al., 2012;Swaminathan and Joachims, 2015;Kitagawa and Tetenov, 2018;Kallus and Zhou, 2018), but their interpretability and scalability has received less attention. Several existing works (Amram et al., 2020;Athey and Wager, 2021;Jo et al., 2021) do construct interpretable tree policies based on counterfactual estimation. However, they focus on constructing exactly optimal trees, which is prohibitively slow in a large scale and unsuitable for a production environment.
Interpretability via distillation and imitation learning. Black-box classification and regression models have been distilled into interpretable models such as trees, rule lists, etc. (Craven and Shavlik, 1995;Frosst and Hinton, 2017;Tan et al., 2018b,a). Craven and Shavlik (1995) distilled a neural net to a decision tree and evaluated the explanation in terms of fidelity, accuracy, and complexity. Namkoong et al. (2020) distilled a Thompson Sampling policy to policies represented using MLPs. Interpreting reinforcement learning methods has also garnered interest; see Alharin et al. (2020) for a survey. These methods can be roughly divided into those that leverage the interpretable representation provided by trees to approximate states and actions (Bastani et al., 2018;Vasic et al., 2019), saliency and attention methods (Wang et al., 2016;Greydanus et al., 2018;Qureshi et al., 2017), anddistillation (Czarnecki et al., 2019). We also use distillation techniques, but in a non-sequential setting.

Problem Setting and Notation
Suppose we have a dataset D of N points and P features with K treatment groups, 1 control group, and J outcomes. Let W i ∈ {0, 1, . . . , K} denote the treatment group indicator for the i-th individual (W i = 0 means the individual is in the control group). Let G k denote the set of individuals in the kth treatment group, k = 0, 1, . . . , K. Following the potential outcome framework (Neyman, 1923;Rubin, 1974), we assume for each individual i and outcome j there are K + 1 potential outcomes Y , the observed potential outcome (out of K + 1 potential outcomes). We suppose we have access to predicted outcomes (Ŷ is the predicted j-th outcome for individual i assuming treatment group W . These could be obtained directly from the HTE model or together with estimates of the baseline surface. We focus on deterministic policies, Π: R P → {0, 1, .., K}, which map feature vectors to treatment groups. An individual with feature vector X i will be given treatment Π(X i ). A segment, S, is a set of indices, representing a set of points, N S , treatment assignments, W S along with their predicted outcomes,Ŷ (W ) ij for i ∈ S. Throughout this work, we refer to lists of segments, {S 1 , S 2 , ..., S n } which form a partition of the input indices, i.e. there is no index overlap between the segments.

Preliminaries: Combining Multiple Outcomes
Throughout this work, we combine multiple outcomes into a single outcome via linear scalarizations, parameterized by a set of weights (c 1 , . . . , c J ). Thus, for each individual i, the multiple observed outcomes becomes a single outcome: We also combine predicted outcomes in the same way:Ŷ iJ .. Our proposed methods take (c 1 , . . . , c J ) as inputs and assume a value 1 (equal weighting on outcomes) if not provided. In practice, it is common to either tune the weights (c 1 , . . . , c M ) (Letham and Bakshy, 2019) or handpick weights to manage particular preferences or tradeoffs.

Explaining Heterogeneous Treatment Effect Models
Multi-task decision trees (MTDT) extend single-task decision trees by combining the prediction loss on each outcome, across multiple outcomes, into a single scalar loss function. This simple representation is effective and suitable for locating segments where individuals present heterogeneity across multiple outcomes. In this tree, each task (label) is the predicted treatment effect for an outcome, and each node is a segment identified to have elevated or depressed treatment effects across multiple outcomes. We learn these trees using distillation techniques. The resulting method is called Distilling HTE model predictions (Distill-HTE). Figure 2 provides a concrete example.
Distillation. To learn an MTDT model in a HTE model-agnostic fashion, we leverage model distillation techniques, taking the HTE model as the teacher and the MTDT model as a student. After obtaining predicted outcomesŶ is the predicted j-th outcome for individual i assuming treatment group W , we construct pairwise predicted treatment effectsT ij , i.e. individual i's predicted treatment effect for outcome j for a particular treatment contrast (analogous definitions will hold for all treatment contrasts). LetF be the prediction function of an MTDT model using the following distillation training objective that minimizes the mean square error between the HTE model predicted treatment effects andF: (1) F j is the MTDT model's prediction for just outcome j, accompanied by c j , a weight that encodes how much outcome j contributes to the overall loss, as reviewed in Section 3.2. When c j = 1, as we will assume by default, each outcome contributes equally to the overall distillation loss.
Figure 2: Distill-HTE: Multitask treatment effect decision tree learned by distilling a black-box HTE model. After the tree is learned, we "unroll" the tree such that each terminal node forms a segment visualized in a barplot.
Counterfactual reasoning. To improve the robustness of the subgroups, we introduce three elements based on counterfactual reasoning ideas: (1) Terminal nodes without sufficient overlap (i.e. enough treatment AND control points) are post-pruned from the tree, and the prediction regenerated; (2) Confidence intervals are provided for treatment effects within each node; (3) Honesty (Athey and Imbens, 2016): splits and predictions are determined on different data splits.
From treatment effects to policies. With K treatment groups, 1 control group, and M outcomes, the proposed method learns K MTDT models, where each model predicts treatment effects of M outcomes for each treatment group compared to control. One MTDT model can then be used to generate a policy that indicates, for each segment, which of the pair of treatment groups should be assigned to maximize treatment effects. However, generating a single policy still requires combining multiple MTDT, and it is not immediately obvious the best way to do so. This motivates our proposal of directly learning a single policy that applies to multiple treatment groups and outcomes.

Interpretable Policy Learning
To generate a single interpretable policy Π that applies to multiple treatment groups and outcomes, we now describe different approaches, all using if-else rules representations.

APPROACHES BASED ON HTE MODEL PREDICTIONS
An already trained HTE model from which predicted outcomes can be obtained is a prerequisite for the approaches described in this subsection. We describe two ways of leveraging these predictions: GreedyTreeSearch-HTE and Distill-Policy.
Greedy tree search from HTE model predictions (GreedyTreeSearch-HTE). In this approach, we directly solve the following optimization problem in the space of trees with pre-determined maximum tree depth, assuming without loss of generality that higher outcomes are better: This can be seen as a cost-sensitive classification problem, where the cost of assigning a point to a group is the negative value of the predicted potential outcome (Elkan, 2001). To achieve a scalable implementation, we solve the optimization problem greedily instead of resorting to exact tree search over the whole space of possible trees. DefineM j (S) = u i ∈SŶ (j) i andM (S) = max jMj (S). The detailed implementation can be found in Algorithm 1.

Algorithm 1 Greedy Tree Search from HTE Model Predictions (GreedyTreeSearch-HTE)
Output: A list of segments: segments. 1. Set depth := 0 and segments = {D}. 2. Add segments to segments using greedy search: while depth ≤ m do Set nodes := {} for S ∈ segments do Split S into S * l and S * r such that: If such S * l and S * r exist, we add them to nodes; otherwise we add S to nodes. end Let segments := nodes and depth := depth + 1. end This approach is similar to  and Amram et al. (2020) which respectively solve an exact tree search problem over predicted outcomes and find the optimal tree using coordinate ascent. Due to scalability concerns and the need for faster test-time inference, we take a greedy approach. On a sample of 100k data points with 12 features and three treatment groups, the Zhou et al. (2018) method took 2.5 hours to run. With a quadratic runtime in the number of data points, and the amount of data at large internet companies sometimes exceeding 100 million data points, it is impossible to run it at a scale suitable for internet companies.
Distilling policy derived from HTE model predictions (Distill-Policy). In this approach, we start from the naive policy Π implied by the outcome predictionsŶ . Then, we train a decision tree classifier on the dataset (X i , Π HT E (X i )) N i=1 , and output the segments from the decision tree classifier. The resulting policy assigns all individuals in the same segment to the majority treatment group of individuals in that segment.

APPROACHES NOT BASED ON HTE MODEL PREDICTIONS
The non-HTE-model based methods proposed in this section are useful when we are not able to train an accurate HTE model. Here we aim to generate an interpretable policy without training a HTE model, by leveraging average treatment effects. Concretely, we define the segment average treatment effect for segment S and treatment j as: where G j denote the set of individuals in the j-th treatment group. We then define a splitting metric A(S) = max j=0,1,...,K A j (S), which considers only the most responsive treatment in that segment. We build a tree and consider binary splits that improve this splitting metric. Algorithm 3 and Algorithm 4 describe two different implementations for this. In the greedy implementation, a split is only considered if both the left and right child segment have an improvement over the parent segment. In the iterative implementation, we run several iterations of splitting. For each iteration we keep the best segment, excluding it in the next run. A split is considered as long as one of the child segments improves the outcome compared to the parent node. One pitfall with this iterative implementation is that segments are not necessarily disjoint in feature space, so one individual could appear in several segments. We resolve this by always assigning the individual to the first found segment.
For both methods, if no segments are found because no eligible split exists, the policy defaults to assigning all individuals to the treatment group with highest average treatment effect. The resulting policy is: Π * = argmax j A j (S) if S ∈ segments and individual i∈ S and if i ∈ S for ∀S ∈ segments, we assign the individual to the control or default treatment. These non-model based methods can be viewed as attempts to segment the feature space using empirical conditional average treatment effects; similar efforts have appeared in the literature (e.g. Dwivedi et al. (2020)).

APPROACHES BASED ON ENSEMBLING DIFFERENT INTERPRETABLE POLICIES
We can already generate interpretable tree-based policies using the methods described above. However, different policies may exhibit different strengths in different regions, and simply training these trees with deeper depth does not necessarily improve the resulting policy in our experiments. We leverage ensemble learning to identify such regions, with the hope of generating a better policy. Concretely, suppose we have access to Q policies Π 1 , Π 2 , . . . , Π Q . We want to train an ensemble policy Π that uses all or a subset of the policies Π 1 , Π 2 , . . . , Π Q while still remaining interpretable. We introduce two ways of doing so. For ease of notation, we assume the trained policies Π 1 , ..., Π Q were obtained from another split of the dataset and we can safely use D as the validation set on which we learn the ensemble policy.
Ensemble based on Uniform Exploration (GUIDE-UniformExplore): This method is inspired by the explore-exploit paradigm in the contextual bandits literature (see Slivkins (2021) for a review). To perform this offline, we use HTE outcome predictions when the observed outcome is not revealed in the dataset. Algorithm 2 provides the implementation. The ensemble policyπ is generated using Algorithm 1 but with the HTE predictions (Ŷ Ensemble using Offline Evaluation (GUIDE-OPE): We aim to find one feature split such that the left and right children uses a different candidate policy. We do so using offline policy evaluation (OPE; see Algorithm 5). Suppose we split dataset D into D l and D r , letΠ be the ensemble policy that applies Π k to D l and Π q to D r . To createΠ, we use exhaustive search to find the optimal feature Figure 3: Example guidance tree learned using GUIDE-UniformExplore or GUIDE-OPE methods, to ensemble multiple policies while still remaining interpretable.
Algorithm 2 Ensemble based on Uniform Exploration (GUIDE-UniformExplore) Output: A list of segments: segments. 1. Generate a dataset with randomly selected policies on individuals for i = 1 to N do Randomly select a policy from Π 1 , .., Π Q . Let A i be the index of this policy, A i ∈ 1, ..., Q.
Apply this policy to individual i, assume the policy assigns treatment group k ∈ {0, 1, . . . , K}. split and candidate policies (D * l , D * r , Π * k , Π * q ) that solves max D l ,Dr max 1≤k =q≤Q OPE(D,Π). We assign to individual i policy Π * k if individual i ∈ D * l and Π * q otherwise. Both methods return trees that we call guidance trees; see Figure 3 for an example. While there exists many other ways to ensemble policies, such as SuperLearner (Montoya et al., 2021), we do not consider these methods as they result in non-interpretable policies.

Comparing Explanations of HTE Models
We compare the Distill-HTE method proposed in Section 3.3 against several other subgroup finding methods that take the predictions of HTE models as input: (1) Virtual Twins (VT) (Foster et al., 2011); (2) R2P (Lee et al., 2020). We deliberately restrict to such post-hoc methods, due to shared motivation of explaining already-trained HTE models. We also compare to a black-box model: a T-Learner that does not find subgroups but rather provides one prediction per individual.
Setup: The first dataset, Synthetic COVID, was proposed by Lee et al. (2020) and uses patient features as in an initial clinical trial for the Remdesivir drug, but generates synthetic outcomes where the drug reduces the time to improvement for patients with a shorter period of time between symptom onset to starting the trial. The second dataset, Synthetic A, was proposed in Athey and  Imbens (2016). For each dataset we generate ten train-test splits, on which we compute the mean and standard deviation of estimates. For R2P, we use the implementation of R2P provided by the authors; For VT we used our own implementation.
Evaluation: On synthetic data where ground truth treatment effects are available, we report the Precision in Estimation of Heterogeneous Effect (PEHE) (Hill, 2011), defined as P EHE = Unlike bias and RMSE metrics computed relative to ground truth treatment effects, PEHE requires accurate estimation of both counterfactual and factual outcomes (Johansson et al., 2016). We also compute between-and within-subgroup variance.
Results: Table 1 presents the results. We make a few observations: (1) Black-box vs. whitebox: As expected, GBDT T-Learners perform well as they do not have interpretability constraints, unlike all the other methods (VT, Distill-HTE, R2P), all of which modify standard decision trees while still remaining visualize-able as a tree. Yet, Distill-HTE tends to be far more accurate than other decision tree approaches, in terms of PEHE. (2) Optimization criterion: R2P, the only method of those presented here, that considers not only homogeneity within subgroups but also heterogeneity between subgroups, has the lowest between-subgroup variance. Other methods that do not try to reduce heterogeneity between subgroups do not fare so well on this metric. On the other hand, they fare better than R2P at minimizing within subgroup variance, because they do not consider the tradeoffs between minimizing within and between subgroup variance, unlike R2P. However, R2P does this at the expense of PEHE.
(3) HTE model class: The choice of HTE model matters, with Virtual Twins not performing as well when using an RF HTE model compared to GBDT HTE model. Similarly, GBDT T-Learners perform better than DT T-Learners. (4) The impact of distillation: While the T-Learner DT model did not perform well, being worst in terms of PEHE on all datasets, Virtual Twins RF, Virtual Twins GBDT and Distill-HTE that train modified decision trees have a marked improvement over T-Learner DT, suggesting that distilling a complex GBDT or RF teacher rather than learning a tree directly is beneficial, which agrees with the existing distillation literature. 33976.605 ± 705.053 * denotes no segments were found, and the resulting policy assigned all individuals to one treatment group. Table 2: Regret (lower is better) of policy-generation methods on synthetic and semi-synthetic datasets. Methods colored in red do not personalize. Methods colored in black are blackbox policies. Methods colored in blue are interpretable policies. Methods colored in green are ensembles that still remain interpretable. Best method's regret in bold. If the best policy is no personalization (assigning all individuals to one treatment group), the next best policy is also bolded. If the best policy is an ensemble policy, the constituent policies are also bolded.

Comparing Policy Learning Methods
We compare the methods proposed in Section 3.4 to (1) a black-box policy: training a T-Learner HTE model, then assigning each individual to the treatment group with best predicted treatment effects; (2) a policy which chooses the treatment for each unit uniformly at random.
Setup: Besides the synthetic datasets described in 4.1, we use other publicly-available datasets. The IHDP dataset (Hill, 2011) studied the impact of specialized home visits on infant cognition using mother and child features. A multiple-outcome dataset, Email Marketing (Hillstrom, 2008) has 64k points, and visits, conversions, and money spent outcomes (details on how we constructed potential outcomes in Appendix). We combine multiple outcomes as explained in Section 3.2. The ensemble policies (GUIDE-UniformExplore, GUIDE-OPE) are based on GreedyTreeSearch-HTE and Distill-Policy -selected because of their individual performance.
Evaluation: To evaluate the different policies Π, we define a notion of regret: where Π(X i ) is the treatment group prescribed by that particular policy for individual i and Y i is individual i's potential outcome under treatment k.
Results: Table 2 presents the results. We make a few observations: (1) In many cases, at least one interpretable policy (GreedyTreeSearch-HTE, Distill-Policy, No-HTE, ensembled policies) improved over or is not far behind the black-box policy.
(2) Personalization: In datasets like IHDP where the benefit from personalization was limited, as the majority of points have positive treatment effects, No-HTE-Greedy and No-HTE-Iterative were able to pick up on this best policy and assign all points to one treatment. (3) The impact of ensembling: In general, ensembled policies GUIDE-UniformExplore and GUIDE-OPE improved over their constituent policies GreedyTreeSearch-HTE and Distill-Policy. However, ensemble methods with deeper guidance trees are not necessarily better. Hence, the tradeoff between one more level of guidance, and reduced interpretability, should be considered.

Ensembling Interpretable Policies
We now display several of the interpretable policies, and show how they can be ensembled while remaining interpretable. On Synthetic A, a simple dataset with only two features, the following policies were obtained: GreedyTreeSearch-HTE: If feature0 < -0.019, assign control. Otherwise, assign treatment. Distill-Policy: If feature0 <= -0.02, assign control. Otherwise, assign treatment. We applied the GUIDE-UniformExplore and GUIDE-OPE approaches to learn guidance trees: trees that ensemble these two policies while still remaining interpretable. While the GUIDE-OPE guidance tree suggested simply following Distill-Policy, the GUIDE-UniformExplore guidance tree suggested: GUIDE-UniformExplore: If feature1 < -0.426, follow Distill-Policy. Otherwise, follow GreedyTreeSearch-HTE. By leveraging feature1 which had so far not been used in the individual policies, GUIDE-UniformExplore returned the lowest regret of all methods, including a regret lower than that of the individual policies it ensembled, GreedyTreeSearch-HTE and Distill-Policy (Table 2). Interestingly, even if the individual policies GreedyTreeSearch-HTE and Distill-Policy were grown to greater depth, feature1 was not still not selected. While we do not always see gains from ensembling, in some datasets guidance trees can correct for the greediness of tree learning. While typical ways of ensembling trees (e.g. random forest) reduces interpretability, depth-constrained guidance trees, added to the top of an interpretable policy tree, merely makes the policy tree a bit deeper.

Bridging Explanations and Policies
We now present the explanations and resulting interpretable policies on the Synthetic COVID data. Interpreting HTE model: Figure 4.4 presents segments found by Distill-HTE on the synthetic COVID data. The segment with the most negative predicted treatment effect (first segment; red color), at -0.0304 +-0.0010, covers individuals who started taking the drug between 4.5 and 10.5 days after the onset of symptoms and had Aspartete Aminotransferase and Lactate Dehydrogenase levels within normal ranges (MedicineNet, 2021), suggesting that they were not extremely sick. It is therefore not surprising that they were not predicted to benefit as much from the drug.
On the other hand, the segment with the most positive predicted treatment effect (last segment; green color) covers individuals who started taking the drug soon (<= 4.5 days) after exhibiting COVID symptoms. These individuals were predicted to benefit the most from the drug, with treatment effect 0.0187 +-0.0006. This agrees with the finding that the Remdesivir drug results in a faster time to clinical improvement for the patients with shorter time from symptom onset to starting trial (Wang et al., 2020).
Policy Generation: A simple policy can be derived from the Distill-HTE segments, with all red segments assigned to not receiving the drug (as they are predicted to not benefit from it) and all green segments assigned to receiving the drug. This policy is different from the interpretable policies learned directly (Table 2). For example, the GreedyTreeSearch-HTE policy is extremely simple: if days from symptoms onset to starting trial < 11, assign to receive the drug. Otherwise, assign to not receive the drug. Another example is No-HTE-Greedy and No-HTE-Iterative whose rather stringent splitting criteria (Equation 3) did not find enough heterogeneity worth assigning different treatment groups to, and assigned all individuals to receive the drug.

Discussion and Conclusion
The motivation for this work was three-fold. (1) HTE models sometimes overfit on extremely large datasets, due to having a large number of data points, extremely noisy outcomes, noisy features, etc.
(2) Interpretable policies can avoid much of the tech-debt ML-based policies are known to incur (Sculley et al., 2015b). (3) Scalability constraints excluded many optimal tree-based policy learning algorithms from deployment in production environments, as illustrated in Section 3.4.
The interpretable policy methods proposed in this paper are not theoretically optimal, a shortcoming we acknowledge. Nonetheless, we believe that the comparison to sensible baselines, such as policies derived from T-Learner HTE models, and the ability of some policies to achieve low regret on simple synthetic datasets, demonstrates their merit.
In real, large-scale datasets where the potential benefit from personalization was higher, No-HTE-Greedy and No-HTE-Iterative suggested more personalized policies, especially when feature selection was performed beforehand. This is because the splitting criteria 3 included all features. If feature selection is not performed, this, together with using sample average treatment effects as our splitting criteria makes it easier to find a good split on the training dataset by chance but not on the test set. Further work is needed to realize the full potential of these methods.
In practice, we have seen that the choice between deriving policies from black-box HTE models and directly learning ground-up interpretable policies depends on the ability of the individuals utilizing such policies to maintain HTE models in production, and the need to explain how exactly personalization is happening, especially for applications with where personalization based on sensitive features incurs disparate treatment unfairness (Lal et al., 2020). While HTE models are resource and maintenance intensive, the ability to continuously retrain the model allows for adjustment to a dynamic user population. Conversely, interpretable policies are easy to implement and maintain, but may not perform best over the long-term without policy regeneration as the features of the user population shift.
Algorithm 4 No-HTE-model policy generation: iterative implementation (No-HTE-Iterative) Input : Maximal tree depth d, number of iterations t, dataset D. Output: segments 1. Set iterations := 0 and segments = {}. 2. Add to segments iteratively while iterations ≤ t do Step 1 to 2 in Algorithm 3, but change all min to max. Denote the resulting list of segments as C.
Add to segments the S in C that maximizes A(S) and remove the data points in S from D. Set iterations := iterations + 1. end