text:
# Absorbing boundary conditions for simulation of gravitational waves with spectral methods in spherical coordinates

## 1 Introduction

### 1.1 Wave equations in General Relativity

The determination of numerical solutions of the Einstein equations is the scope of numerical relativity. It is a fundamental issue not only for the determination of gravitational wave signals for detector data analysis, but also for the study of the properties of relativistic astrophysical objects . Within numerical relativity studies, the most commonly used formulation of the Einstein equations is the so-called “3+1” formalism (also called Cauchy formalism ) in which space-time is foliated by a family of space-like hypersurfaces $`\mathrm{\Sigma }_t`$, which are described by their 3-metric $`\gamma _{ij}`$. The 4-metric $`g_{\mu \nu }`$ is then described in terms of $`\gamma _{ij}`$, a 3-vector $`N^i`$ (called shift) and a scalar $`N`$ (called lapse). In this formalism, the Einstein equations can be decomposed into a set of four constraint equations and six second-order dynamical equations. Solving the Einstein equations then turns to be a Cauchy problem of evolution under constraints and there remains the freedom to choose the time coordinate (slicing) and the spatial gauge.

For example, the choice of maximal slicing for the time coordinate (see ) converts the constraint equations to scalar form and a vectorial Poisson-like equation, for which a numerical method for solution has been presented in . As far as evolution equations are concerned, they consist of six non-linear scalar wave equations in curved space-time, with the additional choice of the Dirac gauge . The whole system is a mixed initial value-boundary problem, and this paper deals with boundary conditions for the time evolution equations. Indeed, a simpler problem is considered: the initial value-boundary problem for a linear and flat scalar wave equation:

$$\mathrm{}\varphi (t,r,\theta ,\phi )=\sigma (t,r,\theta ,\phi )$$
(1)
where

$$\mathrm{}\varphi =\frac{^2\varphi }{t^2}\frac{^2\varphi }{r^2}\frac{2}{r}\frac{\varphi }{r}\frac{1}{r^2}\left(\frac{^2\varphi }{\theta ^2}+\frac{1}{\mathrm{tan}\theta }\frac{\varphi }{\theta }+\frac{1}{\mathrm{sin}^2\theta }\frac{^2\varphi }{\phi ^2}\right)$$

is the usual flat scalar d’Alembert operator in spherical coordinates $`(r,\theta ,\phi )`$ and $`\sigma `$ is a source. To solve a more general problem in curved space-time, like for example:

$$\frac{^2\varphi }{t^2}\mu ^2(t,r)\mathrm{\Delta }\varphi =\varphi ^2,$$
(2)
one can put non-linear terms to the source $`\sigma `$ and represent at each time-step the metric function $`\mu ^2`$ by a polynomial (semi-implicit scheme, see  for an example in spherical symmetry).

### 1.2 Motivations for general quadrupolar absorbing boundary conditions

The study of the simple wave equation and its properties concerning quadrupolar waves is more than a toy-model for numerical relativity. There are many degrees of freedom in the formulation of the Einstein equations and in the gauge choice. It is not clear which of these formulations are well-posed or numerically stable . It is therefore important to have numerical tools that are general in the sense that they can be used within the framework of various formulations and gauges. Still, in many cases, the dynamical degrees of freedom of the gravitational field can be described by wave-like propagation equations in curved space-time. On the other hand, since we are mainly interested in the gravitational wave signal, which has a quadrupolar dominant term, we have to make high precision numerical models (including boundary conditions) to study this mode, as well as lower multipoles.

These statements can be illustrated as follows. One of the main sources we want to study are binaries of two compact objects (neutron star or black hole) orbiting around each other. Gravitational waves take away angular momentum and the system coalesces. In some perturbative approach, the terms corresponding to this “braking force” result from a subtle cancellation between terms of much higher amplitude . In numerical non-perturbative studies, the same phenomenon may happen and, if the dominant modes of the wave are not computed with enough precision, the angular momentum loss may be strongly overestimated. Moreover, the time-scale for coalescence is much larger than the orbital period and the system is almost stationary.

There has been many interesting developments concerning absorbing boundaries in the last years, with the Perfectly Matched Layers (PML, see  and ) which consist in surrounding the true domain of interest by an absorbing layer where the wave is damped. These methods may not be the best-suited for our problems since, as stated above, we might have to change the formulation of the equations we want to solve. Moreover, the main problem we want to address is the simulation of quadrupolar waves and, as it will be shown later in this paper, with our formulation it is possible to have a clear control on the behavior of these quadrupolar waves. Finally, this formulation is straightforward to implement and very little CPU time consuming in the context of spectral methods and spherical coordinates, which we are already using to solve elliptic partial differential equations (PDE) arising in numerical relativity (scalar and vectorial ones, see ). The development and implementation of the PML techniques for our problem would require much more work and computing time, whereas it is not guaranteed at all it would give better results. For all these reasons we chose to develop a new formulation of the Bayliss and Turkel  boundary conditions, particularly well suited for using with spectral methods and spherical coordinates.

The paper deals with this new formulation as well as numerical tests. It is organized as follows. First, Sec. 2 presents boundary conditions: it briefly recalls main results from Bayliss and Turkel (2.1) and we then derive the formulation adapted up to quadrupolar modes of the wave (2.2). Then, Sec. 3 briefly describes spectral methods in spherical coordinates that were used (3.1) and details the numerical results (3.2). Finally, Sec. 4 gives a summary and some concluding remarks.

## 2 Absorbing boundary conditions

An important difference between the solution of the wave equation and that of the Poisson equation (as in ) is the fact that boundary conditions cannot be imposed at infinity, since one cannot use “compactification”, i.e. a change of variable of the type $`u=1/r`$. This type of compactification is not compatible with an hyperbolic PDE, see . One has to construct an artificial boundary and impose conditions on this surface to simulate an infinite domain. These conditions should therefore give no reflection of the wave, that could spuriously act on the evolution of the system studied inside the numerical grid. The boundary conditions have to absorb all the waves that are coming to the outer limit of the grid. The general condition of radiation is derived e.g. in , and defined as

$$\underset{r\mathrm{}}{lim}\left(\frac{}{r}+\frac{}{t}\right)(r\varphi )=0.$$
(3)
At a finite distance $`r=R`$ the condition, which is then approximate, reads

$$\frac{\varphi }{t}+\frac{\varphi }{R}+\frac{\varphi }{r}|_{r=R}=0,$$
(4)
which will be hereafter referred as the “Sommerfeld condition” and is exact only for pure monopolar waves. A completely general and exact boundary condition for the wave equation on an artificial spherical boundary has recently been derived by Aladl et al.  and involves an infinite series of inverse Fourier transforms of the solution. This condition may not be suitable for direct numerical implementation for which Aladl et al. derived a truncated approximate condition.

### 2.1 Asymptotic expansion in terms of multipolar momenta

A rather general method to impose non-reflecting boundary conditions is to construct a sequence of boundary conditions that, for each new term, are in some sense giving better results. Some of the possibilities to define “better” are when the reflected wave decreases:

* as the incident wave approaches in a direction closer to some preferred direction(s) (see e.g. ),
* for shorter wavelengths,
* as the position of artificial boundary goes to infinity.

This last approach is the most relevant to the problem of solving the Einstein equation for isolated systems. It is also a way of expanding condition (3) in terms of asymptotic series, which has been studied in , where a sequence of recursive boundary conditions is derived. Let us recall here some of their results.

A radiating solution of (1) with the source $`\sigma =0`$ can be written as the following expansion:

$$\varphi (t,r,\theta ,\phi )=\underset{k=1}{\overset{\mathrm{}}{}}\frac{f_k(tr,\theta ,\phi )}{r^k}.$$
(5)
The operators acting on a function $`f(t,r,\theta ,\phi )`$ are recursively defined by:

$`B_1f`$ $`=`$ $`{\displaystyle \frac{f}{t}}+{\displaystyle \frac{f}{r}}+{\displaystyle \frac{f}{r}},`$ (6)
$`B_{n+1}f`$ $`=`$ $`\left({\displaystyle \frac{}{t}}+{\displaystyle \frac{}{r}}+{\displaystyle \frac{2n+1}{r}}\right)B_nf.`$ (7)
The family of boundary conditions then reads:

$$B_n\varphi |_{r=R}=0.$$
(8)
In , it is shown that, following from (5), a radiating solution of the wave equation verifies:

$$B_n\varphi =O\left(\frac{1}{r^{2n+1}}\right),$$
(9)
which in particular means that condition (8) is an asymptotic one in powers of $`1/r`$. The condition $`B_1\varphi =0`$ is same as the Sommerfeld condition (4) and the same as the first approximation in terms of the angle between the direction of propagation of the wave and the normal to the boundary, derived in .

Finally, using expression (5) one can verify that the operator $`B_n`$ annihilates the first $`n`$ terms of the expansion. Thinking in terms of spherical harmonics, this means that condition (8) is exact if the wave carries only terms with $`ln1`$. In other words, the reflection coefficients for all modes lower than $`n`$ are zero. Since we are interested in the study of gravitational wave emission by isolated systems, it is of great importance to have a very accurate description of the quadrupolar part of the waves, which is dominant. Therefore, if the $`l=2`$ part of the gravitational wave is well described, higher-order terms may not play such an important role in the dynamical evolution of the system. The situation then is not so bad even if only an approximate boundary condition is imposed for those terms with $`l3`$. Moreover, the error on the function scales like $`1/R^{n+1}`$ so, if we impose

$$B_3\varphi |_{r=R}=0,$$
(10)
we have an exact boundary condition for the main contribution to the gravitational wave and an error going to zero as $`O(1/R^4)`$. When developing this expression, one gets:

$`({\displaystyle \frac{^3}{t^3}}+3{\displaystyle \frac{^3}{t^2r}}+9{\displaystyle \frac{1}{r}}{\displaystyle \frac{^2}{t^2}}+3{\displaystyle \frac{^3}{tr^2}}+18{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{t}}+18{\displaystyle \frac{1}{r}}{\displaystyle \frac{^2}{tr}}+`$
$`{\displaystyle \frac{^3}{r^3}}+9{\displaystyle \frac{1}{r}}{\displaystyle \frac{^2}{r^2}}+18{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}+6{\displaystyle \frac{1}{r^3}}\left)\varphi \right|_{r=R}=0.`$ (11)
### 2.2 New formulation for quadrupolar terms

Starting from (11) and considering that $`\varphi `$ is a solution of the wave equation (1), we replace second radial derivatives with:

$$\frac{^2\varphi }{r^2}=\frac{^2\varphi }{t^2}\frac{2}{r}\frac{\varphi }{r}\frac{1}{r^2}\mathrm{\Delta }_{\text{ang}}\varphi ,$$
(12)
where:

$$\mathrm{\Delta }_{\text{ang}}\varphi =\frac{^2\varphi }{\theta ^2}+\frac{1}{\mathrm{tan}\theta }\frac{\varphi }{\theta }+\frac{1}{\mathrm{sin}^2\theta }\frac{^2\varphi }{\phi ^2}$$
(13)
is the angular part of the Laplace operator. We are making here the assumption that, at the outer boundary of the grid ($`r=R`$), the source term $`\sigma `$ of (1) is negligible. This is a very good approximation for our studies of isolated systems and is also the assumption made when writing a solution to the wave equation in the form (5). For example, the third order radial derivative is replaced with

$$\frac{^3\varphi }{r^3}=\frac{^3\varphi }{t^2r}+\frac{2}{r^3}\mathrm{\Delta }_{\text{ang}}\varphi \frac{1}{r^2}\mathrm{\Delta }_{\text{ang}}\frac{\varphi }{r}+\frac{2}{r^2}\frac{\varphi }{r}\frac{2}{r}\frac{^2\varphi }{r^2};$$
(14)
and the second-order radial derivatives of the last term (combined with its counterpart term in (11) ) is replaced once more using (12). The boundary condition is then written as:

$`B_3\varphi =(4{\displaystyle \frac{^3}{t^3}}+4{\displaystyle \frac{^3}{t^2r}}+16{\displaystyle \frac{1}{r}}{\displaystyle \frac{^2}{t^2}}+18{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{t}}+12{\displaystyle \frac{1}{r}}{\displaystyle \frac{^2}{tr}}+6{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}+`$
$`{\displaystyle \frac{3}{r^2}}\mathrm{\Delta }_{\text{ang}}{\displaystyle \frac{}{t}}{\displaystyle \frac{1}{r^2}}\mathrm{\Delta }_{\text{ang}}{\displaystyle \frac{}{r}}{\displaystyle \frac{5}{r^3}}\mathrm{\Delta }_{\text{ang}}+{\displaystyle \frac{6}{r^3}})\varphi .`$ (15)
We use the auxiliary function $`\xi `$:

$$(t,\theta ,\phi ),B_1\varphi |_{r=R}=\left(\frac{}{t}+\frac{}{r}+\frac{1}{r}\right)\varphi (t,r,\theta ,\phi )|_{r=R}=\xi (t,\theta ,\phi ),$$
(16)
which is defined on the sphere at $`r=R`$. Inserting this definition into the boundary condition $`B_3\varphi |_{r=R}=0`$, with Eq. (15), one gets:

$$\frac{^2\xi }{t^2}\frac{3}{4R^2}\mathrm{\Delta }_{\text{ang}}\xi +\frac{3}{R}\frac{\xi }{t}+\frac{3\xi }{2R^2}=\frac{1}{2R^2}\mathrm{\Delta }_{\text{ang}}\left(\frac{\varphi }{R}\frac{\varphi }{r}|_{r=R}\right);$$
(17)
which is a wave-like equation on the outer boundary of the grid, with some source term, equal to zero if the solution $`\varphi `$ is spherically symmetric. The boundary condition (10) is now equivalent to the system (16)-(17). Written in this way, this formulation can be regarded as a perturbation of the Sommerfeld boundary condition ($`B_1\varphi =0`$) given by (16). The main advantages are that it can be very easily implemented using spectral methods and spherical coordinates (see Sec. 3.1) and that mixed derivatives have almost disappeared: there is only one remaining as a source of (17).

## 3 Numerical experiments

### 3.1 Implementation using spectral methods and spherical coordinates

Spectral methods (, , for a review see ) are a very powerful approach for the solution of a PDE and, in particular, they are able to represent functions and their spatial derivatives with very high accuracy. As presented in , we decompose scalar fields on spherical harmonics $`Y_l^m(\theta ,\phi )`$, for the angular part:

$$\varphi (t,r,\theta ,\phi )=\underset{l=0}{\overset{L}{}}\underset{m=l}{\overset{l}{}}\varphi _{lm}(t,r)Y_l^m(\theta ,\phi ),$$
(18)
and on even Chebyshev polynomials $`\left(T_{2k}(x=r/R)\right)`$ for the radial part of each $`\varphi _{lm}(t,r)`$. Time derivatives are evaluated using finite-difference methods. Since Chebyshev collocation points are spaced by a distance of order $`1/N^2`$, (where $`N`$ is the highest degree of the Chebyshev polynomials used for the radial decomposition) near grid boundaries, the Courant condition on the time step for explicit integration schemes of the wave equation (1) also varies like $`1/N^2`$. This condition is very restrictive and it is therefore necessary to use an implicit scheme. We use the Crank-Nicholson scheme, which is unconditionally stable, as shown by various authors (see e.g. ). This scheme is second-order in time and the smoothing of the solution due to implicit time-stepping remains lower than the other errors discussed hereafter. This implicit scheme results in a boundary-value problem for $`\varphi `$ at each time-step. The solution to this problem is obtained by inverting the resulting spatial operator acting on $`\varphi `$ using the tau method. Its matrix (in Chebyshev coefficient space) has a condition number that is rapidly increasing with $`N`$. This can be alleviated by the use of preconditioning matrices, obtained from finite-differences operators (see ).

At the beginning of time integration, we suppose that $`\varphi `$ satisfies the Sommerfeld boundary condition (4), that is $`(\theta ,\phi )\xi (t=0,\theta ,\phi )=0`$. $`\xi `$ is then calculated at next time-step using (17). This is done very easily since the angular parts of $`\varphi `$ and $`\xi `$ are decomposed on the basis of spherical harmonics; each component $`\xi _{lm}(t)`$ is the solution of a simple ODE in time, which is integrated using the same Crank-Nicholson scheme as for the main wave equation (1), with boundary conditions such that $`\xi `$ is periodic on the sphere. This is already verified by the $`Y_l^m`$ (Galerkin method). We get, with $`\delta t`$ being the time-step, $`\varphi _{lm}^J(r)=\varphi _{lm}(t+J\delta t,r)`$ and $`\xi _{lm}^J=\xi _{lm}(J\delta t)`$:

$`{\displaystyle \frac{\xi _{lm}^{J+1}2\xi _{lm}^J+\xi _{lm}^{J1}}{\delta t^2}}`$ $`+`$ $`{\displaystyle \frac{3}{8}}{\displaystyle \frac{l(l+1)}{R^2}}\left(\xi _{lm}^{J+1}+\xi _{lm}^{J1}\right)+{\displaystyle \frac{3}{R}}{\displaystyle \frac{\xi _{lm}^{J+1}\xi _{lm}^{J1}}{2\delta t}}`$
$`+`$ $`{\displaystyle \frac{3}{4R^2}}\left(\xi _{lm}^{J+1}+\xi _{lm}^{J1}\right)={\displaystyle \frac{l(l+1)}{2R^2}}\left({\displaystyle \frac{\varphi _{lm}^J(R)}{R}}{\displaystyle \frac{\varphi _{lm}^J}{r}}|_{r=R}\right).`$
This equation in $`\xi _{lm}^{J+1}`$ is solved and, for each pair $`(l,m)`$, we impose for $`\varphi _{lm}^{J+1}`$

$$\left(\frac{}{t}+\frac{}{r}+\frac{1}{r}\right)\varphi _{lm}^{J+1}(r)|_{r=R}=\xi _{lm}^{J+1},$$

which looks like a modification of the condition (4).

### 3.2 Tests on outgoing waves

The Sommerfeld boundary condition (4) is an exact condition, even at finite distance from the source, when only considering monopolar waves. In order to test our implementation of absorbing boundary condition (8), we compared its efficiency in being transparent to waves carrying only monopolar, dipolar and quadrupolar terms, to the efficiency of the Sommerfeld boundary condition for monopolar waves. We started with $`\varphi =0`$ at $`t=0`$ and then solved Eq. (1) with

$`\sigma (t,r,\theta ,\phi )`$ $`=`$ $`S(r,\theta ,\phi )e^{1/t^2}e^{1/(t1)^2}0t1`$ (19)
$`\sigma (t,r,\theta ,\phi )`$ $`=`$ $`0\text{ otherwise},`$
with $`S(r,\theta ,\phi )`$ null for $`r>R`$.

In all cases, we performed a first calculation with a very large grid (considered as infinite, we checked with various values of the radius that the result in the interval $`0rR`$ would be the same), so that in the time interval $`[0,2R+1]`$ the wave would not reach the boundary, on which we imposed an homogeneous boundary condition111results obtained here did not depend on the nature of boundary conditions. This gave us the reference solution crossing the $`r=R`$ sphere without any reflection. We then solved again the same problem, but on a grid of radius $`R`$, imposing Sommerfeld boundary conditions $`B_1\varphi =0`$ (4), or our quadrupolar boundary conditions $`B_3\varphi =0`$ through the system (16)-(17). The $`L_1`$ norm of the relative difference between the functions obtained on the small grid and the reference solution was taken as the error.

#### 3.2.1 $`l2`$ case

First, we took

$$S(r,\theta ,\phi )=\left(e^{r^2}e^{R^2}\right)\left(r^2\mathrm{cos}^2\theta +r\mathrm{sin}\theta \mathrm{cos}\phi \right),$$
(20)
which contains only $`l2`$ modes. Figure 1 shows the relative efficiency of $`B_3\varphi =0`$ (16) condition to $`B_1\varphi =0`$ (4) for all three modes present in the wave generated by (20). For the monopolar ($`l=0`$) mode, the evolution of the error would be the same for both types of boundary conditions, within one percent of difference on the error. As far as the discrepancy for dipolar and quadrupolar modes is concerned, one can see that it drops from $`10^4`$ with Sommerfeld boundary condition, to $`10^{12}`$ with $`B_3\varphi =0`$ (16). This lower level is the same as for the monopolar mode with the Sommerfeld boundary condition. We have checked that all solutions had converged with respect to the number of spectral coefficients and to the time-step. The error level at $`10^{12}`$ is then mainly due to the condition number of the matrix operator we invert (see Sec. 3.1 above). We here conclude that our formulation of $`B_3\varphi =0`$ (16) is as efficient for waves containing only $`l2`$ modes as the Sommerfeld boundary condition (4) for monopolar waves.

#### 3.2.2 Waves containing higher multipoles

The study has been extended to a more general source $`\sigma `$ which contains a priori all multipolar terms:

$$S(r,\theta ,\phi )=\left(e^{r^2}e^{R^2}\right)\left(e^{4(x0.7)^2}+e^{3(x+0.5)^2}\right).$$
(21)
Of course, in numerical implementation, only a finite number of these terms are represented. The geometry of this source can be related to the distribution of mass in the case of a binary system of gravitating bodies, which is one of the main astrophysical sources of gravitational radiation we try to model. Let us make a comparison between the errors obtained, on the one hand with the condition $`B_1\varphi =0`$ (Figure 2), and on the other hand with $`B_3\varphi =0`$ (Figure 3).

As in the case in Figure 1, the error in the monopolar component remains roughly the same, regardless of whether one uses boundary condition (4) or (16). The errors for the dipolar and quadrupolar components also exhibit similar properties: the use of condition (16) causes these errors to be of the same magnitude as the error in the monopolar term. In the case of Figure 3, this level is higher than on Figure 1 because a longer time-step has been used. Finally, we have also plotted the discrepancies between the reference and test solutions for the $`l=3`$ multipole. Following , the boundary condition $`B_3\varphi =0`$ is not exact for this component. Nevertheless, one can see a reduction in the error for this component. This can be understood using the result of  which shows that the condition $`B_3\varphi =0`$ cancels the first 3 terms in the asymptotic development in powers of $`1/r`$ of the solution $`\varphi `$ (9). Then, since a given multipolar term $`l_0`$ is present in terms like $`1/r^n`$ with $`nl_0`$ (see e.g. ), it is clear that the condition $`B_3\varphi =0`$ is supposed to cancel all terms decaying slower than $`1/r^4`$ in the $`l3`$ mode. Thus, the error displayed on Figure 4 is three orders of magnitude lower with the condition $`B_3\varphi =0`$ than with $`B_1\varphi =0`$.

We have checked this point, namely that the maximal error over the time interval would decrease like $`1/R^4`$, where $`R`$ is the distance at which the boundary conditions were imposed. We have also checked that the error decreased both exponentially with the number of coefficients used in $`r,\theta `$ or $`\phi `$, as one would expect for spectral methods, and like $`\delta t^2`$ (second-order time integration scheme). Figure 4 shows the overall error as a function of time for both boundary conditions used. Comparing Figure 4 with figures 2 and 3, one can see that most of the error comes from the $`l=1`$ term when using $`B_1\varphi =0`$ boundary condition, and from the $`l=3`$ term when using $`B_3\varphi =0`$. Finally, the computational cost of this enhanced boundary condition is very low with this new approach. For the tests presented here, the difference in CPU time would be of about 10%. This is linked with the fact that our formulation (16) is a perturbation of the Sommerfeld boundary condition (4), where the quantity $`\xi _{lm}(t)`$ is obtained by simple (ordinary differential equation) integration.

## 4 Conclusions

The purpose of this paper has been to provide a boundary condition that is well-adapted for the simulation of astrophysical sources of gravitational radiation, whose dominant modes are quadrupolar. We took the series of boundary conditions derived by Bayliss and Turkel , truncated at quadrupolar order, and derived a new formulation of that third-order condition in terms of a first-order condition (resembling the classical radiation one), combined with a wave-like equation on the outer boundary of the integration domain. This formulation is simple in the sense that mixed derivatives are (almost) absent.

The numerical implementation using spectral methods and spherical coordinates is straightforward and this formulation of high-order boundary conditions requires only a little more CPU time (less than 10% in our tests) than the simplest first-order condition (4). We have verified that our implementation of this boundary condition had the same efficiency with respect to transparency for dipolar and quadrupolar waves as the Sommerfeld condition (4) for monopolar waves. The precision increases very rapidly (like $`1/R^4`$) as one imposes the boundary condition further from the source of radiation. These two points are of great interest for the simulation of gravitational radiation from isolated astrophysical sources.

As an alternative, one can cite that more accurate results may be obtained using the so-called 2+2 formalism in the wave zone  and matching it to the results in 3+1 formalism222this is the characteristic-Cauchy matching, CCM near the source. Our approach is different, much simpler to implement and should give accurate enough results for the Einstein equations.


format: text

text:
# Model Selection for Support Vector Machine Classification

## 1 Introduction

Support Vector Machines (SVMs) have emerged in recent years as powerful techniques both for regression and classification. One of the central open questions is model selection: how does one tune the parameters of the SVM algorithm to achieve optimal generalization performance? We focus on the case of SVM classification, where these “hyperparameters” include any parameters appearing in the SVM kernel, as well as the penalty parameter $`C`$ for violations of the margin constraint.

Our aim in this paper is twofold. First, we extend our work on probabilistic methods for SVMs to the case of quadratic slack penalties; we also develop a “baseline” algorithm which can be used to find in principle exact maxima of the evidence. Second, we perform numerical experiments on a selection benchmark data sets to compare the model selection criteria derived from the probabilistic view of SVMs with alternatives that directly try to optimize estimates of test error. Our focus in these experiments is less on computational efficiency, but rather on the relative merits of the methods in terms of the resulting generalization performance.

We begin in Sec. 2 with a brief review of SVM classification and of its probabilistic interpretation; the setup will be such that the extension of the probabilistic point of view to the quadratic penalty case requires only small changes compared to linear penalty SVMs. In Sec. 3 we review some criteria for model selection that have been proposed based on approximations to the test error. We also describe previous approximations to the evidence for the linear penalty SVM, and then give an analogue for quadratic penalty SVMs. Exact expressions for gradients of the evidence with respect to the hyperparameters are then derived in terms of averages over the posterior. Sec. 4 has a description of the methods we use in our numerical experiments on model selection, including the Hybrid Monte Carlo algorithm which we use to calculate evidence gradients numerically. The results of our experiments on benchmark data sets are discussed in Sec. 5; we conclude in Sec. 6 with a summary and an outlook towards future work.

## 2 SVM classification

In this section, we give a very brief review of SVM classification; for details the reader is referred to recent textbooks or review articles such as . We also sketch the probabilistic interpretation of SVMs, from which we later obtain Bayesian criteria for SVM model selection.

Suppose we are given a set $`D`$ of $`n`$ training examples $`(x_i,y_i)`$ with binary outputs $`y_i=\pm 1`$ corresponding to the two classes. The basic SVM idea is to map the inputs $`x`$ to vectors $`\mathit{\varphi }(x)`$ in some high-dimensional feature space; ideally, in this feature space, the problem should be linearly separable. Suppose first that this is true. Among all decision hyperplanes $`𝐰\mathit{\varphi }(x)+b=0`$ which separate the training examples (i.e. which obey $`y_i(𝐰\mathit{\varphi }(x_i)+b)>0`$ for all $`x_iX`$, $`X`$ being the set of training inputs), the SVM solution is chosen as the one with the largest margin, i.e. the largest minimal distance from any of the training examples. Equivalently, one specifies the margin to be equal to 1 and minimizes the squared length of the weight vector $`𝐰^2`$ , subject to the constraint that $`y_i(𝐰\mathit{\varphi }(x_i)+b)1`$ for all $`i`$. The quantities $`y_i(𝐰\mathit{\varphi }(x_i)+b)`$ are again called margins, although for an unnormalized weight vector they no longer represent geometrical distances . This leads to the following optimization problem: Find a weight vector $`𝐰`$ and an offset $`b`$ such that $`\frac{1}{2}𝐰^2`$ is minimized, subject to the constraint that $`y_i(𝐰\mathit{\varphi }(x_i)+b)1`$ for all training examples.

If the problem is not linearly separable, or if one wants to avoid fitting noise in the training data, ‘slack variables’ $`\xi _i0`$ are introduced which measure how much the margin constraints are violated; one thus writes $`y_i(𝐰\mathit{\varphi }(x_i)+b)1\xi _i`$. To control the amount of slack allowed, a penalty term $`(C/p)_i\xi _i^p`$ is then added to the objective function $`\frac{1}{2}𝐰^2`$, with a penalty coefficient $`C`$. Common values for the exponent parameter are $`p=1`$ and $`p=2`$, giving linear and quadratic slack penalties respectively. Training examples with $`y_i(𝐰\mathit{\varphi }(x_i)+b)1`$ (and hence $`\xi _i=0`$) incur no penalty; the others contribute $`(C/p)[1y_i(𝐰\mathit{\varphi }(x_i)+b)]^p`$ each. This gives the SVM optimization problem: Find $`𝐰`$ and $`b`$ to minimize

$$\frac{1}{2}𝐰^2+C\underset{i}{}l_p(y_i[𝐰\mathit{\varphi }(x_i)+b])$$
(1)
where $`l_p(z)`$ is the loss function

$$l_p(z)=\frac{1}{p}(1z)^pH(1z)$$
(2)
The Heaviside step function $`H(1z)`$ (defined as $`H(a)=1`$ for $`a0`$ and $`H(a)=0`$ otherwise) ensures that this is zero for $`z>1`$. For $`p=1`$, $`l_p(z)`$ is called (shifted) hinge loss or soft margin loss.

In the following we modify the basic SVM problem by adding the quadratic term $`\frac{1}{2}b^2/B^2`$ to (1), thus introducing a penalty for large offsets $`b`$. A discussion of why this is reasonable, certainly within a probabilistic view, can be found in ; at any rate the standard formulation can always be retrieved by making the constant $`B`$ large. We can now define an augmented weight vector $`\stackrel{~}{𝐰}=(b/B,𝐰)`$ and augmented feature space vectors $`\stackrel{\mathbf{~}}{\mathit{\varphi }}(x)=(B,\mathit{\varphi }(x))`$ so that the modified SVM problem is to find a $`\stackrel{~}{𝐰}`$ which minimizes

$$\frac{1}{2}\stackrel{~}{𝐰}^2+C\underset{i}{}l_p(y_i\stackrel{~}{𝐰}\stackrel{\mathbf{~}}{\mathit{\varphi }}(x_i))$$
(3)
This statement of the problem is useful for the probabilistic interpretation of SVM classification, of which more shortly. For a practical solution, one uses Lagrange multipliers $`\alpha _i`$ conjugate to the constraints $`y_i\stackrel{~}{𝐰}\stackrel{\mathbf{~}}{\mathit{\varphi }}(x_i)1\xi _i`$ and finds in the standard way (see e.g. ) that the optimal (augmented) weight vector is $`\stackrel{~}{𝐰}^{}=_iy_i\alpha _i\mathit{\varphi }(x_i)`$. For the linear penalty case $`p=1`$, the $`\alpha _i`$ are found from

$$\underset{0\alpha _iC}{\mathrm{max}}\left(\underset{i}{}\alpha _i\frac{1}{2}\underset{i,j}{}\alpha _i\alpha _jy_iy_jK_{ij}\right)$$
(4)
Here $`K_{ij}=K(x_i,x_j)`$ are the elements of the Gram matrix $`𝐊`$, obtained by evaluating the kernel $`K(x,x^{})=\stackrel{\mathbf{~}}{\mathit{\varphi }}(x)\stackrel{\mathbf{~}}{\mathit{\varphi }}(x^{})=\mathit{\varphi }(x)\mathit{\varphi }(x^{})+B^2`$ for all pairs of training inputs. The corresponding optimal latent or discrimination function is $`\theta ^{}(x)=\stackrel{~}{𝐰}^{}\stackrel{\mathbf{~}}{\mathit{\varphi }}(x)=_iy_i\alpha _iK(x,x_i)`$. Only the $`x_i`$ with $`\alpha _i>0`$ contribute to this sum; these are called support vectors (SVs). SVs fall into two groups: If $`\alpha _i<C`$, one has $`y_i\theta _i^{}y_i\theta ^{}(x_i)=1`$; we will call these the “marginal SVs” because their margins are exactly at the allowed limit where no slack penalty is yet incurred. For $`\alpha _i=C`$, on the other hand, $`y_i\theta _i^{}1`$, and these “hard SVs” are the points at which the slack penalty is active. Non-SVs have large margins, $`y_i\theta _i^{}1`$.

For the quadratic penalty case $`p=2`$, the $`\alpha _i`$ are obtained as the solution of (see e.g. )

$$\underset{0\alpha _i}{\mathrm{max}}\left(\underset{i}{}\alpha _i\frac{1}{2}\underset{i,j}{}\alpha _i\alpha _jy_iy_jK_{ij}^C\right)$$
(5)
where $`K_{ij}^C=K_{ij}+C^1\delta _{ij}`$. Apart from the replacement of $`K`$ by $`K^C`$, this maximization problem is the same as (4) for the linear penalty case in the limit $`C\mathrm{}`$ where no violations of the margin constraints are allowed. There is now only one kind of SV, identified by $`\alpha _i>0`$. It follows by differentiating (5) that for a SV one has $`y_i_j\alpha _jy_jK_{ij}^C=1`$. The margin for a SV is thus $`y_i\theta _i^{}=y_i_j\alpha _jy_jK_{ij}=1\alpha _i/C`$, so that all SVs incur a nonzero slack penalty. Non-SVs again have $`y_i\theta _i^{}1`$.

We now turn to the probabilistic interpretation of SVM classification (see Refs.  and the works quoted below). The aim of such an interpretation is to allow the application of Bayesian methods to SVMs, without modifying the basic SVM algorithm which already has a large user community. (An alternative philosophy would be to consider similar inference algorithms which share some of the benefits of SVMs but are constructed directly from probabilistic models; Tipping’s Relevance Vector Machine  is a successful example of this.) One regards (3) as defining a negative log-posterior probability for the parameters $`\stackrel{~}{𝐰}`$ of the SVM, given a training set $`D`$. The conventional SVM classifier is then interpreted as the maximum a posteriori (MAP) solution of the corresponding probabilistic inference problem. The first term in (3) gives the prior $`Q(\stackrel{~}{𝐰})\mathrm{exp}(\frac{1}{2}\stackrel{~}{𝐰}^2)`$. This is a Gaussian prior on $`\stackrel{~}{𝐰}`$; the components of $`\stackrel{~}{𝐰}`$ are uncorrelated with each other and have unit variance. Because only the ‘latent function’ values $`\theta (x)=\stackrel{~}{𝐰}\stackrel{\mathbf{~}}{\mathit{\varphi }}(x)`$—rather than $`\stackrel{~}{𝐰}`$ itself—appear in the second, data dependent term of (3), it makes sense to express the prior directly as a distribution over these. The $`\theta (x)`$ have a joint Gaussian distribution because the components of $`\stackrel{~}{𝐰}`$ do, with covariances given by

$$\theta (x)\theta (x^{})=(\stackrel{\mathbf{~}}{\mathit{\varphi }}(x)\stackrel{~}{𝐰})(\stackrel{~}{𝐰}\stackrel{\mathbf{~}}{\mathit{\varphi }}(x^{}))=K(x,x^{})$$

The SVM prior is therefore simply a Gaussian process (GP) over the functions $`\theta `$, with zero mean and with the kernel $`K(x,x^{})`$ as covariance function. This link between SVMs and GPs has been pointed out by a number of authors, e.g. . It can be understood from the common link to reproducing kernel Hilbert spaces , and can be extended from SVMs to more general kernel methods . For connections to regularization operators see also . A nice introduction to inference with Gaussian processes can be found in Ref. .

The second term in (3) similarly becomes a (negative) log-likelihood if we define the probability of obtaining output $`y`$ for a given $`x`$ (and $`\theta `$) as

$$Q(y=\pm 1|x,\theta )=\kappa (C)\mathrm{exp}[Cl_p(y\theta (x))]$$
(6)
The constant factor $`\kappa (C)`$ is determined from $`\kappa ^1(C)=\mathrm{max}_z\{e^{Cl_p(z)}+e^{Cl_p(z)}\}`$ to ensure that $`_{y=\pm 1}Q(y|x,\theta )1`$. In the linear penalty case this gives $`\kappa (C)=1/[1+\mathrm{exp}(2C)]`$; for the quadratic penalty SVM, the maximum in the definition of $`\kappa ^1(C)`$ is assumed at a value of $`\theta `$ obeying $`\theta =\mathrm{tanh}(C\theta )`$ and so $`\kappa (C)`$ can easily be found numerically. The likelihood for the complete data set (more precisely, for the training outputs $`Y=(y_1\mathrm{}y_n)`$ given the training inputs $`X`$) is then

$$Q(Y|X,\theta )=\underset{i}{}Q(y_i|x_i,\theta )$$

With these definitions eq. (3) is, up to unimportant constants, equal to the log-posterior111In (7) the unrestricted sum over $`x`$ runs over all possible inputs, and $`K^1(x,x^{})`$ are the elements of the inverse of $`K(x,x^{})`$, viewed as a matrix. We assume here that the input domain is discrete. This avoids mathematical subtleties with the definition of determinants and inverses of operators (rather than matrices), while maintaining a scenario that is sufficiently general for all practical purposes.

$$\mathrm{ln}Q(\theta |X,Y)=\frac{1}{2}\underset{x,x^{}}{}\theta (x)K^1(x,x^{})\theta (x^{})C\underset{i}{}l_p(y_i\theta (x_i))+\text{const}$$
(7)
By construction, the maximum of $`\theta ^{}(x)`$ gives the conventional SVM classifier, and this is easily verified explicitly .

The probabilistic model defined above is not normalized, since $`_{y=\pm 1}Q(y|x,\theta )<1`$ for generic values of $`\theta (x)`$. The implications of this have been previously discussed in detail . The normalization of the model is important for the theoretical justification of tuning hyperparameters via maximization of the data likelihood or “evidence”. Nevertheless, experiments in  showed that promising results for hyperparameter optimization could be obtained also with the unnormalized version of the model. We will therefore, in common with other work on probabilistic interpretations of SVMs , disregard the normalization issue from now on. We will also focus on SVM classifiers constructed from radial basis function (RBF) kernels

$$K(x,x^{})=k_0\mathrm{exp}\left[\underset{a}{}\frac{(x^ax^a)^2}{2l_a^2}\right]+k_{\mathrm{off}}$$
(8)
where the $`x^a`$ are the different input components, $`k_0`$ is the kernel amplitude and $`k_{\mathrm{off}}`$ the kernel offset; $`k_{\mathrm{off}}`$ corresponds to the term $`B^2`$ discussed above that arises by incorporating the offset $`b`$ into the kernel. Each input dimension has associated with it a length scale $`l_a`$. Since in the probabilistic interpretation $`K(x,x^{})`$ is the prior covariance function of the latent function $`\theta (x)`$, each $`l_a`$ determines the distance in the $`x^a`$-direction over which $`\theta (x)`$ is approximately constant; large $`l_a`$ correspond to an input component of little relevance (see e.g. ).

## 3 Model selection criteria

### 3.1 Error bounds and approximations

Model selection aims to tune the hyperparameters of SVM classification (the penalty parameter $`C`$ and any kernel parameters) in order to achieve the lowest test error $`ϵ`$, i.e. the lowest probability of misclassification of unseen test examples. The test error is not observable directly, and so one is lead to use bounds or approximations as model selection criteria. The simplest such bounds  which have been applied as model selection criteria  are expressed in terms of the quantity

$$\frac{R^2}{n}\underset{i}{}\alpha _i$$
(9)
Here $`R`$ is the radius of the smallest ball in feature space containing all training examples, while $`_i\alpha _i`$ can be shown to equal the inverse square of the distance between the separating hyperplane and the closest training points. For RBF kernels, $`R`$ is bounded by a constant since every input point has the same squared distance $`\stackrel{\mathbf{~}}{\mathit{\varphi }}(x)\stackrel{\mathbf{~}}{\mathit{\varphi }}(x)=K(x,x)`$ from the origin.

More recent work has shown that better bounds and approximations can be obtained for the leave-out-out error $`ϵ_{\mathrm{loo}}`$. If $`\theta ^i(x)`$ is the latent function obtained by training the SVM classifier on the data set with example $`(x_i,y_i)`$ left out, then $`ϵ_{\mathrm{loo}}`$ is the probability of misclassification of the left-out example if this procedure is applied to each data point in turn,

$$ϵ_{\mathrm{loo}}=\frac{1}{n}\underset{i}{}H(y_i\theta _i^i)$$
(10)
where we have abbreviated $`\theta _i^i\theta ^i(x_i)`$. Averaged over data sets this is an unbiased estimate of the average test error that is obtained from training sets of $`n1`$ examples. This says nothing about the variance of this estimate; nevertheless, one may hope that $`ϵ_{\mathrm{loo}}`$ is a reasonable proxy for the test error that one wishes to optimize. (This is in contrast to the training error, i.e. the fraction of all $`n`$ training examples misclassified when training on the complete data set, which is general a strongly biased estimate of test error.) For large data sets, $`ϵ_{\mathrm{loo}}`$ is time-consuming to compute and one is driven to look for cheaper bounds or approximations. Since removing non-SVs from the data set does not change the SVM classifier, a trivial bound on $`ϵ_{\mathrm{loo}}`$ is the sum of the training error and the fraction of support vectors, both obtained when training on all $`n`$ examples. To get better bounds, one writes

$$ϵ_{\mathrm{loo}}=\frac{1}{n}\underset{i}{}H\left(y_i[\theta _i^{}\theta _i^i]y_i\theta _i^{}\right)$$

which shows that an upper bound on $`y_i[\theta _i^{}\theta _i^i]`$ will give an upper bound on $`ϵ_{\mathrm{loo}}`$. Jaakkola and Haussler proved a bound of this form, $`y_i[\theta _i^{}\theta _i^i]\alpha _iK_{ii}`$; as before, the $`\alpha _i`$ are those obtained from training on the full data set. More sophisticated bounds were given by Chappelle and Vapnik  in terms of what they called the “span”. We focus on the case of quadratic penalty SVMs, where the span estimates are simplest to state. In the simplified version of Ref. , and adapting to our formulation with the offset $`b`$ incorporated into the kernel, the span $`S_i`$ for a support vector can be defined as

$$S_i^2=\underset{\lambda }{\mathrm{min}}\underset{j,k}{}\lambda _j\lambda _kK_{jk}^C$$
(11)
where the minimum is over all $`\lambda =(\lambda _1\mathrm{}\lambda _n)`$ with $`\lambda _i=1`$, and $`\lambda _j=0`$ whenever $`\alpha _j=0`$. With this definition, one can calculate $`y_i[\theta _i^{}\theta _i^i]`$ exactly under the assumption that dropping the point $`x_i`$ from the training set leaves the “SV set” unchanged, in the sense that no new SVs arise in the new classifier and that all old SVs remain. One thus finds $`y_i[\theta _i^{}\theta _i^i]=\alpha _i(S_i^21/C)`$. $`S_i^2`$ can also be worked out explicitly as $`S_i^2=1/[(𝐊_{\mathrm{SV}}+𝐈/C)^1]_{ii}`$, where $`𝐊_{\mathrm{SV}}`$ is the Gram matrix $`𝐊`$ restricted to the SVs and $`𝐈`$ is the unit matrix. (The same result was obtained by Opper and Winther  using a slightly different approach.) Using finally that $`y_i\theta _i^{}=1\alpha _i/C`$ for quadratic penalty SVMs, one thus has

$$ϵ_{\mathrm{loo}}ϵ_{\mathrm{span}}=\frac{1}{n}\underset{i}{}H\left(\alpha _iS_i^21\right),S_i^2=1/[(𝐊_{\mathrm{SV}}+𝐈/C)^1]_{ii}$$
(12)
This is only an approximation because the assumption of an unchanged SV set will not hold for every SV removed from the training set.

The span estimate (12) of leave-one-out error has the undesirable property of being discontinuous as hyperparameters are varied, making numerical optimization difficult. The discontinuity arises from the discontinuity in the Heaviside step function $`H`$, and from the fact that the size of the matrix $`𝐊_{\mathrm{SV}}`$ changes as training examples enter or leave the set of SVs. To get around this , one can approximate $`H(z)`$ by a sigmoidal function $`1/[1+\mathrm{exp}(c_1x+c_2)]`$ and smooth the span by adding a penalty that forces any nonzero $`\lambda _j`$ to go to zero when $`\alpha _j0`$. This gives the modified span definition

$$S_i^2=\underset{\lambda }{\mathrm{min}}\underset{j,k}{}\lambda _j\lambda _kK_{jk}^C+\eta \underset{ji}{}\frac{\lambda _j^2}{\alpha _j}$$

with the minimum taken over the same $`\lambda `$ as in (11). Explicitly, one finds

$$S_i^2=\frac{1}{[(𝐊_{\mathrm{SV}}+𝐈/C+\eta 𝐀_{\mathrm{SV}}^1)^1]_{ii}}\frac{\eta }{\alpha _i}$$

where $`𝐀_{\mathrm{SV}}`$ is the diagonal matrix containing the nonzero $`\alpha _i`$. This is easily seen to be continuous even when the set of SVs changes as hyperparameters are varied. For $`\eta 0`$ one recovers the original span definition (11); for $`\eta \mathrm{}`$, on the other hand, $`S_i^2K_{ii}^C=K_{ii}+1/C`$ and one recovers the Jaakkola and Haussler bound. Overall, the smoothed span estimate for $`ϵ_{\mathrm{loo}}`$ contains three smoothing parameters $`c_1`$, $`c_2`$ and $`\eta `$.

For linear penalty SVMs, Wahba  considered a modified version of $`ϵ_{\mathrm{loo}}`$, obtained by replacing the Heaviside step function $`H(z)`$ in (10) by the hinge loss $`l_1(z)=(1z)H(1z)`$; since $`l_1(z)H(z)`$, this actually gives an upper bound on $`ϵ_{\mathrm{loo}}`$. Wahba’s GACV (generalized approximate cross-validation) estimate for this modified $`ϵ_{\mathrm{loo}}`$ is

$$ϵ_{\mathrm{gacv}}=\frac{1}{n}\underset{i}{}\left[l_1(y_i\theta _i^{})+\alpha _iK_{ii}f(y_i\theta _i^{})\right]$$
(13)
where

$$f(z)=\{\begin{array}{cc}2\hfill & x<1\hfill \\ 1\hfill & 1x1\hfill \\ 0\hfill & x>1\hfill \end{array}$$

The first term in the sum in (13) would just give the naive estimate of the (modified) $`ϵ_{\mathrm{loo}}`$ from the performance on the training set; the second term effectively corrects for the bias in this estimate. Because of the nature of the function $`f`$, $`ϵ_{\mathrm{gacv}}`$ can exhibit discontinuities as hyperparameters are varied and this has to be taken into account when minimizing it numerically.

### 3.2 Approximations to the evidence

From a probabilistic point of view, it is natural to tune hyperparameters to maximize the likelihood of the data $`Q(Y|X)`$. By definition, $`Q(Y|X)=𝑑\theta Q(Y|X,\theta )Q(\theta )`$ where the integration is over the values $`\theta (x)`$ of the latent function $`\theta `$ at all different input points $`x`$. The likelihood $`Q(Y|X,\theta )`$ only depends on the values $`\theta _i\theta (x_i)`$ of $`\theta `$ at the training inputs; all other $`\theta (x)`$ can be integrated out trivially, so that

$$Q(Y|X)=𝑑𝜽Q(Y|X,𝜽)Q(𝜽)$$

where now the integral is over the $`n`$-dimensional vector $`𝜽=(\theta _1\mathrm{}\theta _n)`$. Because $`Q(\theta )`$ is a zero mean Gaussian process, the marginal $`Q(𝜽)`$ is a zero mean Gaussian distribution with covariance matrix $`𝐊`$. The evidence is therefore

$$Q(Y|X)=|2\pi 𝐊|^{1/2}\kappa ^n(C)𝑑𝜽\mathrm{exp}\left[\frac{1}{2}𝜽^\mathrm{T}𝐊^1𝜽\underset{i}{}Cl_p(y_i\theta _i)\right]$$
(14)
This n-dimensional integral is in general impossible to carry out exactly. But it can be approximated by expanding the exponent around its maximum $`𝜽^{}`$, the solution of the SVM optimization problem. For the linear penalty case, this requires some care: because of the kink in the hinge loss $`l_1(z)`$, the $`\theta _i`$ corresponding to marginal SVs have to be treated separately. The result for the normalized log-evidence, suitably smoothed to avoid spurious singularities, is

$$E(Y|X)\frac{1}{n}\mathrm{ln}Q(Y|X)=\frac{1}{n}\mathrm{ln}Q^{}(Y|X)\frac{1}{2n}\mathrm{ln}det(𝐈+𝐋_m𝐊_m)$$
(15)
where $`𝐊_m`$ is the sub-matrix of the Gram matrix corresponding to the marginal SVs, $`𝐋_m`$ is the diagonal matrix with entries $`2\pi [\alpha _i(C\alpha _i)/C]^2`$ and

$$\frac{1}{n}\mathrm{ln}Q^{}(Y|X)=\frac{1}{2n}\underset{i}{}\alpha _iy_i\theta _i^{}\frac{C}{n}\underset{i}{}l_p(y_i\theta _i^{})+\mathrm{ln}\kappa (C)$$
(16)
The approximation (15) is computationally efficient because it only involves calculating the determinant of a single matrix of size equal to the number of marginal SVs. A related approximation was proposed by Kwok . He suggested to smooth the kink in the hinge loss by using a sigmoidal approximation for the Heaviside function, giving $`l_1(z)s(z)=(1z)/[1+\mathrm{exp}[c(1z)]`$ with a smoothing parameter $`c`$. This has the disadvantage that the SVM solution $`\theta ^{}`$ is no longer a maximum of the smoothed posterior. The analogous result to (15) also involves, instead of $`𝐊_m`$ and $`𝐋_m`$, the whole Gram matrix and a diagonal matrix with entries $`Cs^{\prime \prime }(y_i\theta _i^{})`$, respectively. One thus needs either to evaluate a large determinant of size $`n`$, or—somewhat arbitrarily—to truncate small values of $`s^{\prime \prime }(y_i\theta _i^{})`$ to zero to reduce the size of the problem. We therefore do not consider this approach further.

An approximation similar to (15) can easily be derived for the quadratic penalty case. The loss function $`l_2(z)`$ now has a continuous first derivative and so all $`\theta _i`$ can be treated on the same footing. The derivation of the Laplace approximation is thus standard, and involves the Hessian of the log-likelihood at the maximum $`𝜽^{}`$; the resulting approximation to the (normalized log-) evidence is

$$E(Y|X)=\frac{1}{n}\mathrm{ln}Q^{}(Y|X)\frac{1}{2n}\mathrm{ln}det(𝐈+𝐌_{\mathrm{SV}}𝐊_{\mathrm{SV}})$$
(17)
where $`𝐌_{\mathrm{SV}}`$ is a diagonal matrix containing the second derivatives $`Cl_2^{\prime \prime }(y_i\theta _i^{})`$ of the loss function evaluated for all the SVs. The calculation of $`E(Y|X)`$ according to (17) requires only the determinant of a matrix whose size is the number of SVs, again expected to be manageably small. The matrix $`𝐌_{\mathrm{SV}}`$ as defined above is just a multiple of the unit matrix, since $`l_2^{\prime \prime }(z)=H(1z)`$ and $`z=y_i\theta _i^{}<1`$ for SVs. However, the step discontinuity in $`l_2^{\prime \prime }(z)`$ at $`z=1`$ has the undesirable consequence that the Laplace approximation to the evidence will jump discontinuously when one or several of the $`\alpha _i`$ reach zero as hyperparameters are varied. We therefore smooth the result by replacing $`l_2^{\prime \prime }(z)=H(1z)`$ in the definition of $`𝐌_{\mathrm{SV}}`$ by the approximation $`l_2^{\prime \prime }(z)\mathrm{exp}[a/(1z)]`$ for $`z<1`$ (and $`0`$ for $`z1`$). This is smooth at $`z=1`$ and also has continuous derivatives of all orders at this point. The value of $`a`$ determines the range of values of $`z=y_i\theta _i`$ around 1 for which the smoothing is significant, with $`a0`$ recovering the Heaviside step function.

### 3.3 Evidence gradients

Beyond the relatively simple approximations to the evidence derived above, it is difficult to obtain accurate numerical estimates of the evidence. This is a well-known general problem: while averages over probability distributions are straightforward to obtain, normalization constants for such distributions – such as the evidence, which is the normalization factor for the posterior – require much greater numerical effort (see e.g. ). To avoid this problem, one can estimate the gradients of the evidence with respect to the hyperparameters and use these in a gradient ascent algorithm, without ever calculating the value of the evidence itself. As we show in this section, these gradients can be expressed as averages over the posterior distribution, which one can then estimate by sampling as explained in Sec. 4.2.

Starting from eq. (14) we can find the derivative of the normalized log-evidence $`E(Y|X)=n^1\mathrm{ln}Q(Y|X)`$ w.r.t. the penalty (or noise) parameter C:

$`{\displaystyle \frac{}{C}}E(Y|X)`$ $`=`$ $`{\displaystyle \frac{\mathrm{ln}\kappa (C)}{C}}{\displaystyle \frac{𝑑𝜽Q(𝜽)_il_p(y_i\theta _i)\mathrm{exp}\left[C_il_p(y_i\theta _i)\right]}{𝑑𝜽Q(𝜽)\mathrm{exp}\left[C_il_p(y_i\theta _i)\right]}}`$ (18)
$`=`$ $`{\displaystyle \frac{\mathrm{ln}\kappa (C)}{C}}{\displaystyle \frac{1}{n}}{\displaystyle \underset{i}{}}l_p(y_i\theta _i)`$
where the average is, as expected, over the posterior $`Q(𝜽|D)Q(Y|X,𝜽)Q(𝜽)`$. Similarly, the derivative of the log-evidence w.r.t. any parameter $`\lambda `$ appearing in the kernel is

$`{\displaystyle \frac{}{\lambda }}E(Y|X)`$ $`=`$ $`{\displaystyle \frac{1}{2n}}{\displaystyle \frac{}{\lambda }}\mathrm{ln}|2\pi 𝐊|{\displaystyle \frac{𝑑𝜽\frac{1}{2}𝜽^\mathrm{T}\frac{}{\lambda }𝐊^1𝜽\mathrm{exp}\left[\frac{1}{2}𝜽^\mathrm{T}𝐊^1𝜽_iCl_p(y_i\theta _i)\right]}{𝑑𝜽\mathrm{exp}\left[\frac{1}{2}𝜽^\mathrm{T}𝐊^1𝜽_iCl_p(y_i\theta _i)\right]}}`$ (19)
$`=`$ $`{\displaystyle \frac{1}{2n}}\mathrm{tr}\left({\displaystyle \frac{𝐊}{\lambda }}𝐊^1\right)+{\displaystyle \frac{1}{2n}}𝜽^\mathrm{T}𝐊^1{\displaystyle \frac{𝐊}{\lambda }}𝐊^1𝜽`$
$`=`$ $`{\displaystyle \frac{1}{2n}}\mathrm{tr}\left({\displaystyle \frac{𝐊}{\lambda }}𝐊^1𝐈𝜽𝜽^\mathrm{T}𝐊^1\right)`$
Numerical evaluation of this expression as it stands would be unwise, since the difference $`𝐈𝜽𝜽^\mathrm{T}𝐊^1`$ can be much smaller than the two contributions individually; in fact, for $`n=0`$ we know that it is exactly zero. It is better to rewrite (19), using the fact that the elements of the matrix $`𝐈𝜽𝜽^\mathrm{T}𝐊^1`$ can be obtained as

$$\delta _{ij}\theta _i(𝐊^1𝜽)_j=\mathrm{exp}\left[\frac{1}{2}𝜽^\mathrm{T}𝐊^1𝜽\right]\frac{}{\theta _j}\theta _i\mathrm{exp}\left[\frac{1}{2}𝜽^\mathrm{T}𝐊^1𝜽\right]$$

The posterior average can thus be worked out using integration by parts, giving

$$\delta _{ij}\theta _i(𝐊^1𝜽)_j=Cl_p^{}(y_j\theta _j)y_j\theta _i$$

If we define the matrix $`𝐘`$ as the diagonal matrix with entries $`y_i`$ so that $`(𝐘𝜽)_i=y_i\theta _i`$, and denote by $`l_p^{}(𝐘𝜽)`$ the vector with entries $`l_p^{}(y_i\theta _i)`$, then this can be written in the compact form

$$𝐈𝜽𝜽^\mathrm{T}𝐊^1=C𝜽[l_p^{}(𝐘𝜽)]^\mathrm{T}𝐘)$$

Combining this with equation (19), one has finally

$$\frac{}{\lambda }E(Y|X)=\frac{C}{2n}[l_p^{}(𝐘𝜽)]^\mathrm{T}𝐘\frac{𝐊}{\lambda }𝐊^1𝜽$$
(20)
This expression appears to require the inverse $`𝐊^1`$ of the Gram matrix, which for large $`n`$ would be computationally expensive to evaluate; however, as described in Sec. 4.2 the sampling from the posterior using Hybrid Monte Carlo can be arranged so that samples of both $`𝜽`$ and $`𝐊^1𝜽`$ are obtained without requiring explicit matrix inversions.

## 4 Numerical methods

In Sec. 5 below, we report results from SVM model selection experiments using the criteria described above. Specifically, for linear penalty SVMs we compare maximizing the Laplace approximation to the evidence $`E(Y|X)`$, eq. (15), with minimizing Wahba’s $`ϵ_{\mathrm{gacv}}`$, eq. (13); for quadratic penalty SVMs we again use the relevant approximation to the evidence, eq. (17), and contrast with minimization of the span estimate $`ϵ_{\mathrm{span}}`$ of the leave-one-out error, eq. (12). These four model selection criteria are “simple” in the sense that they can be evaluated explicitly at moderate computational cost. In order to be able to compare the different criteria directly, and because one of them ($`ϵ_{\mathrm{gacv}}`$) has possible discontinuities as a function of the hyperparameters, we use a simple greedy random walk algorithm for optimization that is described in Sec. 4.1. For the other three criteria, more efficient gradient-based optimization algorithms can be designed  but since our focus here is not on computational efficiency we do not consider these.

For linear penalty SVMs we also studied evidence optimization using numerical estimates of the evidence gradients (18,20). The Monte Carlo method used to perform the necessary averages over the posterior is outlined in Sec. 4.2, while Sec. 4.3 describes the details of the gradient ascent algorithm. Note that our use of evidence gradients provides a baseline for model selection methods based on approximations to the evidence since it locates, up to small statistical errors from the Monte Carlo sampling of posterior averages, a local maximum of the exact evidence.

In all experiments using approximations to the test error ($`ϵ_{\mathrm{span}}`$ and $`ϵ_{\mathrm{gacv}}`$) as model selection criteria, the hyperparameters being optimized were the parameters of the RBF kernel (8), i.e. the amplitudes $`k_0`$, $`k_{\mathrm{off}}`$ and the logarithms of the length scales $`l_a`$ (one per input dimension). The penalty parameter $`C`$ can be fixed to e.g. $`C=1`$ since these criteria depend only on the properties of the SVM solution (i.e. the maximum of the posterior, rather than the whole posterior distribution). This SVM solution only depends on the product $`CK(x,x^{})`$ rather than $`C`$ and the kernel individually, as one easily sees from (4,5) or equivalently from (7). In contrast, the evidence (14) takes into account both the position of the posterior maximum and the shape of the posterior distribution around this maximum; the latter does depend on $`C`$. We therefore include $`C`$ as a hyperparameter to be optimized in evidence maximization. The SVM predictor of the final selected model will of course again be dependent only on the product $`CK(x,x^{})`$; but the value of $`C`$ itself would be important e.g. for the determination of predictive class probabilities. This issue, which we do not pursue here, is discussed in detail in Ref. .

### 4.1 Optimization of “simple” model selection criteria

We optimized the simple model selection criteria using a simple greedy random walk, or “zero temperature Monte Carlo” search. This is a simple adaptation of the common Metropolis Algorithm (see e.g. ) used to sample from a probability distribution; in the zero temperature limit the algorithm reduces to repeatedly adding a small step (which we take to be Gaussian) to each parameter, recalculating the quantity being optimized, and moving to the new point if and only if the new point yields a better value (higher for the evidence, and lower for error estimates). The randomness in the algorithm may appear disadvantageous in terms of computational efficiency; but for our purposes, it is actually helpful since it allows us to assess whether the model selection criteria in question have a number of local optima or a single (global) optimum. It also made further randomization over the initial hyperparameter values unnecessary, and so the experiments with the simple model selection criteria were all started with a fixed set of initial values for the SVM hyperparameters. A few preliminary trials were used to choose initial values with an appropriate order of magnitude, and all results reported were initialized with the hyperparameters $`C=1`$, $`l_a=1`$ for all length scales, $`k_0=1`$ and $`k_{\mathrm{off}}=0.1`$.

The span error estimate and the Laplace evidence for quadratic loss each have additional smoothing parameters that had to be selected ($`c_1`$, $`c_2`$ and $`\eta `$ in the span estimate, $`a`$ for the Laplace evidence; see Secs. 3.1 and 3.2 respectively). Appropriate values for these parameters were found by a simple (log) line search in the parameter values. These tests were not done extensively, but the results for SVM model selection did not seem to depend strongly on the values of these parameters as long as they were of a reasonable order of magnitude. For all tests presented here the values used were $`c_1=5`$, $`c_2=0`$, $`\eta =1`$ for the span estimate, and $`a=0.1`$ for the Laplace evidence with quadratic loss.

In the greedy random walk algorithm, the step size used for each hyperparameter was adapted separately by measuring the acceptance rate for proposed changes in the parameter and scaling the step size up or down to keep the acceptance rate close to 50%. Thus a decreasing step size can be taken as one measure of how well the process has converged to an optimum. The search is terminated when either the step size has become very small, or the change to the criterion being optimized becomes very small. It was also found during experimentation that a useful addition to the basic algorithm was to enforce minimum and maximum values of the hyperparameters. Without such bounds the algorithm would occasionally get “stuck” in a plateau region of the model selection criterion where one or more hyperparameters were either very large or very small. Note that for the kernel hyperparameters steps in the random walk were taken in the natural logarithm of the hyperparameter values, as these scale parameters were expected to show a significant range of variation. Steps for $`C`$ were taken in a linear scale, reflecting the smaller range of variation.

### 4.2 Estimating evidence gradients

We used Hybrid Monte Carlo (HMC, see e.g. ) to estimate the posterior averages required in the expressions (18) and (20) for the exact evidence gradients. The HMC algorithm is a standard technique from statistical physics that works by simulating a stochastic dynamics with a Hamiltonian “energy” defined by the target distribution plus a “momentum”, or kinetic energy term. Denoting the momentum variables $`𝐩`$, the Hamiltonian we choose for our case is

$$(𝜽,𝐩)=\frac{1}{2}𝐩^T\mathrm{𝐊𝐩}+\frac{1}{2}𝜽^\mathrm{T}𝐊^1𝜽+V(𝜽),V(𝜽)=C\underset{i}{}l_p(y_i\theta _i)$$
(21)
and the corresponding “Boltzmann” distribution $`P(𝜽,𝐩)\mathrm{exp}[(𝜽,𝐩)]\mathrm{exp}(\frac{1}{2}𝐩^\mathrm{T}\mathrm{𝐊𝐩})Q(𝜽|D)`$ factorizes over $`𝜽`$ and $`𝐩`$, so that samples from $`Q(𝜽|D)`$ can be obtained by sampling from $`P(𝜽,𝐩)`$ and discarding the momenta $`𝐩`$. The $`𝐩`$ are nevertheless important for the algorithm, since they help to ensure a representative sampling of the posterior. An update step in the HMC algorithm consists of two parts. First, one updates a randomly chosen momentum variable $`p_i`$ by Gibbs sampling according to the Gaussian distribution $`\mathrm{exp}(\frac{1}{2}𝐩^\mathrm{T}\mathrm{𝐊𝐩})`$; this will in general change the value of the Hamiltonian. Second, one changes both $`𝜽`$ and $`𝐩`$ by moving along a Hamiltonian trajectory for some specified “time” $`\tau `$; the trajectory is determined by solving an appropriately discretized version of the differential equations

$`{\displaystyle \frac{d\theta _i}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{}{p_i}}=(\mathrm{𝐊𝐩})_i`$ (22)
$`{\displaystyle \frac{dp_i}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{}{\theta _i}}=(𝐊^1𝜽)_i{\displaystyle \frac{V(𝜽)}{\theta _i}}`$ (23)
For an exact solution of these equations, $``$ would remain constant; due to the discretization, small changes in $``$ are possible and one accepts the update of $`𝜽`$ and $`𝐩`$ from the beginning to the end of the trajectory with the usual Metropolis acceptance rule. Iterating these steps the algorithm will, after some initial equilibration period, produce samples from $`P(𝜽,𝐩)`$.

The occurrence of $`𝐊^1`$ in (23) is inconvenient. We circumvent this by introducing $`\stackrel{\mathbf{~}}{𝜽}=𝐊^1𝜽`$; $`𝜽`$ is initialized to the SVM solution $`𝜽^{}`$, since then the corresponding $`\stackrel{\mathbf{~}}{𝜽}`$ is obtained trivially as $`\stackrel{~}{\theta }_i=y_i\alpha _i`$ without requiring matrix inversions. The Hamiltonian equations (22,23) simplify to

$`{\displaystyle \frac{d\stackrel{~}{\theta }_i}{d\tau }}`$ $`=`$ $`p_i`$
$`{\displaystyle \frac{dp_i}{d\tau }}`$ $`=`$ $`\stackrel{~}{\theta }_i{\displaystyle \frac{V(𝜽)}{\theta _i}}`$
and the simple form of the first equation is in fact what motivated our choice of the momentum-dependent part of $`H`$, eq. (21). The correspondence between $`\stackrel{\mathbf{~}}{𝜽}`$ and $`𝜽`$ is maintained by updating $`𝜽=𝐊\stackrel{\mathbf{~}}{𝜽}`$ whenever $`\stackrel{\mathbf{~}}{𝜽}`$ is changed. As a by-product, we automatically obtain samples of $`𝐊^1𝜽`$ as required for (20).

Averages over the posterior distribution are taken by sampling after each trajectory step, repeating the procedure over some large number of steps. In practice usually the first half of the steps are discarded to allow for equilibration. We chose a total of 40,000 samples, giving 20,000 “production samples” with which to calculate the averages needed for the calculation of the gradients, eqs. (18) and (20).

### 4.3 Gradient ascent algorithm

The numerical values for the gradient of the evidence, estimated as explained above, were used in a simple gradient ascent algorithm to move the hyperparameters to a local maximum of the evidence. More powerful optimization techniques are not feasible in our case because neither the evidence itself, nor the Hessian of the evidence are available. The conjugate gradient method, for example, incorporates information about the Hessian and also usually employs a line search using the values of the function to be optimized. Approximations to Newton’s method such as the Levenberg-Marquardt algorithm also use second derivative information. Fortunately, using the first derivatives of the evidence with respect to the hyperparameters leads to convergence to some maximum in a reasonable amount of time: typically between 40 and 80 steps of gradient ascent are required before the gradients have shrunk to small values.

For the experiments described the “learning rate” multiplier for the derivative of each parameter is adapted separately throughout the optimization. This is necessary as the gradients vary over several orders of magnitude during a typical simulation. In our case the adaptation of the“learning rate” of the optimization must be based on the change in the gradients only rather than on the change in the evidence itself. We expect gradients to increase only at the start of a simulation, but thereafter they should decrease as the parameters approach a maximum in the evidence. If the gradients do not decline quickly then the learning rate is increased, if the gradients increase sharply then the ascent step is discarded and the learning rate is decreased. For vector parameters (like the length scales in an RBF kernel) the change in gradient diretion can also be used for learning rate adaptation: sudden and large changes in the gradient suggest that the optimization may have passed a maximum and the step should be redone with a smaller learning rate. As in the experiments with zero temperature Monte Carlo search, gradient ascent steps for the kernel hyperparameters were actually taken in the logarithms of these parameters.

As noted above, the HMC simulation calculates averages over the posterior with only a relatively small amount of noise. Consequently, for a given set of starting hyperparameters an optimization based on gradient ascent in the evidence is practically deterministic. So in order to investigate the properties of local maxima in the evidence repeated trials were performed with the SVM hyperparameters initialized to random values. A few preliminary trials were used to choose reasonable orders of magnitude, and unless specified otherwise all results reported begin with uniform random initalization in the ranges $`C[0.4,0.8]`$, $`\mathrm{ln}l_a[1,2]`$ for all length scales, $`\mathrm{ln}k_0[1,1]`$ and $`\mathrm{ln}k_{\mathrm{off}}[2,1]`$.

### 4.4 Computational effort

We conclude this section with a brief discussion of the computational demands of the various model selection methods; though we stress once more that our focus was not on computational efficiency, so that faster algorithms can almost certainly be designed for all of the model selection criteria that we consider.

The computationally cheapest of the simple model selection criteria is $`ϵ_{\mathrm{gacv}}`$, which can be evaluated in time $`𝒪(n)`$ from the properties of the trained SVM classifier. The span estimate $`ϵ_{\mathrm{span}}`$ requires the inversion of a matrix of size equal to the number of SVs; assuming that the number of SVs is some finite fraction of $`n`$ for large $`n`$ this gives a cost of $`𝒪(n^3)`$ for large $`n`$. The Laplace approximations to the evidence, for both linear and quadratic penalty SVMs, are dominated by the evaluation of determinants whose size is also the number of SVs (or, for linear penalty SVMs, the number of marginal SVs), giving again a scaling of approximately $`𝒪(n^3)`$.

Table 2 lists the running times for a single optimization step with each of the different methods. The evidence approximations and the test error approximations showed more or less similar running times, although the span estimate took somewhat longer on average. By far the greatest run time was needed for the gradient ascent on the evidence, due to the HMC sampling involved; a typical optimization run on a single processor HP V-Class took anywhere from 6 hours to 6 days. In comparison, most optimizations based on the simple model criteria were under an hour. The run time of the HMC algorithm should scale relatively benignly as $`O(n^2)`$ in the size of the training set, but our experiments show that the prefactor is large. The $`n^2`$ scaling comes mainly from the conversion from $`\stackrel{\mathbf{~}}{𝜽}`$ to $`𝜽`$ via $`𝜽=𝐊\stackrel{\mathbf{~}}{𝜽}`$ which is necessary during the solution of the Hamiltonian equations. (The length of the Hamiltonian trajectory, i.e. the time $`\tau `$, does not need to be increased with $`n`$; the same is true for the number of samples required to obtain the posterior averages to a given accuracy.)

Note that the theoretical dependence of running time on training set size was not strictly followed in reality. One reason for this is that the average time per step presented includes time spent on discarded steps in the zero temperature Monte Carlo search algorithms. That is, the speed of the simple optimization techniques used here depends on the complexity of the search space.

As stated above, we were interested in the evidence gradient ascent algorithm mainly as a baseline for SVM model selection based on probabilistic criteria. Computational efficiency could however be increased in a number of ways; the Nyström method , for example, could significantly reduce the dimensionality (currently $`n`$) of the space over which the posterior needs to be sampled using HMC.

## 5 Numerical results

### 5.1 Data sets

The model selection methods under consideration were applied to five two-class classification problems that are common in the machine learning literature. Three of these are from real world problems: the Pima Indian Diabetes data set, the Crabs data set and the Wisconsin Diagnostic Breast Cancer (WDBC) data set. The remaining two data sets, Twonorm and Ringnorm, are synthetic. The dimensionality of the inputs $`x`$ and the size of the training and test sets for each data set are given in Table 1. All benchmark data sets are available through the UCI Machine Learning Repository (http://www1.ics.uci.edu/$``$mlearn/MLRepository.html) and/or the DELVE archive (http://www.cs.toronto.edu/$``$delve/). More detailed descriptions are also available on the web. Inputs were standardized so that across each complete data set all input components had zero mean and unit variance. For each data set the training and test sets were held constant for all experiments. The first $`n`$ points in the data set were used for training and the remaining points were used for testing. The one exception is the Crabs data set, where the 6th attribute (color) was not used for classification and the remaining points were sampled to ensure an even distribution of the unused color attribute in the training and test sets. The number of training points, given in Table 1, was the same as that used in previous research.

### 5.2 Model selection using simple criteria

We discuss first the results obtained by optimizing the four simple model selection criteria: the Laplace evidence (LE), eq. (15), and the GACV (13) for linear penalty SVMs, and the Laplace evidence (17) and span error estimate (12) for quadratic penalty SVMs. The experiments with gradient ascent on the evidence, for linear penalty SVMs, are described separately in Sec. 5.3.

A typical example of selecting parameters for a linear penalty SVM by optimizing the Laplace approximation of the evidence is shown in Fig. 1. The data set for the experiment shown is the Twonorm data set. This example is chosen because it shows several typical features that appear in similar forms in all of the optimizations. To what extent optimizations for the other data sets match this example will be noted where appropriate.

The parameters all move more or less stochastically to stable final values as the evidence is optimized and it is clear that maximizing the Laplace Evidence correlates to reducing the error on a test set. Both the Laplace evidence and the GACV are shown alongside the test error, although only the Laplace evidence is used for optimization. Maximizing the evidence generally reduces the GACV, although this correspondence is not strict. Similar behaviour is observed for optimization with the Laplace evidence for the quadratic penalty SVM, and for optimization of the GACV and the span estimate. For quadratic penalty SVMs, the Laplace evidence and the span estimate also have the same qualitative correlation as the Laplace evidence and the GACV for the linear penalty case.

An important issue in all of these methods is the existence of many local optima in the model selection criteria. Starting from the same initialization, the hyperparameters converged to significantly different values in repeated trials. We verified explicitly, e.g. by evaluating the chosen model selection criterion along a line in hyperparameter space connecting different end points of two optimization runs, that the different local optima found were genuine and not artefacts due to incomplete convergence of the optimization algorithms. The search criteria always deteriorated in between the points found by the search, confirming that the latter were in fact local optima.

To analyse the characteristics of the local optima, 25 repeated trials were performed on all data sets for the simple optimization criteria. Comparison of the final SVM hyperparameter values at the local optima showed that they were highly variable. For all methods and all data sets the variance of the final parameter values was always of the same order of magnitude as the average value of the final parameters. Tuning of the length scales is often interpreted as “relevance determination” for the different dimensions of the data because a large length scale indicates that the classification does not vary significantly with changes in that parameter. The results here however indicate that the relevance of each dimension probably depends in a complicated way on the relevance assigned to the other dimensions, and that different assignments of the length scales can yield similar results.

In addition to variance in the final SVM hyperparameter values, the test error also showed significant trial to trial variation. For all methods, many of the trials result in a final test error that is close to the best achieved by any method, but for some methods a large portion of the trials end in a test error that is significantly worse. The average and standard deviation of the test errors achieved with the different methods are shown in Table 3. Table 4 shows the best test errors achieved on any trial for each method and data set. For reference, Table 5 shows the test errors achieved on the same data sets by comparable methods in previous research. Unfortunately, these previous studies do not always include error bars for test error results so it is hard to compare the results for the averages and standard deviations of the error.

To illustrate the variability in the final error resulting from each optimization method histograms of the errors achieved in all trials are shown for the Twonorm, Pima and WDBC data sets in Figs. 23, and 4 respectively. These plots show the difficulty of picking a “best” method from among the simple model selection criteria. For the Twonorm data set all of the methods produce test errors that are around the best for any method in previous research, but with the evidence approximation for linear penalty SVMs around a third of the trials end in errors that are signficantly greater. For the Pima data set, all of the simple methods are inferior to the best methods in previous research, and both of the evidence approximations perform worse than the error estimates; while on the WDBC data set the evidence approximations are superior to the error estimates.

Comparing Tables 3, 4, and 5, it is clear that the high variability in the results achieved by optimizing the four “simple” model selection criteria is undesirable; while the best trials for each method and data set are approximately the same as the best results reported in previous research, the average performance over trials is rather disappointing. One possible productive use of the high variability of classifiers produced by convergence to local optima of the model selection criteria could be to combine the resulting classifiers in some ensemble or voting scheme. Such approaches normally benefit precisely from high variability among the classifiers being combined, so this could be an interesting subject for future research.

### 5.3 Model selection using evidence gradients

Figs. 5 and 6 show a typical run of evidence gradient ascent on the Twonorm data set. Fig. 5 shows the tuning of a subset of the RBF kernel length scales and Fig. 6 shows the tuning of kernel amplitude $`k_0`$, the kernel offset $`k_{\mathrm{off}}`$ and the penalty parameter $`C`$. (Although statistics for the performance of the evidence gradient method were determined by initialization to random parameter values, for the specific sample shown we started all length scales with identical parameters.) Both the gradients of the evidence with respect to each parameter and the parameter values themselves are shown. The gradients typically start at small values, rise to a peak and then decline. Most parameter ultimately arrive at a constant value with small gradients, indicating that the evidence is at a local maximum with respect to that parameter. The optimization is terminated when the gradients have reached a small fraction of their peak magnitude. During this process the error on the test set decreases significantly.

As with the simple model selection criteria analysed in the previous section (Laplace approximations to the evidence and error approximations), repeated trials of gradient ascent in the evidence showed the existence of many local maxima in the evidence at widely varying parameter values. Due to the long run time required for the HMC sampling used to calculate the evidence gradients only 10 trials were performed for each benchmark data set, with the exception of the Crabs data set where 25 trials were performed. (See Section 4.4 for a discussion of the running time of the algorithm on the various data sets.) As explained in Sec. 4.3, the gradient ascent algorithm is essentially deterministic once the initial hyperparameter values have been fixed. Consequently, repeated trials were started from random initial values for the hyperparameters in order to investigate the existence and variability of local maxima in the evidence.

Tables 3 and 4 above list the resulting test errors obtained with gradient ascent optimization of the evidence, along with results obtained from the simpler methods discussed earlier. Comparing with the results found in previous studies (Table 5), one sees that gradient ascent on the evidence for SVMs with radial basis function kernels achieves approximately the same test error as the best methods that have been previously applied. For some data sets the best performance obtained by evidence gradient ascent is superior to the performance previously reported. An interesting point of comparison is with SVM model selection by optimization of a variational approximation of the evidence, as described in . (This comparison is somewhat tentative because of the small number of trials for our evidence gradient ascent method, combined with the lack of information about trial-to-trial variation in .) Still, it is worth noting that although the method desribed here uses gradients of the evidence without further approximation, and also tunes the $`C`$ parameter (which is effectively fixed to unity in the approach of ), it does not seem to achieve systematically superior performance than the variational approximation.

Figs. 2, 3 and 4 above contain the histograms of the test error produced by SVM model selection by evidence gradient ascent, and the comparison with the simple model selection criteria, for the Twonorm, Pima and WDBC data sets. Although strong conclusions cannot be drawn due to the small number of trials, gradient ascent in the evidence seems to produce significantly better performance on all of the data sets than any of the other methods. In all tests the distribution of resulting errors is both closer to the best results found in previous studies, and less variable. The most likely explanation for the superior performance of the gradient ascent method is the fact that it actually maximizes an exact, unapproximated model selection criterion (the evidence), while the simple model selection criteria (Laplace evidence and error estimates) are all to some extent approximate. The poorly performing local optima of these simple criteria may then arise from errors introduced by the approximations.

We comment briefly on the actual values of the hyperparameters found by gradient ascent on the evidence, in particular for the kernel amplitude $`k_0`$ and the offset $`k_{\mathrm{off}}`$. For the Twonorm data set, the maximum in the evidence occurs at a relatively large value of $`k_0`$, with an average of $`k_0180`$ across trials. (Large values of $`k_0`$ were also found with model selection with the approximate evidence, but were much less common when using the error estimates.) This may appear surprising. However, it should be born in mind that from the probabilistic view the prior variance of the latent function $`\theta `$ is $`\theta ^2(x)=K(x,x)=k_0+k_{\mathrm{off}}`$. The typical prior scale for $`\theta (x)`$ is therefore $`\sqrt{k_0}`$ (since $`k_{\mathrm{off}}`$ is small, see below), which equates to around 13 for $`k_0180`$; this is not unreasonably large compared to the scale of 1 set by the SVM margin. Similar final values of $`k_0`$ were obtained for the Crabs and WDBC data sets, while for Pima and Ringnorm $`k_0`$ was rather smaller. Previous experiments with simple synthetic data sets  suggest that an evidence maximum at large $`k_0`$ correlates with small apparent levels of noise in the data set; we have not attempted to verify this correlation for our five benchmark data sets.

The offset hyperparameter $`k_{\mathrm{off}}`$ was typically tuned to very small values by evidence gradient ascent (e.g. around 0.03 for the Twonorm data set). This provides a posteriori justification for our approach of including the offset parameter $`b`$ from the conventional SVM framework into the kernel.

#### 5.3.1 Noise in evidence gradients

In the final portions of the optimization shown in Figs. 5 and 6 it can be observed that there is significant noise in the gradients as the evidence approaches a maximum. This arises from statistical fluctuations in the HMC sampling, which come to dominate when the true gradient values are small. Although the noise could be decreased by increasing the length of the HMC runs, that did not seem to be necessary for the cases considered here: because of the learning rate adaptation the learning rate is quite small by the time the evidence is close to its maximum and the noise in the gradients has little effect on the final results.

In the Twonorm example it can also be seen that when the parameters are nearly at a maximum in the evidence the gradients with respect to the kernel parameters are calculated as zero in some steps. This effect occurred typically for larger values of $`C`$. Regions of $`𝜽`$-space where the potential $`V(𝜽)`$ in the Hamiltonian (21) is zero, i.e. where all $`y_i\theta _i1`$, are then much more probable then regions where $`y_i\theta _i<1`$ for some $`i`$. It is then possible that the HMC sampling only returns samples from the region with $`V(𝜽)=0`$, where $`l_p^{}(y_i\theta _i)=0`$ for all $`i`$ so that (20) gives an estimate of zero for all gradients with respect to kernel parameters. Experiments showed that scaling the trajectory length in the HMC runs proportionally to $`\frac{1}{C}`$ for large $`C`$ could avoid this effect. The rationale is that the shorter trajectories makes the HMC sampling more likely to sample values of $`𝜽`$ which are just outside the boundary of the $`V(𝜽)=0`$ region; these still have appreciable posterior probability but do give nonzero values for some of the $`l_p^{}(y_i\theta _i)`$. We did not explore this issue in detail, however.

#### 5.3.2 “Overfitting” by evidence maximization

Close inspection of progress of the test error in Figs. 5 and 6 shows an interesting aspect of tuning SVM hyperparameters using the evidence. While the overall evolution of the test error shows a large decline as gradient ascent on the evidence progresses, a closer look at the region of small error values (see the lower right plot of Fig. 6) shows that the test error goes through a shallow minimum before a small rise to its final value. Not all data sets show such a clean example of this behaviour as the Twonorm data set, but all except Pima did exhibit the phenomenon to some degree.

One possible explanation for the observed test error minimum is the fact that we are not using the evidence of a properly normalized probabilty model (see Sec. 2). An alternative interpretation, which seems to us more likely, is that we are observing here a kind of overfitting. This takes place not on the level of the “network” parameters ($`𝐰`$ or $`\theta (x)`$) as in conventional overfitting – which is due to a lack of regularization – but on the level of the hyperparameters: If we imagine sampling a number of data sets of size $`n`$ from a given true distribution, then the evidence as a function of the hyperparameters, and hence the position of its maximum, will depend on the particular data set. Only for large $`n`$ would the evidence become independent of the data set (and related to the Kullback-Liebler divergence, or cross-entropy, between the true distribution over data sets and the one predicted by the inference model; see e.g. ). For finite $`n`$, maximization of the evidence for a specific data set is therefore not expected to lead to strict minimization of the error on an independent test set.

This interpretation leads naturally to the idea of using an early stopping mechanism when optimizing the evidence, where the gradient ascent is abandoned when performance on an independent validation set ceases to improve. Note that this is not the same as simply returning to hyperparameter tuning by cross-validation; in fact, a grid search using cross-validation error over the large number of hyperparameters in our examples ($`C`$, $`k_0`$, $`k_{\mathrm{off}}`$ and the length scales $`l_a`$ associated with each of the $`d`$ input dimensions) would be essentially impossible (see also ). To gauge the possible benefits of such an approach, we have included in Table 3 above both the final test error when the optimization is run until the gradients are small, and the minimal value of the test error during the gradient ascent. True early stopping with an independent validation set would be expected to yield a performance in between these two values; the results in Table 3 suggest that this could be useful for some data sets.

## 6 Conclusion

In this paper we have investigated the issue of model selection for SVM classifiers. We have restricted ourselves to model selection in the sense of tuning the parameters of an RBF kernel and the penalty parameter $`C`$, though the general approaches described could also be used for choosing between different functional forms of the kernel.

We reviewed briefly the probabilistic view of SVMs, and extended our previous work on Laplace approximations to the evidence to the case of SVMs with quadratic slack penalties. Exact expressions for the gradients of the evidence in terms of posterior averages were also derived, and we described how these averages can be estimated numerically using Hybrid Monte Carlo techniques and used in a model selection algorithm which performs gradient ascent on the exact (unapproximated) evidence.

In our numerical experiments on five benchmark data sets, we compared optimization of four “simple” model criteria with the evidence gradient descent. Two of the simple criteria were estimates of test error: the generalized approximate cross-validation error (GACV) for SVMs with linear slack penalties, and the span error estimate for SVMs with quadratic penalties. The two other criteria were derived from probabilistic concepts; these were the Laplace approximations to the evidence for the linear and quadratic penalty cases. Our main result is that all the simple model criteria exhibit multiple local optima with respect to the hyperparameters. While some of the resulting “locally optimal” SVM classifiers give test performance that is competitive with results from other approaches in the literature, a significant fraction lead to rather higher test errors. The results for the evidence gradient ascent method show that also the exact evidence exhibits local optima. But these give much less variable test errors, which are also typically lower test errors than for the simpler model selection criteria. In this sense, “you get what you pay for”: the computationally rather more expensive evidence gradient ascent approach gives better and more consistent performance than the cheaper model selection criteria.

There are a number of directions for possible future work. First, our results strongly suggest that the hunt is still on for a model selection criteria for SVM classification which is both simple and gives consistent generalization performance. Alternatively, one could try to cope with the existence of local maxima in the simple model selection criteria by testing the selected models on a validation set and performing repeated optimizations until satisfactory performance is found. A more interesting approach might be to try to exploit the large variability in the locally optimal classifiers by using them in some scheme for combining classifiers. Finally, if evidence gradient ascent turned out in more comprehensive tests to be the model selection method of choice, it would be worth investigating possible speed-ups of the algorithm. We already hinted at the Nyström method  above, but one could also explore running the model selection only on randomly sampled subsets of data, and then possibly combining the resulting classifiers appropriately.

#### Acknowledgements

Access to the Hewlett-Packard V2500 was provided by the Caltech Center for Advanced Computing Research (http://www.cacr.caltech.edu) through the National Partnership for Advanced Computational Infrastructure–A Distributed Laboratory for Computational Science and Engineering, supported by the NSF cooperative agreement ACI-9619020. We also thank Andrew Buchan for assistance with the early stages of the numerical experiments for quadratic penalty SVMs.