Structural Kernel Search via Bayesian Optimization and Symbolical Optimal Transport

Despite recent advances in automated machine learning, model selection is still a complex and computationally intensive process. For Gaussian processes (GPs), selecting the kernel is a crucial task, often done manually by the expert. Additionally, evaluating the model selection criteria for Gaussian processes typically scales cubically in the sample size, rendering kernel search particularly computationally expensive. We propose a novel, efficient search method through a general, structured kernel space. Previous methods solved this task via Bayesian optimization and relied on measuring the distance between GP's directly in function space to construct a kernel-kernel. We present an alternative approach by defining a kernel-kernel over the symbolic representation of the statistical hypothesis that is associated with a kernel. We empirically show that this leads to a computationally more efficient way of searching through a discrete kernel space.


Introduction
In many real-work applications of machine learning, tuning the hyperparameters or selecting the machine learning method itself is a crucial part of the workflow. It is often done by experts and data-scientist. However, the number of possible methods and models is constantly growing, and it is becoming increasingly important to automatically select the right model for the task at hand. Bayesian optimization (BO) is a prominent method that can be used for model selection and hyperparameter tuning. It can handle black-box oracles with expensive function evaluations, which are two characteristics often encountered when doing model selection [17]. Important applications in this context are choosing the hyperparameters in the training process of deep neural networks (DNN) [17], or dealing with discrete and structured problems like choosing the architecture of DNN's [5,11].
Gaussian processes (GP) are another important model-class. They are often utilized as surrogate models in BO [17], for time-series and statistical modeling [7,3] or in active learning loops [23,15]. The properties of GP's are mainly governed by its kernel that specifies the assumptions made on the underlying function. Choosing the right kernel is therefore a crucial part of applying GP's and is often done by the expert. Recent work [9] treated the kernel selection as a black-box optimization problem and used Bayesian optimization to solve it. This allowed searching over a highly structured, discrete space of kernels. However, their proposed kernel-kernel measures the distance between two GP's directly in function space, which is a computationally expensive task itself. This makes the method difficult to apply for the frequent scenarios where the evaluation of the model selection criteria requires only a medium amount of time.
We propose measuring the distance between two kernels via their symbolical representation of their associated statistical hypothesis. We utilize the highly general kernel-grammar, presented in [3], as underlying kernel space, where each kernel is build from base kernels and operators, like e.g.
forming effectively a description of the statistical hypothesis that is modeled by the GP. Our main idea is to build a distance over these symbolical descriptions, rather than measuring the distance between two GP's directly in function space. We employ optimal transport principles, known from neuralarchitecture search (see [11], [5]), to build a pseudo-distance between two hypotheses descriptions and use it to construct a kernel-kernel, which is subsequently utilized in the BO loop. We will show that the induced kernel search method is more efficient, in terms of number of function evaluations and computational time, compared to alternative kernel search methods over discrete search spaces.
The main challenge we encountered is the quantification of dissimilarity between two symbolical representations of kernels. We use the tree representation of each symbolical description and apply optimal transport distances over tree features. Subsequently, we empirically show that the deduced GP over GP's provides a well-behaved meta-model and show its advantages for kernel search. In summary, our contributions are: 1. We construct a pseudo-distance over GP's that acts over the symbolical representations of the underlying statistical hypothesis.
2. We use the pseudo-distance to construct a novel "kernel-kernel" and build it in a BO loop to do model selection for GP's. 3. We empirically show that our meta-GP model is well-behaved and that we outperform previous methods in kernel search over discrete kernel spaces.
Related Work: There is considerable existing work on constructing flexible kernels and learning their hyperparameters [1,20,21,18,22]. In these kinds of works, the structural form of the kernel is predefined and the free parameters of the kernel are optimized, often via marginal likelihood maximization. While being able to efficiently finding the hyperparameter due to the differentiability of the marginal likelihood, one still needs to predefine the structural form of the kernel in the first place. This is a hard task as one need to decide if long-range correlation or nonstationarity should be considered or if dimensions are ignored or not. Our method is build to automatically select the structure of the kernel. Some of the mentioned methods [20,21] are able to approximate any stationary kernel via Bochner's Theorem and therefore consider a broad kernel space themselves. However, even elementary statistical hypotheses require nonstationary kernels, such as linear trends modeled by the linear kernel. Our search space is not restricted to stationary kernels.
We consider a highly general, discrete kernel space that is induced by the kernel grammar [3]. Recent work [9] used BO to search through this space via a kernel that measures similarity in function space. We also utilize BO but employ a fundamentally different principle of measuring the distance, which is computationally more efficient. We dedicate Section 3.2 to a more precise comparison to the method of [9]. Additionally, the original kernel grammar paper [3] suggested greedy search for searching through the kernel grammar. Furthermore, [4] employed a genetic algorithm based on cross-over mutations. We empirically compare against all approaches.
In the area of BO over structured spaces, our method is most similar to the neural-architecture search (NAS) procedures presented in [11,5] who use optimal transport distances over features extracted from the graph-representation of neural networks. Compared to these methods we use OT principles to do model selection for GP's which is a fundamentally different task.

Background and Set-up
Our main task is efficient model selection for Gaussian processes. In order to provide background information, we give a small introduction to GP's and model selection for GP's. Subsequently, we present a review of the kernel-grammar [3] and show how we use it for our approach. We show how a kernel can be represented via a symbolical description. In Section 3, we will present how we use the symbolical description to construct a kernel over kernels and how we use it in the BO loop.

Gaussian Processes
A Gaussian process is a distribution over functions f : X → R over a given input space X which is fully characterized via the covariance/kernel function k(x, x ) = Cov(f (x), f (x )) and the mean function µ(x) := E[f (x)]. We therefore can write as shorthand notation f ∼ GP(µ(·), k(·, ·)). The kernel can be interpreted as a similarity measure between two elements of the input space, with the GP assigning higher correlations to function values whose inputs are more similar according to the kernel. Furthermore, the kernel governs the main assumptions on the modeled function such as smoothness, periodicity or long-range correlations and therefore provides the inductive-bias of the GP. An important property of GP's is that they are not restricted to euclidean input spaces X ⊂ R d , but can also be defined on highly structured spaces like trees and graphs, a property we will later use to define a GP over GP's.
While our method might also be used for model selection in classification, we consider from now on Gaussian processes regression. For regression, a dataset D = (X, y) with X = {x 1 , . . . , x N } ⊂ X and y = (y 1 , . . . , y N ) ∈ R N is given, where we suppose that f ∼ GP(µ(·), k(·, ·)) and . Given the observed data D the posterior distribution f |D is again a GP with mean and covariance functions [13]). Probabilistic predictions can be done via the resulting predictive distribution p(f * |x * , D) = N (µ D (x * ), k D (x * , x * )).

Model selection for GP's
Typically, the kernel k θ comes with a set of parameters θ that can be learned via maximization of the marginal likelihood p(y|X, θ, σ 2 ) = N (y; µ(X), k θ (X, X) + σ 2 I) or via maximum a posteriori (MAP) estimation, in case the parameters θ are equipped with a prior p(θ). This procedure is sometimes also called model selection, as one selects the hyperparameters of the kernel given the data (see [13]). However, we consider selecting the structural form of the kernel itself. The structural form of the kernel determines the statistical hypothesis that is assumed to be true for the data-generating process. Intuitively, the kernel is similar to the architecture in deep neural networks, which induces an inductive bias.
Our goal is to do model selection over a discrete, possibly infinite space of kernels K := {k 1 , k 2 , . . . }. As each kernel comes with its own parameters, we are actually dealing with a space of kernel families. Thus, when mentioning a kernel k we associate it with its whole family over parameters {k θ |θ ∈ Θ}. Once a kernel is selected, predictions are done with learned kernel parameters (that usually are a by-product of calculating the model selection criteria). The parameters θ are potentially equipped with a prior p(θ) depending on the selection criteria. As mean function, we always utilize the zero mean function µ(x) := 0 in case y is centered and a constant mean function otherwise. This is a common choice in GP regression. Given some model selection criteria g : K → R our task is solving k * = arg max k∈K g(k|D). While our method is not restricted to a specific model selection criteria, we focus on the model evidence p(y|X, k), which is a well-known selection criteria for probabilistic models (see [9,14,8]). Given a prior on the kernel parameters p(θ) and the likelihood variance p(σ 2 ), the log-model evidence of the marginalized GP is given as This quantity can be approximated via Laplace approximation of p(θ, σ 2 |D) (see Appendix A for details). Computing this approximation includes performing a MAP estimation of the GP parameters. Thus, once the log evidence has been computed, learned kernel hyperparameters θ MAP are automatically provided. Performing the MAP estimation scales cubically in the data set size N , which renders model selection for GP's computationally intense.

Kernel Grammar
The kernel grammar introduced by [3] specifies a highly general search space over kernels. The grammar is based on the observation that kernels are closed under addition and multiplication, meaning, for kernels k 1 (x, x ) and k 2 (x, x ) also k 1 (x, x ) + k 2 (x, x ) and k 1 (x, x ) × k 2 (x, x ) are kernels. Given some base kernels, such as the squared exponential kernel SE, the linear kernel LIN, the periodic kernel PER or the rational quadratic kernel RQ, different statistical hypotheses can be generated via addition and multiplication. For example, LIN + PER × SE describes a linear trend with a locally periodic component, such that it might be a useful hypothesis for time-series applications. For multidimensional data the base kernels are applied to single dimensions, denoted e.g. with SE i for the squared exponential kernel defined on dimension i. The grammar contains many popular hypotheses such as the ARD-RBF kernel that can be expressed via multiplication over the dimensions d i=1 SE i as well as additive models d i=1 SE i or polynomials of order m with m j=1 LIN. The specification of the complete search space is given as a grammar, where a base kernel is denoted as B and a subexpression is denoted with S. For example, in the expression LIN + (PER × SE) the expression PER × SE is a subexpression. Starting with all base kernels, the search space is defined as all kernels that can be reached via the following operations: 1. Add a base kernel to a subexpression: S → S + B 2. Multiply a subexpression with a base kernel: S → S × B 3. Exchange a base kernel with another base kernel: B → B Starting from the set of base kernels, these operations can lead to all algebraic expressions, showing the expressiveness of the grammar [3]. The authors of [3] suggest using greedy search to search through the grammar, where the kernel with the highest value for the selection criteria is selected and expanded with all possible operations. After the selection criteria is calculated on the neighbors, the search progresses to the next stage by expanding the neighbors of the best kernel found.
We consider a generalized notion of the kernel grammar, where we consider a set of base kernels {B 1 , . . . , B r }, r ∈ N and a set of operators {T 1 , . . . , T l }, l ∈ N where T j : K × K → K, j = 1, . . . , l are closed operators on the space of all kernel functions K. Examples are addition and multiplication, but also the change-point operator which was considered in [7]. Thus, in general the grammar operations are 1. Apply operator T j onto a subexpression and a base kernel 1 : S → T j (S, B) 2. Exchange a base kernel with another base kernel: B → B The considered kernel space can be defined precisely in the following way: Definition 1 For r, l ∈ N, let {k 1 , . . . , k r } be a set of (base-) kernels (symbollically represented as {B 1 , . . . , B r }) and {T 1 , . . . , T l } a set of operators with T j : K × K → K, j = 1, . . . , l. Let for i = 1, . . . , M . We callK := L M the kernel-grammar generated kernel space with depth M .

Representation of Kernels
The base kernels and their combination via the operators describe the structural assumptions of the final kernel and imply the assumptions that are made in function space. Our main hypothesis is that the symbolical representation of the kernel, given by the subexpressions S and base kernels, already contains sufficient information for model selection. Therefore, we will define a kernel-kernel over the symbolical representations and utilize it for Bayesian optimization. Concretely, each kernel k ∈K in our described space can be written as a tree T , for example: 1 In case the operator Tj is not symmetric we also add the operation S → Tj(B, S). However, all operators considered in this work are symmetric.
Each operator T j and each base kernel B i is represented by their respective name, where operators are the nodes of the tree and base kernels are the leafs. The way how the operators and base kernels are connected is represented through the tree structure. For a given expression tree T , we denote the multiset of all subexpressions/subtrees as Subtrees(T ). Furthermore, we consider paths to the leafs of the tree, for example, for the tree above one path to a leaf would be: We denote the multiset of all paths in the tree T as Paths(T ). Lastly, we also consider the multiset of all base kernels that exist in a given expression tree as Base(T ). For each described multiset, we denote the number of occurrences of element E in the multiset as n(E). When building the multisets we also account for two symmetries in the elements that can be applied if an operator is associative and commutative, which we elaborate further in Appendix A.
Depending on the operators that are used, several trees can describe the same kernel k ∈K. For technical reasons, we denote with f :K → T (K) a mapping that maps a given kernel k ∈K to one tree T that induces this kernel, where T (K) denotes the set of all expression trees that can generate a kernel inK (see details in Appendix A).

Kernel-Kernel
Our kernel-kernel will be defined via a pseudo-metric over the expression trees. Optimal transport (OT) principles have proven themselves to be effective in BO methods over structured spaces, as shown in [5] and [11] for neural architecture search or in [6] for BO over molecule structures. We follow this line of work and also rely on optimal transport, although we only use a simple ground metric to allow for closed-form computations. To allow OT metrics to be used, we summarize each expression tree T to discrete probability distributions of their building blocks Base(T ), Paths(T ) and Subtrees(T ) via where δ E is the Dirac delta and ω E is the frequency of the element E in the respective multiset. For example, the frequency of expression SE × PER in the multiset Subtree(T ) is calculated as Each probability distribution represents a different modeling aspect that is induced by a kernel k ∈K and its corresponding expression tree T : 1. ω base specifies which base kernels are present in T , thus, which base assumptions in function space are included such as periodicity, linearity, local smoothness.
2. ω paths specifies how the base kernels in T are used, e.g. whether a periodic component is applied additively or multiplicatively.
3. ω subtrees specifies the interaction between the base kernels in T , for example if T contains an addition of a linear and a periodic component or not.
Our pseudo-metric between kernels k 1 and k 2 (and its associated trees T 1 and T 2 ) uses all three modeling aspects via the optimal transport distances between the respective discrete probability distributions ω base , ω paths and ω subtrees .
In general the OT distance with ground metricd between two discrete probability distribution where π ∈ R(ω 1 , ω 2 ) is a combined probability distribution over Ω × Ω with marginal distribution π(A × Ω) = ω 1 (A) and π(Ω × B) = ω 2 (B) for Borel sets A, B. While for general ground metric d : Ω × Ω → R the optimization problem in (1) can not be solved in closed-form we use as ground metricd(E, E ) = 1 E =E which has as closed-form solution the total variation distance between ω 1 and ω 2 (see [19], p. 22), i.e.
Utilizing this distance over the modeling assumptions ω base , ω paths and ω subtrees allows for fast computation of the final pseudo-metric and for a proper positive semi definite (p.s.d) kernel-kernel in the end [see Appendix C]. This is not guarantueed for general OT distances (see [11]). We define the final distance between two kernels k 1 and k 2 (and its associated trees T 1 and T 2 ) as a sum over the OT distances of their respective modeling components where α i ≥ 0 and i α i = 1 are weighting parameters that will later be learned automatically via marginal likelihood maximization for GPs.
In case the kernel grammar contains base kernels that act on single dimensions such as SE i , we define for each dimension i = 1, . . . , D an individual distribution over base kernels where Base(T , i) is the multiset of all base kernels in T defined on dimension i. We also include the empty-expression E NULL for which ω ENULL = 1 whenever no base kernel of dimension i is contained in T . ω (i) base thus summarizes which base kernels are present that act on dimension i or if the dimension is ignored. The distance in this case is defined as +α 2 Wd(ω 1,paths , ω 2,paths ) +α 3 Wd(ω 1,subtrees , ω 2,subtrees ).

(3)
The reason for this distinction is that two kernels k 1 and k 2 , which have different active dimensions, are considered to be farther apart with this distance, making it easier to allow variable selection, which is an important aspect of model selection.
Our defined function d(T 1 , T 2 ) induces indeed a pseudo-metric in the kernel-grammar generated kernel space, as shown in the following proposition. In particular, it fulfills the triangle inequality.
Proposition 1 LetK be the kernel space generated by a kernel grammar. Let f :K → T (K) be a mapping that maps a kernel k ∈K to one of its expression trees T . Thend(k 1 , is a pseudo-metric onK where d is given by (2) or (3) depending on the base kernels inK.
Given the pseudo-metric, we are now able to define a kernel-kernel with which we call Symbolical-Optimal-Transport (SOT) kernel-kernel. Here, l denotes the lengthscale and σ 2 the variance of the kernel-kernel. Both parameters are learned in combination with the distance weights via marginal likelihood maximization (see Appendix A).

Bayesian Optimization for Model Selection
We utilize the proposed kernel-kernel to do model selection for GP's via Bayesian optimization, which is a similar task to [9]. Compared to [9], we use our proposed kernel-kernel, which is a fundamentally different and computationally more efficient way of measuring similarity in GP space. Given a model selection criteria g D :K → R for a given dataset D, we want to solve k * = arg max k∈K g(k|D) via Bayesian optimization. We thus define a surrogate GP model for g(k|D) via where µ c (k) = c is the constant mean function. As BO acquisition function a(k|D t ) we use Expected-Improvement (EI). Here,D t denotes the set of already queried kernel-selection-criteria pairs (k, g(k|D)) for which the meta-GP posterior f |D t is calculated. Starting with an initial set of pairsD 0 we follow the standard BO iterations where in each iteration the acquisition function is maximized k t = arg max k∈K a(k|D t ) given the current set of already evaluated kernel-selectioncriteria pairsD t . Then the selection criteria is queried at the chosen kernel g(k t |D). A complete description of the method can be found in Algorithm 1 in Appendix A.
A crucial part in the BO cycle is the optimization of the acquisition function max k∈K a(k|D t ). While for Euclidean spaces gradient-based methods or grid-based methods could be used to solve that task, this is not an option for a structured space like the considered grammar-generated kernel space. We propose using an evolutionary algorithm to optimize the acquisition function, which can be seen in Algorithm 2 in Appendix A. While an evolutionary algorithm seems to be a computationally intense procedure that takes place in each BO iteration, we emphasize that the evaluation of the acquisition function is very efficient for our proposed method. The reason for the efficiency is that the evaluations of our kernel-kernel K SOT (k 1 , k 2 ) at two kernels k 1 , k 2 ∈K is very cheap, compared for example to the method in [9] (see Appendix E for computational time comparision).

Comparision to Hellinger Kernel-Kernel
The work of [9] also uses BO for kernel selection but with a different principle of measuring the distance in GP space. In particular, they propose measuring the distance of two GP's M and M via the induced prior distributions p(f |X, M) in function space evaluated on the design matrix X of the dataset. Conditioned on the parameters of the GP θ, the distributions p(f |X, M, θ) are Gaussian and they use the Hellinger distance over Gaussians as base distance between GP models d(M θ , M θ |θ, θ ) given the hyperparameters. They finally take the expectation over the hyperparameter priors to construct their final distance. Thus, for each kernel-kernel evaluation K(M, M ) the integral needs to be computed, where the integrand scales cubically in |X|. Each kernel-kernel evaluation is therefore a computationally hard problem by itself. They side-pass this issue partly by selecting a subset of X as input locations and by using sample-based estimation of the integral. While this could make the method useful for oracles with very long run times, such as model selection for very large data sets, we see conceptual problems in the case of medium oracle run times. For our method, kernel-kernel evaluations K(k 1 , k 2 ) are very efficient as they scale only in the size of the expression tree and no integrals need to be computed.

Experiments
In the following section, we show experimental results for our novel meta-GP model and kernel search method. We evaluate our meta-model on a meta-regression task where we predict test kernellog-evidence pairs based on training pairs. Secondly, we consider kernel search and compare it against the greedy method in [3], the evolutionary algorithm in [4] and against the BO method of [9]. The implementation of our method is available at https://github.com/boschresearch/bosot.

Selection Criteria and Datasets:
In all our experiments, we consider the normalized log-model evidence g(k|D) = log p(y|X, k)/|D| as model selection criteria, as also done in [9]. We always use the Laplace approximation to calculate the log-model evidence, where we use 10 repeats to do MAP estimation of the kernel parameters. Furthermore, we consider the following publicly available datasets: Airline, LGBB, Airfoil, Powerplant, Concrete. Airline is a one dimensional time-series dataset, LGBB is a two-dimensional dataset with low observation noise. Powerplant, Airfoil and Concrete are four, five, and eight dimensional, whereas Powerplant exhibits higher observation noise. We use medium-sized training sets -100 and 150 datapoints for Airline and LGBB and 500 datapoints for the other three datasets. All training sets are uniformly drawn from the full datasets. The outputs are normalized and the inputs are scaled to be in the unit interval. Further details can be found in Appendix D.
Search Spaces: We consider two search spaces. The first consists of the base kernels SE i , LIN i , PER i , RQ i and the operators + and ×. This is the space considered in [3] for time series and low-dimensional datasets. The second search space uses as base kernels SE i and RQ i and also + and × as operators and was considered in [3] for higher dimensional base datasets. We consider the first space for Airline and LGBB and the second for Powerplant, Airfoil, and Concrete. The hyperparameter priors for the base kernels can be found in the Appendix A.
Prediction of Selection Criteria: We evaluate our meta-GP model on a meta-regression task by predicting test kernel-log-evidence pairs based on training pairs. The quality of predictions in kernel space might also be an indicator for good performance in BO for kernel search. Given the base dataset D we create a training setD train = {(k i , g(k i |D))|k i ∈K train ⊂K, i = 1, . . . , n train } and a test set . . , n test } -each containing 500 kernel-log-evidence pairs. We generate the train and tests sets exactly as in [9], where we create one setK complete by first initializing it with all base kernels and iteratively pick one kernel of the current set, apply one random operation of the kernel grammar and add the resulting kernel to the current set and repeat this process until we have n train + n test kernels inK complete . We then divide the set uniformly intõ K train andK test . We compare our model to the mean-predictor as baseline, which just predicts the mean of the train set at each test point and a kNN predictor based on the kernel-grammar operations [see Appendix D or [9]]. Furthermore, we compare against the meta-GP model of [9]. We report root mean squared error (RMSE) scores on the test sets in Table 1. It can be observed that our method leads to more precise predictions on all four meta prediction tasks, indicating that the symbolical representations already contain sufficient information to predict log-evidence values.  Model Selection -Setup: Concerning kernel search, we compare against the BO method that employs the Hellinger kernel-kernel [9], against greedy search [3] and against the evolutionary algorithm presented in [4], referred to as TreeGEP. Both BO methods run for 50 iterations and the kernel-kernel hyperparameters are updated in each iteration via marginal likelihood maximization. Our method applied the evolutionary Algorithm 2 (see Appendix A) to optimize its acquisition function, using a population size of 100. As it is computationally unfeasible to apply the same kind of acquisition function optimization for the Hellinger kernel-kernel we use their method of optimizing the acquisition function where an active set of kernels is updated in each iteration. Both BO methods and TreeGEP start with the same set of initial kernels, for which we apply two random grammar operations for each base kernel. Greedy search by design needs to start from the empty kernel. We give it a head start by the number of initial datapoints (see Appendix D). For each dataset, we display results from 30 independent runs with different seeds, namely medians and quartiles of the log-evidence score over iterations and CPU-time in Figure 1. The implementation of both BO methods is based on GPflow [10]. Further experimental details and parameter settings for all methods can be found in Appendix D.
Model Selection -Results: As shown in Figure 1, we outperform all methods in terms of performance over number of model evaluations. This is not surprising against the two heuristics, as they are not optimized towards keeping the number of model evaluations low. However, it is notable that we are more sample-efficient compared to [9], who solve a much harder problem for computing their kernel-kernel. In terms of performance over CPU-time, we outperform all methods on Airfoil, Powerplant, LGBB and Concrete significantly. On Airline, the advantage over the heuristics is smaller -the reason is the lower oracle time -which benefits the heuristics that do not need to optimize the acquisition function. The reason for the poor performance of the Hellinger kernel-kernel in terms of CPU-time is the high ratio of acquisition function optimization to oracle time, which was as high as 50 : 1 in our experiments (see detailed numbers in Appendix D).
Test Performance: When optimizing a model selection criteria, one expects that this also materializes in a better test performance. In Table 2 we therefore show the predictive negative log-likelihood (NLL) scores on a held out test-set of the selected models at the final time stamp. We observe that the advantage in the model selection value gets transferred to the test performance.
Further Experiments: In Appendix E we include further investigations on the behavior of our method on simulated data coming from a ground-truth kernel. Furthermore, we include a comparison of the selected kernel to the standard RBF kernel and against Functional Kernel Learning (FKL) [2].

Limitations
In case the dataset size is very small, greedy search or evolutionary algorithms might have an advantage in terms of computational time as the optimization of the acquisition function outweighs the computation of the model selection criteria. Thus, in these instances our method might not show a strong benefit. Furthermore, the kernel grammar is often used for downstream applications that use the selected hypothesis for interpretation, such as [7] who build an automatic natural language description of the selected hypothesis. Depending on the dataset and how many steps are employed in the acquisition function optimization, relatively large hypotheses might be found as optimal. This might render the interpretation difficult. However, we show in Appendix E example configurations of our algorithm that can be used to get smaller, well interpretable hypothesis.

Conclusion
We presented a novel way of doing BO for model selection of GP's by measuring the distance between two GP's via the symbolical description of the underlying statistical hypothesis. The main contribution is the deduced pseudo-metric over kernels and the resulting SOT kernel-kernel. We show that our approach leads to a more efficient way of searching through a discrete kernel space compared to other BO methods and search heuristics.

A Method Details
In the following section, we specify further details of the proposed method. We show the symmetries we employ in the multisets and give details on the approximation of the log-evidence. Furthermore, we describe the base kernels and their parameter priors, the marginal likelihood maximization of the kernel-kernel hyperparameters and how the acquisition function optimization is done.

A.1 Symmetries in the Multisets
Our method can be used with general operators T : K × K → K in kernel space. However, depending on the concrete operators, one might also incorporate properties that leave the expression/subexpressions unchanged. In particular, addition and multiplication are commutative and associative. For commutative operators two trees/subtrees describe the same kernel if one rotates the nodes under a commutative operator, e.g.

PER SE
This symmetry is automatically considered in the multisets Base(T ) and Path(T ) as the base kernels and paths stay unchanged when rotating two nodes. For Subtree(T ), we account explicitly for that symmetry via hashing functions that are invariant to rotation of subtrees. This allows counting tree S 1 and tree S 2 as an identical subtreeS in the multiset Subtree(T ).
For the multiset Path(T ), we furthermore consider a symmetry that exists in case an operator is associative and commutative. In this case, one can exchange the nodes of two consecutive operators of the same kind without changing the expression, e.g. We incorporate this symmetry into Path(T ) by only considering the smallest path to that base kernel that exist in an equivalent expression (equivalent under this symmetry). This is realized in Path(T ) by counting identical operators in a row along a path only as one, e.g. the following paths would be considered the same: While one might also integrate more symmetries that stem from, e.g. the distributive property of the multiplication, we found that these two symmetries are particularly easy and efficient to implement. The symmetries are also not restricted to multiplication or addition. The first symmetry can be used whenever T is associative, thus, T (k 1 , k 2 ) = T (k 2 , k 1 ) and the second whenever T is associative and commutative, thus, T (k 1 , T (k 2 , k 3 )) = T (k 2 , T (k 1 , k 3 )) = T (T (k 1 , k 3 ), k 2 ). The motivation of considering symmetries at all is that two trees T 1 and T 2 are considered more similar in case they share the same structure (with respect to the symmetry) and, therefore, similar kernels can be detected more easily.

A.2 Mapping from Kernels to Trees
As we already observed in the previous section, two different trees T 1 and T 2 might describe the same kernel k. For technical reasons, we therefore consider in Proposition 1 a mapping f : K → T (K) that maps a kernel always to the same tree. Concretely, this is done such that the pseudo metric is defined in kernel space rather than in tree space. When using addition and multiplication, one could check if T 1 and T 2 describe the same expression via recursive hashes that follow the same rules as addition and multiplication. One could use this check to implement such a function f : K → T (K). However, in our experiments we directly deal with the trees (the BO and the acquisition function algorithm directly act on the trees anyway) and ignore filtering out potential collisions of two treeswe did not observe any downsides of doing that.

A.3 Approximation of Log-Model-Evidence
In all our experiments, we use the normalized log-model-evidence g(k|D) = log p(y|X, k)/|D| as selection criteria. Similar to [9] we use the Laplace approximation to approximate the log-evidence, which is (see [9]) where γ = (θ, σ 2 ) ∈ R d are the parameters of the kernel and the likelihood variance, d ∈ N is the number of parameters,γ denotes the MAP estimate of γ, and Σ −1 = −∇ 2 log p(γ|D, k)| γ=γ . Creating the MAP estimate of γ scales cubically in the dataset size |D| in each optimization step. We use the LBFGS optimizer to create the MAP estimate. As the loss function is non-convex, we make 10 restarts with random initialization of the initial parameters. Fit Meta GP model f ∼ GP(µ c (·), K SOT (·, ·)) toD t 5 k t ← EvolutionaryAlg(a(·|D t ),L) 6 Query model selection criteria t * ← arg max t=0,...,T −1 g t 10 return k t * Algorithm 1: BO for kernel search via SOT kernel-kernels.

A.6 BO Algorithm and Acquisition Function Optimization
In Algorithm 1, we show the BO steps for kernel search. First, we draw an initial dataset of kernelselection-criteria pairs. Our base setting applies two random grammar operations from each base kernel. The meta-GP parameters are fitted in each BO optimization. The acquisition function is optimized via the evolutionary algorithm in 2. This algorithm searches in the hypotheses space from small to big hypotheses. Given an initial population of kernels, we calculate the acquisition function on all kernels in the population. Then the n survive best kernels in the population survive. Each selected kernel gets n offspring children by generating n offspring new kernels via applying one random grammar operation. The new population is formed via the selected kernel combined with its offspring. L is the number of iterations in the evolutionary algorithm. In combination with the initial population it determines how many base kernels the hypotheses can contain maximally. In each iteration, the number of base kernels in the best hypothesis can maximally grow by one.

B Example Calculation of SOT Kernel-Kernel
In this section, we give a small example calculation of our proposed kernel-kernel. We calculate our proposed pseudo-metric d SOT (k 1 , k 2 ) for the two kernels k 1 and k 2 with the following expression trees In the first step, we extract the tree features and create the multisets. Here, we use the notation {(E 1 ; n(E 1 )), (E 2 ; n(E 2 )), . . . } to denote the existence of an element E in the multiset as well as the cardinality of the element. We obtain the following multisets: (SE; 2), (PER; 1)}. Next, we build the probability vector for each feature multiset, for example, for the base kernels: and calculate the total variation distance between ω 1,base and ω 2,base : Wd(ω 1,base , ω 2,base ) = 1 2 For the paths, the total variation distance results to: Wd(ω 1,path , ω 2,path ) = 1 2 The complete distance given the weights is then d(k 1 , k 2 ) = α 1 3 20 + α 2 11 20 + α 3 11 21 and The parameters are learned via marginal likelihood maximization and are dependent on the dataset.

C Technical Details and Proofs
When referring to a kernel k from the kernel-grammar generated kernel spaceK such as writing k ∈K, we actually refer to the associated kernel family over its parameters {k θ |θ ∈ Θ}. Applying an operator T onto two kernels k 1 and k 2 then results in another kernel family {T (k 1,θ1 , k 2,θ2 )|θ 1 ∈ Θ 1 , θ 2 ∈ Θ 2 } and we refer with the notation T (k 1 , k 2 ) to this kernel family with new parameter space Θ 1 × Θ 2 . In case the same kernel family {k θ |θ ∈ Θ} appears twice in an operator the parameters are not shared, e.g. T (k,k) corresponds to the family {T (k θ1 ,k θ2 )|θ 1 ∈ Θ, θ 2 ∈ Θ} with new parameter space Θ × Θ.

Definition 2
We call a binary tree T whose nodes are associated with the operators {T 1 , . . . , T l } and whose leafs are associated with the base kernels {B 1 , . . . , B r } an expression tree of a kernel k if k is constructed by applying the operators recursively onto the leafs, meaning k is the result of Algorithm 3.
Any kernel k in a kernel-grammar generated kernel space has an expression tree T via the construction ofK. We further denote the set of all trees that can generate a kernel inK as T (K) := {T |∃k ∈K : k is result of GetKernel(T )}.
Proof of Proposition 1: W.l.o.g. we considerK without separate base kernels for dimensions. We denote by g i : T (K) → [0, 1] Li and i ∈ {base, path, subtree} the mappings from expression trees to the probability vectors ω base , ω paths , ω subtrees , where L i denotes the number of different elements of the respective type (e.g. number of different paths to the leafs for expression trees of depth M ). Thend is a pseudo metric as chaining of a metric with a mapping results in a pseudo metric, and the positive-weighted sum of pseudo metrics still is a pseudo metric.
Proposition 2 LetK be the kernel space generated by a kernel grammar. Then is a proper p.s.d. kernel overK.
Proof: As in the proof of Proposition 1 we can writed(k 1 , k 2 ) as a weighted sum of Manhatten metrics, which leads to: Thus, K SOT (k 1 , k 2 ) can be written as a product of Ornstein-Uhlenbeck kernels chained with functions g i • f . Kernels chained with arbitrary mappings are kernels and products of kernels are kernels (see [16], Proposition 3.22). Thus, K SOT (k 1 , k 2 ) is a proper (p.s.d.) kernel.

D Experimental Details
In this section, we give further details on parameter configuration and implementation details of the different methods.  Hellinger Kernel-Kernel: For the Hellinger Kernel-Kernel [9], we use their principle of optimizing the acquisition function, where an active set of kernels is kept in memory over the iterations. In each iteration, 15 random walks are performed in the kernel grammar, where the walk length is drawn randomly from a geometric distribution with p = 1 3 . These kernels are added to the active set. Furthermore, the neighbors of the best kernel found so far are added to the active set (we limit this set to 50 neighbors per iteration, since for large expressions there may be several hundred neighbors, making it computationally almost infeasible to evaluate the kernel-kernel on all of them). The active set is limited to 600 kernels, where only the ones are kept with the highest acquisition function value. We initialize the active set with random walks in the grammar, with one random walk of length 5 from each base kernel (this is an alternative to their computationally very expensive version of using all kernels two edges apart from the base kernels). We stick to their base settings of using 100 hyperparameter samples and a subset of the design matrix X of size of 20 inside the Hellinger distance calculation. We cache all distance calculations d(M, M ) over the iterations, thus, recalculation of the Hellinger distance for kernels in the active set or the current dataset is very cheap. On the other hand, calculating the distance that involves a new kernel, such as a new neighbor of the currently best kernel, is very expensive. Given the distance matrix between kernels, optimization of the kernel-kernel hyperparameters is cheap again.
Greedy Search: The greedy method in [3] starts with the empty kernel, then evaluates all base kernels. It picks the best performing base kernel and determines its immediate neighbors via expanding the base kernels with all possible grammar operations. Then the neighbors are evaluated (we pick the order in which the neighbors are evaluated at random -for each seed a different order). Once all neighbors have been evaluated, the best kernel is determined and the process repeats. We give greedy search a head start of n initial log-evidence evaluations in Figure 1, meaning log-evidence values are shown once greedy search has evaluated as many kernels as the other methods have in their initial dataset.
TreeGEP: For the evolutionary algorithm in [4], we use their base settings in the paper, which is a population size of 200, a reproduction rate of 0.1, and a probability of mutation vs. cross-over of 50% each. They don't specify the size of the mutation subtrees and the tournament fraction in the tournament selection. Here, we make a reasonable choice and generate mutation subtrees of size four and use a tournament fraction of 0.1.
Test Performance -Box-Plots: In Figure 2, we show box-plots of the RMSE and NLL on the held out test-sets for the selected kernels. We see that on both metrics, our search method also finds a final model that often yields the best test-performance or is among the best models. However, we also note that test performance is mainly a property of the model selection criteria and how well models that maximize that criteria generalize to new data points.
Acquisition Function Optimization vs. Oracle Evaluation Ratio: In Table 3, the ratio tAcquisition t Oracle between the CPU-time needed for optimization of the acquisition function and the CPU-time for evaluating the oracle (calculating the log-evidence) is shown for the SOT kernel-kernel and the Hellinger kernel-kernel [9]. As described in the previous section, the approach in [9] uses a different kind of acquisition function optimization, that has fewer evaluations of the acquisition function than the evolutionary algorithm we use. Nevertheless, we are magnitudes of orders faster. This can be understood more clearly when looking at the raw kernel-kernel evaluation times in Figure 3. A faster acqusition to oracle time ratio in the end results in a faster search procedure, measured over CPU-Time (as we show in our main results in Figure 1), which is the metric we are interested in.
kNN for Meta Prediction: For the k-nearast-neighbour approach that was used in Table 1 we employ the same principle as in [9] who also consider kNN for comparision. The set K complete    Figure 3.
Kernel-Index 1 2 3 4 5 # Parameters 2 5 8 10 12 # Base Kernels 1 2 3 4 5 consists of kernel expressions which form a directed graph where two neighbour nodes are one grammar operation apart from each other. For each test expressionk ∈ K test , we search in this graph the k expressionsk 1 , . . . ,k k ∈ K train with shortest path in this directed graph. The prediction of the log-evidence value ofk is the average of the log-evidence values ofk 1 , . . . ,k k . The number k of neighbours is determined via cross validation.
Runs and Time Stamps: In Figure 1, all methods were repeated over 30 seeds, where each seed corresponds to a different initial dataset for the both BO methods and TreeGEP and a different ordering of neighbor evaluations in greedy search. All runs have different run-times. The reasons for this is that the oracle evaluation times differ depending on the kernel that is evaluated. The log-evidence can be calculated faster for smaller hypothesis. The final time stamps in Figure 1 are therefore determined by the shortest run of our method, as we only have log-evidence values for all runs/seeds up to that time point. In rare cases, a run can be interrupted, in case the Laplace approximation returns NaNs. This can happen due to numerical instabilities in the Cholesky decomposition and happened independently of the search method. We filtered out these runs.

E Further Experiments
Computational Time for Single Kernel-Kernel Evaluations: In Figure 3, we show CPU-time for single kernel-kernel evaluations K(k 1 , k 2 ) for the SOT kernel-kernel and the Hellinger kernel-kernel, each evaluated on five kernels generated from the kernel grammar (search space was the same as used for the LGBB dataset). We selected the five kernels with increasing numbers of base kernels and, thus, kernel parameters, as this could affect the computation time (see Table 4). Our kernel-kernel can be evaluated orders of magnitudes faster, which also explains the smaller acquisition function optimization times.
Kernel-Kernel Hyperparameters: In Figure 4, we show the values of the distance weights for the three OT metrics Wd(ω 1,base , ω 2,base ) (on Airfoil summed over dimensions), Wd(ω 1,paths , ω 2,paths ) and Wd(ω 1,subtrees , ω 2,subtrees ) over the BO iterations on Airline, and Airfoil. The weights were fitted via marginal likelihood maximization as described in A.5. First, we observe that all OT metrics are used. Secondly, we see that it is dependent on the dataset which OT metric is used primarily. For example, on Airline the Subtree features are the most important ones (according to the marginal likelihood maximization) whereas on Airfoil the Base features obtain the largest weights.
Interpretability of Final Hypothesis: The kernel grammar is used in a line of work called Automatic Statistician [3,7]. In particular [7], utilize the kernel grammar to automatically generate natural language description of the data from the selected hypothesis (they also employ greedy search).
Here, we show a configuration of our algorithm where the method of [7] can be applied to the final hypothesis of our model selection procedure and can give good descriptions of the dataset. As done in [7], we drop the rational quadratic kernel from the search space, as this captures small and long range correlations at the same time, which can also be modelled with separate squared exponential kernels (see [7]). Furthermore, we take fewer steps in the acquisition function optimizer such that fewer base kernels are maximally possible in the final hypothesis. This renders the final hypothesis smaller and easier to interpret, as fewer components needs to be described. We show two example hypotheses that were selected for the Airline dataset using the described search space and four steps in the acquisition function optimizer. We used the principles presented in [7] to simplify the expression and to generate the sentences.

LIN + SE + LIN × SE + PER × LIN × SE
Desciption -Example 1: The data can be described as a sum of:  Comparision with RBF Kernel: In Table 5 the test-set results (RMSE and NLL) of our proposed kernel search method is shown after 50 iterations compared to test values of a standard RBF kernel.
We note that Table 5 shows different NLL values as Table 2 as it shows values at last iteration and not at last time stamp (see Appendix D for details). In terms of NLL, the RBF kernel shows competitive performance on the LGBB and Concrete dataset. However, considering the RMSE and NLL values on the other three datasets, it appears that kernel selection in general seems to be very important.
Comparision with Nonparametric Kernel Learning Methods: In Table 5 we also compare with the nonparametric kernel learning method, presented in [2], called Function Kernel Learning (FKL). FKL places a GP prior on the spectral density of the kernel -thus utilizing a nonparametric approach to kernel learning/selection. This results in a prior over spectral-mixture kernels, thus a prior over a fixed but highly flexible kernel structure. The results in Table 5 might be an indication that search over a discrete set of structural kernels might be beneficial, compared to placing a prior over a very flexible, but fixed kernel family.
Experiments on Simulated Data & Type-3 Maximum Likelihood Overfitting: In Figure 5, we show experimental results on simulated data. Here, we draw n datapoints from a GP prior with Figure 5: Learning-curves of our kernel selection method (over iterations) for simulated datasets. Upper plot shows log-evidence values, lower plot shows test-set RMSE and NLL values. Blue lines mark log-evidence/test values of the ground-truth kernel, from which the dataset was generated. a given kernel structure and employ our kernel selection method on that dataset (likelihood noise was 0.01 for all datasets). First, we observe that we reach (almost) the same log-evidence values as the ground truth kernel within 50 iterations. For the RBF ground-truth kernel we simulated one big dataset and used n = 100 and n = 300 for model selection in the third and fourth plot. For the smaller dataset, we observe that we find kernels that have even higher log-evidence values than the ground-truth kernel. Considering the test-set NLL we also observe a small increase in the NLL from iteration 10 to 30, indicating a small overfitting. We think that in particular for smaller datasets maximizing the log-evidence over kernel structures can also lead to overfitting. However, this might not be surprising as the same was observed for type-2 maximum likelihood in GP's (see. [12]) whereas maximizing the log-evidence might be interpreted as type-3 maximum likelihood. Possible methods to avoid overfitting could be to use a different model selection criterion such as the Bayesian information criterion or cross-validation error, or to use a smaller search space.