Kernel Functional Optimisation
Year: 2022

Abstract:
Traditional methods for kernel selection rely on parametric kernel functions or a combination thereof and although the kernel hyperparameters are tuned, these methods often provide sub-optimal results due to the limitations induced by the parametric forms. In this paper, we propose a novel formulation for kernel selection using efficient Bayesian optimisation to find the best fitting non-parametric kernel. The kernel is expressed using a linear combination of functions sampled from a prior Gaussian Process (GP) defined by a hyperkernel. We also provide a mechanism to ensure the positive definiteness of the Gram matrix constructed using the resultant kernels. Our experimental results on GP regression and Support Vector Machine (SVM) classification tasks involving both synthetic functions and several real-world datasets show the superiority of our approach over the state-of-the-art.

heading: 1 Introduction
text: Kernel machines (Hofmann et al., 2008) generally work well with low-dimensional and small to medium-scaled data. In most kernel machines, the kernel function is chosen from the standard bag of popular kernels (Genton, 2001, Stein, 2015) such as Squared Exponential kernel (SE), Matérn kernel and Periodic kernel, or a weighted combination thereof (Aiolli and Donini, 2015, Gönen and Alpaydın, 2011, Rakotomamonjy et al., 2007). Recent developments (Jang et al., 2017, Wilson and Adams, 2013) in kernel learning parameterise the kernel function to boost the expressiveness of the kernel. However, the expressiveness of such kernels remains limited by the chosen parametric form and thus they often fall short in providing the best kernel function for complex data distributions.
There have been some early attempts to design an optimal non-parametric kernel to remove the limitations associated with the parametric forms. Ong et al. (2003, 2005) proposed a hyperkernel framework by defining a Reproducing Kernel Hilbert Space (RKHS) on the space of kernels i.e., a kernel on kernels to support kernel learning. They formulate a semidefinite programming (Vandenberghe and Boyd, 1996) based optimisation problem using the representer theorem (Steinwart and Christmann, 2008, Vapnik, 1999) to find the best kernel. However, their method suffers from two key limitations: (i) their way of enforcing the positive definiteness property produces a restrictive search space, resulting in a sub-optimal solution, and (ii) the computational complexity of their method scales with the dataset size, making it infeasible for larger datasets. Benton et al. (2019) proposed Functional Kernel Learning (FKL), which extends the function space view of the Gaussian Process (GP) for kernel learning. FKL uses a transformed GP over a spectral density to define a distribution over kernels. However, the formulation of kernel functionals using the spectral densities induces strong assumptions on the properties such as periodicity, stationarity, etc. and thus are not generally applicable. Malkomes et al. (2016) proposed an automated kernel selection (BOMS) using Bayesian optimisation. The kernel space in BOMS is defined by the base kernels and the associated grammar to combine them. Although the search space is constructed by summing or multiplying the base kernels, the resultant kernel space is restricted in the compositional space of parametric forms.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
In this paper, we propose a generic framework called Kernel Functional Optimisation (KFO) to address the aforesaid shortcomings. First, it provides a flexible form of kernel learning whose computational complexity is decoupled from dataset size. Next, it allows us to use a computationally efficient Bayesian optimisation method to find the best kernel. We incorporate hyperkernels into our Bayesian framework that allows us to search for the optimal kernel in a Hilbert space of kernels spanned by the hyperkernel (Ong et al., 2005). We draw kernel functionals from a (hyper) GP distribution fitted using a hyperkernel. As the kernel drawn from the hyper-GP may be indefinite, we provide ways to ensure positive definiteness by transforming indefinite, or Kreı̆n (Oglic and Gärtner, 2019, Ong et al., 2004) kernel space into a positive definite kernel space. The optimisation of kernel functionals necessitates solving larger covariance matrices and thus adds to the computational burden of the overall process. To speed up the computations, we perform a low-rank decomposition of the covariance matrix. Further, we provide a theoretical analysis of our method showing that it converges efficiently as in its cumulative regret grows only sub-linearly and eventually vanishes.
We evaluate the performance of our method on both synthetic and real-world datasets using SVM classification (Diehl and Cauwenberghs, 2003, Scholkopf and Smola, 2001, Burges, 1998) and GP regression tasks. Comparison of predictive performance against the state-of-the-art baselines demonstrates the superiority of our method. Further, we compare with the state-of-the-art performance reported in the latest survey paper on classifier comparison (Zhang et al., 2017) and find that our method provides the best performance on most of the datasets. Our main contributions in this paper are as follows: (i) we propose a novel approach for finding the best non-parametric kernel using hyperkernels and Bayesian functional optimisation (Section 3), (ii) we provide methods to ensure positive definiteness of the kernels optimised (Section 3), (iii) we derive the convergence guarantees to demonstrate that the regret grows sub-linearly for our proposed method (Section 4), (iv) we provide empirical results on both synthetic and real-world datasets to prove the usefulness (Section 5).

heading: 2 Background
text: Notations We use lower case bold fonts v for vectors and vi for each element in v. vᵀ is the transpose. We use upper case bold fonts M (and bold greek symbols) for matrices and Mij for each element in M. | · | for the absolute value. Nn = {1, 2, · · · , n}. R for Reals. X is a non-empty (index) set and x ∈ X . X̃ is a non-empty (compounded index) set and x̃ ∈ X̃ , X̃ = X 2. (·)+ clips a negative value to zero. J·K is the Iverson bracket (Iverson, 1962) defined for any boolean value I as JIK = 1, if I is True, 0 otherwise. Matrix M = [Mij ]i,j∈N and ‖M‖F is the Frobenius Norm of M.

heading: 2.1 Bayesian Optimisation
text: Bayesian Optimisation (BO) (Brochu et al., 2010, Shahriari et al., 2015, Frazier, 2018) offers an elegant framework for finding the global extrema of an unknown, expensive and noisy function f(x), represented as x∗ = argmaxx∈X f(x), where X is a compact search space. Bayesian optimisation is comprised of two main components: (i) a Gaussian Process (GP) (Williams and Rasmussen, 2006) model of f , and (ii) an acquisition function (u) (Kushner, 1964, Močkus, 1975, Wilson et al., 2018) to guide optimisation. Let D = {x1:t,y1:t} denote a set of observations of f , where y = f(x) + ′ is the noisy observation corrupted with white Gaussian noise ′ ∈ N (0, σ2noise). Then the predictive distribution at any point x∗ is given as f(x∗)|D ∼ N (µ(x∗), σ2(x∗)), where µ(x∗) = kᵀ[K + σ2noiseI]
−1y1:t, σ2(x∗) = k(x∗,x∗)− kᵀ[K + σ2noiseI]−1k, k =[k(x∗,x1) · · · k(x∗,xt)], k : X × X → R and K = [k(xi,xj)]i,j∈Nt . The negative log-likelihood for a GP distribution is
− logP(y∗|D,x∗)= 12 log(2πσ2(x∗)) + (y∗−µ(x∗))2 2σ2(x∗) (1)
The acquisition function (u) guides the search by balancing between exploitation (searching known high-value regions) and exploration (searching high-variance regions). Gaussian Process - Upper Confidence Bound (GP-UCB) acquisition function (Srinivas et al., 2012, Brochu et al., 2010) is the commonly used acquisition function to find the next best candidate for the evaluation, given as
ut(x) = µ(x) + √ βt σ(x) (2)
where βt grows as O(log t) with iteration t. Further, it can be shown that the average regret (R , 1t ∑t t′=1 |f(x∗)− f(xt′)|) grows as O( √ log t/t), and hence the average regret vanishes as t→∞. An algorithm for standard Bayesian optimisation is provided in the supplementary material.
The aforementioned standard Bayesian optimisation procedure often suffers from scaling issues originating from the curse of dimensionality. Wang et al. (2016) proposed REMBO - Random EMbedding Bayesian Optimisation - to address these scaling issues. REMBO works by projecting the objective function onto a lower-dimensional subspace prior to optimisation. LINEBO (Kirschner et al., 2019) builds on the same idea but instead of a fixed subspace, it decomposes the given black-box optimisation problem into a sequence of one-dimensional subproblems. Further, our method builds upon the principles of Bayesian functional optimisation methodologies (Vien et al., 2018, Vellanki et al., 2019, Shilton et al., 2020) in the literature to find a function to optimise the given process.

heading: 2.2 RKHS and Hyper-RKHS
text: The kernel functions used in the Gaussian process uniquely define an associated Reproducing Kernel Hilbert Space (RKHS) (Aronszajn, 1950). Formally:
Definition 1: LetHk be a Hilbert space of functions f : X → R on a non-empty set X . A function k : X × X → R is a reproducing kernel of Hk, and Hk a Reproducing Kernel Hilbert Space (RKHS), if the following properties are satisfied.
• k spansHk i.e.,Hk = span{k(·,x)|x ∈ X} • ∀x ∈ X , ∀f ∈ Hk, 〈f(·), k(·,x)〉Hk = f(x) (the reproducing property) • ∀x, x′ ∈ X , k(x,x′) = 〈k(·,x), k(·,x′)〉Hk
Next, we consider the Reproducing Kernel Hilbert Space (RKHS) of kernels by introducing a compounded index set X̃ : X × X and a hyperkernel κ (Ong and Smola, 2003, Ong et al., 2003). Analogous to the RKHS (Aronszajn, 1950) associated with the kernel function, a hyperkernel defines an associated Hyper-Reproducing Kernel Hilbert Space (Hyper-RKHS) (Ong et al., 2003).
Definition 2: Let X be a non-empty set and X̃ denote X × X . The Hilbert space Hκ of functions k : X̃ → R is called a Hyper-Reproducing Kernel Hilbert Space (Hyper-RKHS), if there exists a hyperkernel κ : X̃ × X̃ → R that satisfies the following properties:
• κ spansHκ i.e.,Hκ = span{κ(·, x̃) | x̃ ∈ X̃} • ∀x̃ ∈ X̃ , ∀k ∈ Hκ, 〈k(·), κ(·, x̃)〉Hκ = k(x̃) (the reproducing property) • ∀x̃, x̃′ ∈ X̃ , κ(x̃, x̃′) = 〈κ(·, x̃), κ(·, x̃′)〉Hκ • κ(x′,x′′,x′′′,x′′′′) = κ(x′′,x′,x′′′,x′′′′) ∀x′,x′′,x′′′,x′′′′∈X
The GP distribution defined by a hyperkernel κ is a distribution on the space of kernels. This Hyper-RKHS is a Hilbert space comprised of positive definite, negative definite and indefinite kernels. A Kreı̆n kernel k (Oglic and Gärtner, 2018, Ong et al., 2004) is an indefinite kernel with a positive decomposition i.e., there exist positive kernels k+ ∈ H+ and k− ∈ H−, such that k = k+ − k−. From Definition 2, we see that κ(x̃, x̃′) = κ(x′,x′′,x′′′,x′′′′) is a kernel, where x̃ = (x′,x′′). Generally, the samples drawn from GP(0, k) do not lie in the corresponding RKHS Hk, but in a larger RKHSHk′ 6=k (see discussion in Kanagawa et al. (2018), Remark 3.8 and Section 4). We also note that the posterior mean of GP(0, k) lies in the RKHS Hk. Similarly, with hyperGP, the samples drawn from GPκ(0, κ) lie in RKHS Hκ′ 6=κ, whereas its posterior mean (µ) lies in Hκ. Further, µ can be decomposed with positive and negative weights as µ = µ+ − µ− =∑ i αi+κ(·, x̃i+) − ∑ i αi−κ(·, x̃i−), where αi+ , αi− > 0; and µ± = ∑ i αi±κ(·, x̃i±) is a kernel (Definition 2 and Ong et al. (2004)). Thus, µ = µ+−µ− is a Kreı̆n kernel (Oglic and Gärtner, 2019).

heading: 3 Framework
text: In this paper, we address the global optimisation problem formulated as K∗ = argmaxK∈Hκf(K), where f : Hκ → R is an expensive objective functional and κ is a hyperkernel. In particular, we are interested in finding the best kernel K∗ ∈ Hκ to maximise the model performance represented by the objective functional f (for example, f can be the leave-one-out classification performance of a SVM classifier). First, we describe the construction of valid kernel functionals using hyperkernel, followed by a discussion on the kernel functional optimisation using Bayesian optimisation. A flowchart
describing the overall optimisation process of kernel functionals is shown in Figure 1. A complete algorithm for the Kernel Functional Optimisation (KFO) is given by Algorithm 1.

heading: 3.1 Construction of Kernel Functionals from Hyper-Gaussian Process
text: Ong and Smola (2003) and Ong et al. (2003, 2005) have discussed the general guidelines to design a hyperkernel. We follow the same strategy to formulate Matérn Harmonic Hyperkernel (κ):
κ(x,x′,x′′,x′′′) = 1− λh 1− (λh c1 c2 exp ( − √ 3 l (r1 + r2) ) (3) where λh and l correspond to the hyperparameters of the hyperkernel, r1 = ‖x − x′‖, r2 = ‖x′′−x′′′‖, c1 = ( 1+ √ 3 l r1 ) , and c2 = ( 1+ √ 3 l r2 ) . The derivation of Matérn Harmonic Hyperkernel is provided in the supplementary material. In our proposed method, we use the draws from a (hyper) Gaussian process GPκ(0, κ) to construct finite-dimensional subspaces of our kernel space on which we perform optimisation. As discussed in Section 2.2, the kernel samples drawn from GPκ(0, κ) do not lie inHκ, hence we approximate the draws using the posterior mean of GPκ(0, κ) lying inHκ. In practice, when sampling from GPκ(0, κ), we assume a grid G with Ng points {x̃1, x̃2, · · · |x̃i ∈ X̃ : X × X ,∀i ∈ NNg} for placing a GP distribution on kernels using a hyperkernel κ mentioned in Eq. (3). The sample set k ∼ GPκ(0, κ) is essentially a set of noiseless observations of the kernel K on the grid-points x̃1, x̃2, · · · lying inHκ′ 6=κ. The number of points in the grid is chosen such that the resulting grid is sufficiently fine to represent the kernel K everywhere on X̃ . Therefore, for any point x̃i ∈ X̃ , the posterior variance of the kernel K given the observations {(x̃i, ki) | i ∈ NNg} is negligible and thus the kernel K can be approximated using the posterior mean of GPκ(0, κ) as
K(x̃) ≈ [κ(x̃, x̃1) κ(x̃, x̃2) κ(x̃, x̃3) · · · ] κ−1 k = ∑ i αi κ(x̃, x̃i),whereα = κ−1 k (4)
A very fine resolution grid ensures that we can capture small-scale patterns in the kernel. However, a large grid size comes with large computational costs. Therefore, the choice of Ng is a trade-off between the overall computational cost and the accuracy of kernel optimisation expected. We discuss the computational complexity and the associated memory demands pertaining to Ng in Section 4.4.

heading: 3.2 Kernel Functional Optimisation
text: We adopt the ideas from Bayesian optimisation method - LINEBO (Kirschner et al., 2019) for the optimisation of non-parametric kernel functionals via a sequence of one-dimensional projections. First, we discuss the construction of low-dimensional subspaces. The key challenge here is to address the computational burden with the use of large grid. Next, we describe the Bayesian functional optimisation for each of the subspace and across many such subspaces. Since the best kernel obtained is a Kreı̆n kernel, we apply transformations to ensure the positive definiteness of the Gram matrix.
Construction of Low-dimensional Spaces We start with the construction of low-dimensional search space spanned by randomly chosen basis vectors drawn from the hyper-GP GPκ(0, κ). The hyper-GP surrogate modelling requires the computation of covariance matrix κ ∈ RNg×Ng using κ for the predefined grid G. Further, the accuracy of the kernel functional to represent the kernel K is directly proportional to the assumed grid size Ng. To avoid the computational burden arising
from the larger grid size Ng, we perform Principal Component Analysis (PCA) (Wold et al., 1987) and choose N ′ principal components. Mathematically, we represent κ = (E √ Λ)(E √ Λ)ᵀ, where ith column ei in E ∈ RNg×N ′
corresponds to the ith principal component and Λ ∈ RN ′×N ′ is the diagonal matrix containing top N ′ eigenvalues. The outer-loop in Algorithm 1 iterates through a sequence of S d-dimensional subspaces by drawing d random basis vectors in each subspace from GPκ(0, κ) i.e., k(1),k(2), · · · ,k(d) ∼ GPκ(0, κ), where k(·) = E √ Λ · β(·) and β(·) ∼ N (0, IN ′).
Kernel Optimisation Observation Model As discussed earlier, we construct kernel functionals K(·, ·) from the hyper-GP distribution GPκ(0, κ) as per Eq. (4) using
k = K# + λ(1)k(1) + · · ·+ λ(d)k(d) (5) where λ(·) ∈ [0, 1], k(·) are the random basis vectors drawn and K# corresponds to the best kernel found across all the previous subspaces. The optimal kernel in the given subspace s is obtained by optimising λ using a Bayesian optimisation procedure with another GP distribution GP(0, kSE). The observation model for GP(0, kSE) is considered as D ′
s = {(K, y = f(K))}, where K is the kernel functional constructed and y is a measure signifying the ability of the latent kernel to represent the given data. For example, log-likelihood can be used as the measure y in our observation model.
Building GP for Kernel Optimisation We fit a GP distribution GP(0, kSE) on the observed kernel functionals using the Squared Exponential (SE) kernel (kSE) given by
kSE(K1,K2) = σ 2 f exp ( −1 2Υ 2 ∥∥K1 −K2∥∥2Hκ′ 6=κ )
(6)
where σ2f and Υ correspond to the signal variance and lengthscale parameters of kSE. Although there is no restriction on the kernel choice here, we consider the commonly used SE kernel. As mentioned earlier, we approximate K using the posterior mean (µ), therefore we compute the similarity between kernel functionals using the RKHS norm (‖ · ‖Hκ ) estimated as
‖K1 −K2‖Hκ′ 6=κ ≈ ‖µ1 − µ2‖Hκ = √ αᵀ1κα1 +α ᵀ 2κα2 − 2αᵀ1κα2 (7)
where µ1 and µ2 are the posterior mean approximations of K1 and K2, respectively. We refer to the supplementary material for the details of similarity formulations using L2−Norm.
Kernel Optimisation We find the best kernel functional in the given low-dimensional subspace using GP-UCB acquisition function (Eq. (2)) with βt = 2 log(t2+ ñ 2 π2/3δ̃), where ñ corresponds to the total number of kernel functional observations and δ̃ is a value in [0, 1]. The best kernel found (K#) across all the previous subspaces acts as a subspace bias guiding the optimisation in the subsequent subspaces as per Eq. (5). The selection of S d-dimensional subspaces (outer-loop) and optimising the kernel (for T iterations) in each of the subspace (inner-loop) continues until the search budget is exhausted. The hyperparameters θ = {σ2f ,Υ} in kSE are tuned by maximising the log marginal likelihood. In addition to that, the hyperparameters of the hyperkernel (Θ = {λh, l}) mentioned in Eq. (3) are tuned using another standard Bayesian optimisation procedure. The observation model for the hyperparameter tuning of hyperkernel is constructed as D = {(Θ, y′ = Γ(Θ))}, where Γ maps the model performance y′ with the corresponding hyperparameter set Θ. We refer to the supplementary material for the detailed discussion on tuning the hyperparameters of both kernel and hyperkernel.
From Kreı̆n kernels to Positive Definite Gram Matrix
As the kernel approximated by Eq. (4) is an indefinite, or Kreı̆n kernel (K), the Gram matrix (C) constructed for the datapoints using K is also indefinite. We use the following matrix post-processing methods to ensure the positive definiteness of the Gram matrix constructed.
The Eigen Value Decomposition (EVD) based matrix post-processing involves the decomposition of the Gram matrix C as C = Z∆Zᵀ, where Z is the square matrix containing eigenvectors corresponding to the eigenvalues in the diagonal matrix ∆. The Eigen spectrum clip (∆ii = (∆ii)+) ensures positive definiteness of the given training and test covariance matrix, but in isolation, without considering the transformation of the underlying kernel function, thus resulting in inconsistency
Algorithm 1 Kernel Functional Optimisation Input: Ng - Number of points in the grid, S - Number of subspaces search, T - Number of iterations
1. Initialise (K#, ybest)← (0, 0), D0 ← ∅ 2. Compute κ for Ng grid points x̃1, x̃2,· · · using Eq. (3) 3. Perform PCA of κ as κ = (E √ Λ)(E √ Λ)ᵀ 4. for Subspace s = 1, 2, · · · , S do (outer-loop) 5. Sample k(1),k(2), · · · ,k(d) ∼ GPκ(0, κ) 6. Generate random initial observations in the current subspace s
D′s = {(K, y) |K Eq. (4)←−−−− K#+λ(1)k(1)+ · · ·+λ(d)k(d), y = f(K), λi∈Nd ∼ U(0, 1)}
7. for each iteration t = 1, 2, · · · , T do (inner-loop) 8. Solve λ∗ = argmax
λ∈[0,1]d ut(µ(K(λ)) +
√ βt σ(K(λ)))
9. Compute the new kernel Knew as Knew Eq. (4)←−−−− K# + λ(1)∗ k(1) + · · ·+ λ(d)∗ k(d)
10. Use the kernel Knew and Ĉ to measure the fitting quality y as ynew = f(Knew) 11. D′s ← D ′
s ∪ {(Knew, ynew)} 12. end for 13. Ds ← Ds−1 ∪ D ′
s
14. (K#, ybest) = argmax (K,y) ∈Ds y 15. end for 16. K∗ ← K# 17. return (K∗, ybest)
(see discussion 2.2 in Chen et al. (2009)). Therefore, to consistently transform both the training and test points, the Eigen spectrum clip is treated as a linear transformation on the training points first i.e., Ĉtrain = ϑclipCtrain, where ϑclip is the spectrum transformation matrix and then, apply the same transformation on ctest = [K(xtest,x1)K(xtest,x2) · · · ]ᵀ as ĉtest = ϑclipctest , whereϑclip = Z∆clipZᵀ and ∆clip = diag(J∆11 ≥ 0K, J∆22 ≥ 0K, · · · ). The magnitude of change in the transformed matrix (Ĉ) from the given indefinite kernel matrix (C) is minimum with the spectrum clip transformations i.e., Ĉclip = argminĈ<0 ‖C− Ĉ‖F. We note that, it is possible to use the original optimised kernel for specialised SVMs (Ying et al., 2009), but we consider this as part of the future work.
For GPs, there is a strong requirement that the covariance matrix is positive definite as it needs to generate positive definite covariances. Ayhan and Chu (2012) have demonstrated the vulnerabilities of GP with indefinite kernels. The aforestated EVD based post-processing gets complicated for GP. The GP predictive distribution involves the calculation of mean µ(·) and variance σ2(·) for the test samples. The variance requires the computation of [K(xtest,xtest)]. Although the linear transformation ϑclip on Ctrain ensures positive definiteness of ctest = [K(xtest,x1)K(xtest,x2) · · · ]ᵀ, it does not consistently transform [K(xtest,xtest)]. Therefore, we need ways to enforce positive definiteness before we compute predictive variances. To ensure positive definiteness in GPs, we clip the values of α i.e., α = [(αi)+] in the posterior mean approximation of kernels by visualising the kernel approximation (Eq. (4)) in terms of the representer theory mentioned in Ong et al. (2005).

heading: 4 Theoretical Analysis

heading: 4.1 Inner-loop
text: The cumulative regret for the optimisation in the inner-loop is given as RT = ∑T t=1 f(K
∗)− f(Kt), where K∗ is the best kernel found across all the subspaces. In the inner-loop, our goal is to derive the upper bound for the cumulative regret (RT ) in terms of the total number of iterations T .
In conventional BO algorithms, the variables being optimised are directly used in the model construction. In contrast, the inner-loop in our proposed method constructs the model using the projection of the variables (λ∗) being optimised in the functional space i.e., k = K# + ∑ i λ (i)k(i).
Proposition 1: Let Ss be the subspace constructed in each instance s of the outer-loop. Then, at each iteration t of the inner-loop, the maximum information gain (γt) of the kernel k : Ss × Ss → R is same as that of the information gain of the standard kernel in Euclidean space k : X × X → R. The proof of proposition 1 is deferred to the supplementary material.
It is important to note that the model for f in the inner-loop is constructed with the observations obtained from the current and previous subspaces search and not just the observations from the current search. Therefore, the bounds on the overall regret for the inner-loop can be derived as follows.
Theorem 1: Let f(K)|Ds−1 be the posterior of f in subspace s before entering the inner-loop and f(K)|Ds−1 ∪ D ′
s be the posterior at iteration t of the inner-loop. Then, the updated posterior f(K)|Ds−1 ∪D ′
s is equivalent to the posterior of the biased GP with prior covariance k̂Ds−1 and the inner-loop regret grows sub-linearly asO∗( √ dtγDs−1,t), where γDs−1,t is the maximum information gain for the prior covariance k̂Ds−1 andO∗ notation is a variation ofO with log factors suppressed. The proof of Theorem 1 is provided in the supplementary material.

heading: 4.2 Outer-loop
text: We provide a theoretical analysis of the outer-loop based on the notion of effective dimension (Kirschner et al., 2019, Wang et al., 2016). As we deal with the functionals in our proposed method, the standard definition of effective dimension is slightly modified as follows:
Definition 3: A function f : Hκ → R is said to have effective dimensionality d′ ∈ N, if there exists k(1),k(2), · · · ,k(d′) ∈ Hκ , such that ‖f(K + K⊥) − f(K)‖ = 0,∀K ∈ K,∀K⊥ ∈ K⊥, where K = span(k(1),k(2), · · · ,k(d′)) and K⊥ = {K̃ ∈ Hκ | 〈K, K̃〉Hκ = 0,∀K ∈ K}. Following Kirschner et al. (2019), we derive the regret bounds for the outer-loop.
Theorem 2: Given a twice Frechet-differentiable kernel k : Hκ × Hκ → R, let 0 < δ < 1, f ∼ GP(0, k) with effective dimension d′ and maxima K∗ = argmaxK∈Hκ f(K). Then, after s subspaces search (s outer-loop iterations), with probability at least 1−δ, the regret f(K∗)−f(K#) ∈ O(Jd < d′K( 1s log( 1δ )) 2 d′−d + d,δ), where K# is the best kernel found across all the previous subspace searches and d,δ is the regret bound for the inner-loop and J·K is the Iverson bracket. The proof of Theorem 2 is provided in the supplementary material.

heading: 4.3 Overall Convergence
text: In LINEBO, one-dimensional subspaces (or the lines) are optimised up to err(K+) < for some fixed (Lemma 4 of Kirschner et al. (2019)) and K+ = argmaxKi∈K1:t f(Ki). In our method, for a given subspace s, we terminate after T iterations with accuracy err(K+) ≤ d,δ. In our setup with d = 1, given a fixed budget (T iterations) for the inner-loop, we get 1,δ ∈ O(T c− 1 2 ), where c ∈ (0, 0.5) (Assumption 2 in Kirschner et al. (2019)). On the other hand, if the number of vectors (d) spanning the random basis is same as the effective dimensionality (d′), then our convergence is analogous to REMBO (Wang et al., 2016), with the regret imposed only by d′,δ . Further, the order of regret bound in such cases remains unchanged even if we consider only one subspace search (S=1).
Alternatively, simple regret measure implemented as a terminating condition in the inner-loop results in the regret bound d,δ = . If we consider one-dimensional spaces (d = 1) and use err(K+) < as the terminating condition for the inner-loop, the convergence guarantee of our algorithm is exactly same as that of LINEBO with d,δ = . Thus, the inner-loop of our algorithm is expected to complete in T ∈ O( 21−2c ) iterations for some c ∈ (0, 0.5) (see discussion around Assumption 2 in Kirschner et al. (2019)), resulting in O(S 21−2c ) total number of function evaluations overall.

heading: 4.4 Computational Analysis
text: The computational complexity of our approach is in the order of O(STN3g ), where S is the number of subspace searches, T is the number of iterations in each subspace and Ng is the number of points in the grid, without including the complexity of the downstream class (as it would be different for
different kernel machines). The main bottleneck of our method is the computation of the covariance matrix κ ∈ RNg×Ng . To avoid the computational burden resulting from the large covariance matrix κ for the given Ng , we perform Principal Component Analysis (PCA) of κ. Here, we do not perform a full PCA, rather we choose only top N ′ principal components (N ′ Ng). The computational complexity of finding top N ′ principal components is O(N ′N2g ), which is much lower than O(N3g ). Moreover, we perform PCA only once, prior to entering the outer and inner optimisation loops. Thus, we incur a cost on startup but are rewarded with significant computational savings in the main optimisation loop where the computational burden is proportional to N ′ rather than N2g . The memory complexity for optimising the kernel functionals using our proposed method is in the order ofO(N2g ). Further, as we deal with a kernel selection problem, we are only concerned with the complexity of the observed search (kernel) space. Theoretically, the optimality of our method is not limited to any dataset-specific characteristics such as the number of dimensions (n) or the number of target classes in the given problem. Such characteristics do not have a significant role in the kernel optimisation, but the complexity of the given search (kernel) space plays a vital role in the optimisation performance.

heading: 5 Experiments
text: We evaluate the performance of our proposed algorithm (KFO) on synthetic benchmark functions and also apply our method on real-world datasets for SVM classification and GP regression tasks. We have considered the following experimental settings for KFO. We have used Matérn Harmonic Hyperkernel (Eq. (3)) to define the space of kernel functionals. To express the kernel as kernel functional in Hyper-RKHS, we consider Ng & 10 × n for a given n dimensional problem. The outer-loop representing the number of low-dimensional subspace searches (S) to find the best kernel function is restricted to S = 5 and the number of iterations (T ) in each of the subspace (inner-loop) is restricted to T = 20. We use GP-UCB acquisition function to guide the search for optimum in all our experiments and at all levels. The hyperparameters λh and l of the hyperkernel (Eq. (3)) are tuned in the interval (0, 1] using a standard BO procedure mentioned in the supplementary material.

heading: 5.1 Synthetic Experiments
text: In this experiment, we test our algorithm (KFO) with the following synthetic functions: (i) Triangular wave, (ii) a mixture of three Gaussian distributions (Gmix), and (iii) SINC function. We compare with the following stationary and non-stationary kernels: (i) SE kernel, (ii) Matérn kernel with ν = 3/2 (Mat3/2), and (iii) Multi-Kernel Learning (MKL) as a linear combination of SE, Mat3/2 and Linear kernel. The hyperparameters Υ, σ2f and weights w (in the case of MKL) of the baseline kernels are tuned by maximising the log-likelihood. We compute the posterior distributions for the aforesaid synthetic functions. We report the mean and the standard deviation of the maximum log-likelihood computed over 10 random runs. We show the posterior distribution and the maximum log-likelihood estimates obtained for Triangular wave function in Figure 2. We refer to the supplementary material for the results on other synthetic functions. It is evident that the posterior distribution computed using the standard kernels has poor predictions in the held-out test region. By contrast, the kernel suggested by KFO has better predictive mean and variance in the held-out test region. Especially note that the KFO optimised kernel was able to find the correct periodicity even without explicit enforcement.

heading: 5.2 Real-world Experiments
text: We compare the performance of our proposed algorithm in SVM classification and GP regression tasks against the state-of-the-art baselines. In our classification and regression experiments, we use the publicly available multi-dimensional real-world datasets from the UCI repository (Dua and Graff, 2017). In SVM classification problems, we use C-SVM in conjunction with KFO to minimise the test classification error (Er). We perform 10-fold cross-validation on the training data set containing 80% of the total instances and tune the cost parameter (C) of the SVM in the exponent space of [−3, 3]. We compare our results with Radial Basis Function (RBF) based traditional C-SVM classifier (SVMRBF) and MKL based SVM classifier (SVM-MKL). We also compare with ν parameterised Linear SVM (ν−SVM) adhering to the definition of the hyperkernel optimisation problem using the results mentioned in Ong and Smola (2003). The classification error (in %) obtained for the test set consisting of 20% of the total instances using different classifiers over 10 random runs are shown in Table 1. To demonstrate the efficiency of our approach, we also present the best test classification error (last column of Table 1) reported by state-of-the-art classifiers in the literature (Zhang et al., 2017). To the best of our knowledge, Zhang et al. (2017) is the most recent work that surveyed numerous classifiers and reported their performance on UCI datasets. Additionally, we also construct a SVM classifier (KFO-MKL) with its kernel formulated as a weighted combination of KFO tuned kernel and standard kernels (analogous to MKL), we refer to the supplementary material for the results with KFO-MKL.
In GP regression tasks on UCI datasets, we compute the negative log-likelihood (Eq. (1)) on the test set as a measure of performance. We compare our results with the standard parametric kernels such as RBF and Automatic Relevance Determination (ARD) Matérn kernel and the non-parametric kernels such as Functional Kernel Learning based kernels (FKL-Shared and FKL-Separate) mentioned in Benton et al. (2019). In FKL-Separate, the functional kernel learning is achieved by formulating a product of one-dimensional kernels, each of which has its own GP and hyperparameters. In contrast, FKL-Shared uses a GP with unique set of hyperparameters to draw one-dimensional kernels. The results of our GP regression tasks are shown in Table 2, with each cell containing the mean negative log-likelihood and the standard deviation computed over 10 repeated runs with random 80/20 train/test splits. Evidently, our method outperformed the state-of-the-art baselines in both the SVM classification and GP regression experiments, demonstrating the significant improvement in generalisation performance. We refer to the supplementary material for the experimental details and the additional results. The code base used for the experiments mentioned above is available at https://github.com/mailtoarunkumarav/KernelFunctionalOptimisation
To provide brief insights on the computational time, we have reported the average CPU time (in %) spent optimising (or searching) the kernel and the average CPU time (in %) spent evaluating the kernel by our approach in Table 3. We observe that the percentage of time spent optimising the kernel is no more than 10% of the whole model fitting time. Thus, the proposed method does not add much overhead to the model fitting process. We have also measured the total runtime (in seconds) required for an instance of KFO tuned SVM to complete S × T iterations, where S = T = 5. The total runtime also includes the runtime required for generating 4 random observations in each subspace. The aforesaid runtimes are measured on a server with Intel Xeon processor having 16 GB of RAM.
Furthermore, we ideally expect our proposed method to at least achieve the generalisation performance demonstrated by any standard parametric kernel for the reason that we find the optimum kernel in the whole space of kernels composed of a plethora of kernels including parametric kernels. Although our proposed approach is able to find the global optimal kernel in most cases, we do occasionally observe that our method does not provide the optimal kernel. A possible reason for this could be the insufficient computational budget allocated or the substandard approximations and optimisations. Our empirical results have demonstrated that we can achieve a good generalisation performance even with smaller grids (smaller Ng) using Kernel Functional Optimisation (KFO) framework.

heading: 6 Conclusion
text: We present a novel formulation for kernel selection via the optimisation of kernel functionals using Bayesian functional optimisation. The kernel functional learnt is a non-parametric kernel capable of capturing the intricate stationary and non-stationary variations. Our algorithm iteratively searches through a sequence of random kernel functional subspaces where the best kernel obtained from all the previous subspace searches biases the next search. The resultant kernel is an indefinite, or Kreı̆n kernel, thus we use matrix post-processing techniques to ensure the positive definiteness of the resulting Gram matrix. The theoretical analysis shows a fast convergence rate of our algorithm. The experimental results show that our method outperforms the other state-of-the-art baselines.

heading: Acknowledgments
text: This research was partially funded by the Australian Government through Australian Research Council (ARC). Prof. Venkatesh is the recipient of an ARC Australian Laureate Fellowship (FL170100006).

Review summaries:
This work is interested in kernel selection. In contrast with most of the literature they consider non-parametric class of kernels using hyper-kernels. Nevertheless, some previous works have already consider this setting but the authors argue that they had several limitations which motivates a new method. The authors first recall the background on Bayesian optimization and hyper-kernels then explain their method to mix the two while avoiding some computational complexity issues. They also show how to fix a naturally arising problem : the posterior mean of the hyper-GP may be positive. Finally, some theoretical and experimental results are provided.

This article proposes to use the hyperkernel formalism coupled with Gaussian Processes to derive the ‘best’ kernel for regression. The latter kernel depends mostly on the observations rather than on a pre-selected kernel family. The proposed procedure relies on several stages of Bayesian optimization as well as a clipping of the coefficients to avoid dealing with (indefinite) Krein kernels. The performance is presented on both synthetic and benchmark datasets, to assess visually the fit of the confidence intervals as well as the numerical performance.

The paper proposes a zero-order optimization method where the optimized variable is a kernel function in Hyper-RKHS induced by a selected hyper-kernel. The method is an instance of the Bayesian Optimization method LINEBO [Kirschner 2019]. The contribution is in the adaptation of the LINEBO algorithm for efficient optimization w.r.t. positive definite kernel functions. The algorithm is applied to the optimization of kernel functions for C-SVM and GP regression. Experiments show significant improvement over existing methods. 

This paper introduces a novel approach which constructs the best fitting kernel function with adaptive tuned parameters. This kernel function is obtained by the linear combination of multiple kernels sampled from a prior Gaussian Process. The experimental results of two different kernel-based algorithms show the superiority of the proposed method.