Bandwidth Selection for Gaussian Kernel Ridge Regression via Jacobian Control

Most machine learning methods depend on the tuning of hyper-parameters. For kernel ridge regression (KRR) with the Gaussian kernel, the hyper-parameter is the bandwidth. The bandwidth specifies the length-scale of the kernel and has to be carefully selected in order to obtain a model with good generalization. The default method for bandwidth selection is cross-validation, which often yields good results, albeit at high computational costs. Furthermore, the estimates provided by cross-validation tend to have very high variance, especially when training data are scarce. Inspired by Jacobian regularization, we formulate how the derivatives of the functions inferred by KRR with the Gaussian kernel depend on the kernel bandwidth. We then use this expression to propose a closed-form, computationally feather-light, bandwidth selection method based on controlling the Jacobian. In addition, the Jacobian expression illuminates how the bandwidth selection is a trade-off between the smoothness of the inferred function, and the conditioning of the training data kernel matrix. We show on real and synthetic data that compared to cross-validation, our method is considerably more stable in terms of bandwidth selection, and, for small data sets, provides better predictions.

However, kernelization introduces hyper-parameters, which need to be carefully tuned in order to obtain good generalization of the inferred model. Bandwidth, σ, is a hyper-parameter used by many kernels, including the Gaussian, or radial basis function, kernel. The bandwidth specifies the length-scale of the kernel. A kernel with a too small bandwidth will treat most new data as far from any training observation, while a kernel with a too large bandwidth will treat each new data point as basically equidistant to all training observations. None of these situations will result in good generalization.
The problem of bandwidth selection has been extensively studied for kernel density estimation, KDE, which is the basis for KDE-based kernel regression, such as the Nadaraya-Watson estimator (Nadaraya, 1964;Watson, 1964) and locally weighted regression (Cleveland and Devlin, 1988). Köhler et al. (2014) review existing methods for bandwidth selection for KDE-based kernel regression, methods that all make varyingly strong assumptions on the underlying data density and smoothness of the non-parametric regression model, and how the latter can be approximately estimated. On one end of the spectrum, cross-validation makes almost no assumptions on the underlying data structure, resulting in a very flexible, but computationally heavy, estimator, with high variance. On the other end, with strong assumptions on the underlying data structure, Silverman's rule of thumb, originally from 1986, (Silverman, 2018) is a computationally light estimator with low variance, but possibly far from optimal. Other approaches on the spectrum include Park and Marron (1990); Sheather and Jones (1991); Fan and Gijbels (1995).
Although similar in name and usage, KRR and KDE-based kernel regression are not the same. While KDE-based kernel regression estimates the probability density of the data, and use this density to estimate E[y|x], KRR takes a functional perspective, directly estimatingŷ =f (x), similarly to how is done in neural networks.
For neural networks, Jacobian regularization, which penalizes the Frobenius norm of the Jacobian, ∂f (x) ∂x 2 F , has recently been successfully applied to improve generalization (Jakubovitz and Giryes, 2018;Chan et al., 2019;Hoffman et al., 2019;Finlay et al., 2020;Bai et al., 2021). The Jacobian penalty is a non-linear generalization of the linear ridge penalty. To see this, consider the linear modelf (x) = x β, for which the Jacobian penalty becomes exactly the ridge penalty, β 2 2 . Thus, both Jacobian and ridge regularization improve generalization by constraining the derivatives of the inferred function.
This connection motivates our investigation into how Jacobian constraints can be applied for bandwidth selection in KRR: If we knew how the kernel bandwidth affects the Jacobian of the inferred function, then we could use Jacobian control as a criterion for selecting the bandwidth.
Our main contributions are: • We derive an approximate expression for the Jacobian of functions inferred by KRR with the Gaussian kernel.
• We propose a closed-form, computationally feather-light, bandwidth selection method for KRR with the Gaussian kernel, based on controlling the approximate Jacobian. • We show on synthetic and real data that Jacobian based bandwidth selection outperforms cross-validation in terms of stability and, when data are scarce, performance.

Bandwidth Selection through Jacobian Control
Consider the left panel of Figure 1. When large (absolute) derivatives of the inferred function are allowed, the function varies more rapidly between observations, while a function with constrained derivatives varies more smoothly, which intuitively improves generalization. However, the derivatives must not be too small as this leads to an overly smooth estimate, as seen in the right panel.
The functions in Figure 1 are all constructed using kernel ridge regression, KRR, with the Gaussian kernel. For training data X ∈ R n×p and y ∈ R n , the objective function of KRR is H k denotes the reproducing kernel Hilbert space corresponding to the symmetric, positive semi-definite kernel function k(x, x ), and λ ≥ 0 is the regularization strength. Solving Equation 1, a predictionf (x * ) ∈ R, where x * ∈ R p , is given byf where K(x * , X) ∈ R 1×n and K(X, X) ∈ R n×n are two kernel matrices; K(X, X ) ij = k(x i , x j ). In the absence of regularization, regardless of the bandwidth, the inferred function perfectly interpolates the training data, i.e. it hits all training observations. When the bandwidth is too small, the kernel considers most new observations as far away from any training data, and quickly resorts to its default value, 0. A too large bandwidth on the other hand, results in extreme predictions between some of the observations. The addition of regularization affects larger bandwidths more than smaller. A too large bandwidth, in combination with regularization, produces a function too simple to capture the patterns in the training data.
The Gaussian kernel is given by where the bandwidth, σ, specifies the length scale of the kernel, i.e. what is to be considered as "close".
Returning to our question above, how can we select σ in Equation 3 to control ∂f (x * ) , withf (x * ) given by Equation 2? In general, there is no simple expression for ∂f (x * ) ∂x * 2 , but in Definition 1, we state an approximation based on derivations presented in Section 2.1.
where l max denotes is the maximum distance between two training observations, and C(n, y 2 ) is a constant with respect to σ.
Remark 1: Since we are only interested in how J a 2 depends on σ, we will henceforth, with a slight abuse of notation, omit the constant C(n, y 2 ).
Remark 2: Technically, since we use a univariate response, ∂f (x * ) ∂x * is a special case of the Jacobian, the gradient. We chose however to use the word Jacobian, since nothing in our derivations restricts us to the univariate case. Proposition 1 below characterizes how the approximate Jacobian norm, J a 2 , depends on σ. We see that, depending on λ, it can behave in three different ways: In the absence of regularization, J a 2 becomes arbitrarily large for σ small or large enough, and enjoys a global minimum at a value σ 0 , which is consistent with the left panel of Figure 1. However as soon as regularization is added, J a 2 goes to 0 as bandwidth goes to infinity, as indicated in the right panel. As long as the regularization parameter λ ≤ 2ne −3/2 ≈ 0.45n, there is still a local minimum at σ 0 . This is further illustrated in Figure 2, where we plot J a 2 together with its components j a and j b for three different values of λ, reflecting the three types of behaviour. Since j b is bounded by 1/λ, for λ large enough it is negligible compared to j a . Proposition 1. Let J a 2 (σ) be defined according to Definition 1, and let, for k ∈ {−1, 0}, where W k denotes the k-th branch of the Lambert W function. Then Sketch of Proof. We evaluate J a 2 (0), J a 2 (∞) and points where ∂J a 2 (σ) ∂σ = 0 to obtain extreme point candidates. To avoid infeasible solutions, we have to consider the domain of the Lambert W function. For the full proof, see Appendix A.
Based on Proposition 1 we can now propose a bandwidth selection scheme based on Jacobian control: For λ ≤ 2ne −3/2 , we choose σ 0 as our Jacobian based bandwidth. For λ > 2ne −3/2 , σ 0 is not defined; in this case we choose our bandwidth as if λ = 2ne −3/2 , the largest allowed value of λ. Note σ 0 is quite stable to changes in λ. The square root expression in Equation 5 increases from 1 for λ = 0 to √ 3 for λ = 2ne −3/2 .

Theoretical Details
We first use Proposition 2 to approximate the norm of the two kernel matrices in Equation 2 with a product of two matrix norms. We then use Propositions 3 and 4 to estimate these two norms for the case of the Gaussian kernel. Note that Proposition 2 holds for any kernel, not only the Gaussian. Proposition 2.
Sketch of Proof. We first use submultiplicativity to split the norm it into three factors, and equivalence of matrix norms to replace the 2-norm with the 1-norm in the first factor. Finally, we show that ∂f (x * ) ∂x * = ∂f (x * ) ∂d i to obtain the forms for the two factors. For the full proof, see Appendix A. , and its two factors (right y-axis) as functions of the bandwidth for three different values of λ, and n = 10, l max = 1, p = 1. Note that the scales of the axes differ between the three panels and that j a (σ) is scaled by a factor 100. For Sketch of Proof. Since the Gaussian kernel is rotationally invariant we only need to consider the radial coordinate, d i , when calculating the gradient. The value of d i that maximizes the gradient is calculated by setting the derivative (of the gradient) to zero. For the full proof, see Appendix A.
Sketch of Proof. Bermanis et al. (2013) provide an estimate for the number of singular values larger than δ · s 1 , for a Gaussian kernel matrix, where δ > 0 and s 1 denotes the largest singular value. Using this expression we can upper bound the smallest singular value, or equivalently, lower bound the largest singular value of the inverse matrix. The addition of λI shifts all singular values by λ. For the full proof, see Appendix A.
In the absence of regularization, j b (σ) is a bound of the spectral norm of the inverse training data kernel matrix. With increasing σ, the elements in K(X, X, σ) become increasingly similar, and K(X, X, σ) becomes closer to singular, which results in an ill-conditioned solution, wheref (x * ) is very sensitive to perturbations in X. Introducing regularization controls the conditioning, as seen in Figure  With the Jacobian approach, the contribution of an ill-conditioned of matrix, j b (σ), is balanced by how quickly the inferred function decays in the absence of training data, j a (σ). Returning to the left panel of Figure 1, for the small bandwidth, the inferred function quickly decays to 0, j a (σ) is large; while for the large bandwidth, K(X, X, σ) is almost singular, j b (σ) is large, resulting both in large derivatives and in extreme predictions between observations.
• 1D Temperature Data: The temperature at 248 different times in January 2020 at the weather station at Toulouse-Blagnac The French temperature data was obtained from Météo France 1 and processed following the setup by Vanwynsberghe (2021).
In all experiments the Jacobian based bandwidth was compared to 10-fold cross-validation with 100 logarithmically spaced values between 0.01 and l max . We also included Silverman's rule of thumb, σ = 4 n(p+2) 1/(p+4) ·σ, whereσ is the estimated standard deviation of the observations. Note however that Silverman's method is developed with KDE, rather than KRR, in mind and thus does not take λ into account. We chose however to include it as a reference since, just like the Jacobian method, it is a computationally light, closed-form solution. To avoid singular matrices, a small regularization of λ = 10 −3 was used, unless otherwise is stated. Each experiment required less than one minute on a laptop computer, but since we used multiple values of n, λ and random seeds, we used a cluster with Intel Xeon Gold 6130 CPUs to run experiments in parallel. Figures 3 and 4 show the results of jackknife resampling for the two temperature data sets. For each data set, the experiments were repeated n times, omitting a different observations in each experiment. For the 1D temperature data, half of the observations were set aside as reference data. Apart from plotting the jackknife mean and standard deviations of the prediction, we also state the mean and standard deviations of the selected bandwidths. In the 1D case, the Jacobian method and cross-validation perform very similarly. In the 2D case, when data are more scarce, the Jacobian method is more stable than cross-validation.
In Figure 5 we alternately vary n for fixed λ, and λ for fixed n, for the three data sets. For the synthetic data 1000 test observation were generated, while the real data was randomly split into training and testing data. When varying n, 15 % of the real data was saved for testing, resulting in splits of size n/1000, n/37 and n/6 for the synthetic, 1D temperature and 2D temperature data, respectively. When varying λ, the training/testing splits were 40/1000, 100/148 and 25/15. In both cases 1000 random splits were used to estimate the variance of R 2 , i.e. the proportion of the variation inf (x * ) that is explained by x * , and the selected bandwidth σ. It is again confirmed that the Jacobian method is much more stable in terms of bandwidth selection. We also note that, for small sample sizes, the Jacobian method does better than cross-validation in terms of R 2 . This becomes increasingly clear when we consider the worst case scenario (lower bound). In the bottom panel of Figure 5, we can see where λ > 2ne −3/2 , since this is where the Jacobian selected bandwidth stops increasing with λ. We note that this is close to where R 2 = 0, i.e. the regularization is too strong to provide a useful model.
The superior stability of the Jacobian method opens up for combining it with cross-validation, using the former as a starting-point for the latter, thus reducing the variance and/or computational cost of cross-validation. We include some promising preliminary results of such a scheme in Appendix B . However, this should be investigated further both in terms of computational trade-off and performance, but this is left for future work.

Conclusions
We proposed a computationally low-cost method for choosing the bandwidth for kernel ridge regression with the Gaussian kernel. The method was motivated by the observed improved generalization properties of neural networks trained with Jacobian regularization. By selecting a bandwidth with Jacobian control, we implicitly constrain the derivatives of the inferred function. To achieve this, we derived a simple, closed-form approximate expression for the Jacobian of Gaussian KRR as a function of the bandwidth and were thus able to find an optimal bandwidth for Jacobian control. We demonstrated how selecting the optimal bandwidth is a trade-off between utilizing a well-conditioned training data kernel matrix and a slow decay of the inferred function towards the default value. Compared to crossvalidation, apart from being much simpler to evaluate, our method is much more stable and exhibits better performance for small data. Even though we only considered Jacobian bandwidth selection for the Gaussian kernel, the principle holds for any kernel. That however, requires new, kernel-specific, estimates of the kernel matrix norms. Similarly, in the Gaussian case the estimate of the norm of the inverse kernel matrix could potentially be further improved. Finally, preliminary results indicated that Jacobian bandwidth selection could be used to guide cross-validation, reducing computational burden and improve on bandwidth selection stability. These research problems are left for future work.
Code is available at https://github.com/allerbo/jacobian_bandwidth_selection.  Cross-validation Silverman Figure 5: Mean and 5-th and 95-th percentiles of R 2 of the predictions, and selected bandwidth, σ, as functions of training sample size, n (top two panels), or regularization, λ (bottom two panels). The Jacobian method outperforms cross-validation in terms of prediction when n is small, especially in terms of lower bound. It is also the most stable method in terms of bandwidth selection.
where W denotes the Lambert W function. Since this function has real outputs only if its argument is greater than −e −1 , in order to obtain stationary points we need ≥ −e −1 ⇐⇒ λ ≤ 2ne −3/2 which gives us the two stationary points 0), which implies that σ 0 < σ −1 . Combined with the limits above, this implies that, when existing, σ 0 is a local minimum and σ −1 is a local maximum.

Proof of Proposition 2.
We first note that for where in (a) we used the chain rule together with ∂d i ∂x * = I.

Proof of Proposition 3.
In spherical coordinates, where the sum is over the angular coordinates. Since the Gaussian kernel is rotationally invariant, this sum is 0 and To find the d i that maximizes the derivative, we look where the second derivative is zero.
Setting the second derivative to zero amounts to Plugging this into the first derivative we obtain 1 σ exp − 1 2 which is greater than Proof of Proposition 4.
To alleviate notation, from now on we do not explicitly state that K depends on X. We first note that (K + λI) −1 2 = 1 sn(K+λI) , where s n denotes the smallest singular value. Since K is symmetric and positive semi-definite, it is diagonalizable as K = U ΣU , while λI = λU U , which means that K + λI = U (Σ + λI) U , i.e. the singular values of K + λI are the singular values of K, shifted by λ.
According to Bermanis et al. (2013), for x ∈ R p , where each x i is restricted to an interval of length l i , i = 1, 2, . . . p, for a Gaussian kernel matrix K ∈ R m×n , with singular values s 1 , . . . s n , the number of singular values larger than δ · s 1 for some δ > 0, R δ (K), is bounded according to Solving for δ, we obtain Now, if R δ (K) = n, then all singular values (including s n ) are larger than or equal to δ · s 1 . If R δ (K) = n − 1, then all but one (namely s n ) of the singular values are larger than or equal to δ · s 1 . So for R δ (K) = n − 1, s n < δ · s 1 , which implies where we used σ 1 (K)· ≤ n · K max = n · 1. Thus

B Preliminary Results for Jacobian Seeded Cross-validation
For our Jacobian seeded cross-validation we used 10-fold cross-validation with 100 logarithmically spaced values between 1 5 σ 0 and 5σ 0 . In Figure 6, we add the results of this method to those already presented in Figure 5. Compared to standard cross-validation, the bandwidth is much more stable, and predicted R 2 gets slightly better, especially for small data sets.

C More on Kernel Ridge Regression
Below follows a slightly more detailed introduction to kernel ridge regression than in Section 2

C.1 Linear Ridge Regression
In linear ridge regression, predictions,ŷ * ∈ R n * , are made from new data, X * ∈ R n * ×p , according tô whereβ ∈ R p is the solution to arg min β y − Xβ 2 2 + λ β 2 2 , with X ∈ R n×p and y ∈ R n denoting training data, and where λ ≥ 0 is the strength of the regularization. Including the ridge penalty, λ β 2 2 , reduces the complexity of the inferred model by constraining β 2 . A larger value of λ results in a smaller value of β 2 .
The analytical solution toβ, is given bŷ Cross-validation Silverman Jacobian Seeded Cross-validation Figure 6: Mean and 5-th and 95-th percentiles of predicted R 2 and selected bandwidth, σ, as functions of training sample size, n, or regularization, λ. Compared to standard cross-validation, Jacobian seeded cross-validation is more stable and performs slightly better when data is scarce.
where the right hand side is sometimes called the dual form. The fact that can be proven using the singular value decomposition of X. Unless X is invertible, for λ = 0 only one of the inverses in Equation 12 exists, and the existing one can be used to define the Moore-Penrose pseudo-inverse of X: When X is invertible so are both matrices above; in this case the two definitions of X + coincide and X + = X −1 .
For λ = 0, we can combine Equations 9 and 11 with the definition of the pseudo-inverse to obtain β = X + ŷ y * = X * X + y.

C.2 Kernel Ridge Regression
Kernel ridge regression is a non-linear generalization of Equation 10 into where f (X) := [f (x 1 ) . . . f (x n )] and H k is the reproducing kernel Hilbert space corresponding to the symmetric, positive semi-definite kernel function k(x, x ). For kernel ridge regression, new predictions are modelled as functions of new observations according tô y * = f (X * ) = K (X * , X) · (K (X, X) + λI) −1 · y, where K(X * , X) ∈ R n * ×n ∈ R and K(X, X) ∈ R n×n are two kernel matrices, K(X, X ) ij = k(x i , x j ). We note that for kernel ridge regression with the linear kernel, k(x i , x j ) = x i x j , we obtain K(X, X ) = XX and kernel ridge regression becomes exactly linear ridge regression, on the dual form.

C.3 Kernel Ridge Regression as "Feature Space Ridge Regression"
Every kernel can be written as the inner product of a feature expansion of x, k(x, x ) = ϕ(x) ϕ(x ) =⇒ K(X, X ) = Φ(X)Φ(X ) , where ϕ(x) ∈ R p1 and Φ(X) ∈ R n×p1 . In general p 1 p and often the feature space is infinite-dimensional. Using the feature expansion of the kernel, we can think of kernel regression as "feature expansion regression", which is simply linear regression in feature space: For λ = 0 we obtain, similarly to in the linear case, β = Φ(X) + ŷ y * = Φ(X * )Φ(X) + y.

This means that we can express Equation 2 aŝ
f (x * ) = ϕ(x * ) Φ(X) + y, which implies This formulation makes it perhaps even clearer how the Jacobian can be split into two parts. The first part is related to how quickly the feature map changes on new data. The second part is related to how close to singular the feature expansion matrix of the training data is. However, since the feature expansion of the Gaussian kernel is infinite-dimensional, this is not a form we can easily use in practice.