An Empirical Approach For Probing the Definiteness of Kernels**preprint submitted to Soft Computing journal, Springer (2024)

2.1 Distance Measures

Distance measures compute the dissimilarity of two objects $x,x^{\prime}\in\mathcal{X}$,where we do not assume anything about the nonempty set $\mathcal{X}$.Such objects can be, e.g., permutations, trees, strings, or vectors of real values.Thus, a distance measure${d}:~{}\mathcal{X}\times\mathcal{X}\to\mathbb{R}_{+}$expresses a scalar, numerical value ${d}(x,x^{\prime})$ that should becomelarger the more distinct the objects $x$ and $x^{\prime}$ are.For a set of $n\in\mathbb{N}$ objects, the distance matrix $D$ of dimension $n\times n$collects all pairwise distances $D_{ij}={d}(x_{i},x_{j})$ with $i=1,\ldots,n$ and $j=1,\ldots,n$.Intuitively, distances can be expected to satisfy certain conditions, e.g., theyshould be zero when comparing identical objects.

A more formal definition is implied by the term distance metric.A distance metric ${d}(x,x^{\prime})$ is symmetric ${d}(x,x^{\prime})={d}(x^{\prime},x)$, non-negative ${d}(x,x^{\prime})\geq 0,$preserves identity ${d}(x,x^{\prime})=0\iff x=x^{\prime},$ and satisfies the triangle inequality${d}(x,x^{\prime\prime})\leq{d}(x,x^{\prime})+{d}(x^{\prime},x^{\prime\prime}).$Distance measures that do not preserve identity are often called pseudo-metrics.

An important class of distance measures are edit distance measures.Edit distance measures can be defined to count the minimal number of edit operationsrequired to transform one object into another. An edit distancemeasure may concern one specific edit operation (e.g., only swaps)or a set of different operations (e.g., Levenshtein distance with substitutions, deletions, insertions).Edit distances usuallysatisfy the metric axioms.

2.2 Kernels

2.3 Definiteness of Matrices

One important property of kernels and distance measures is their definiteness.We refer to the literature for more in-detail descriptions and proofs, which arethe basis for the following paragraphs(Berg etal., 1984; Schölkopf, 2001; Camastra and Vinciarelli, 2008).

2.4 Definiteness of Kernel Functions

An important special case are conditionally definite functions.Analogous to the matrix case, they observe the respective condition in equation(2).One example of a CNSD function is the Euclidean distance.The importance of CNSD functions is also due to the fact that the distance measure${d}(x,x^{\prime})$ is CNSDif and only if the kernel $k(x,x^{\prime})=\exp(-\theta{d}(x,x^{\prime}))$ is PSD $\forall~{}\theta>0$ (Schölkopf, 2001, Proposition 2.28).

It is often stated that kernels must be PSD (Rasmussen and Williams, 2006; Curriero, 2006).Still, even an indefinite kernel function may yield a PSD kernel matrix.This depends on the specific data set used to train the model(Burges, 1998);(Li and Jiang, 2004, Theorem 2)as well as the parameters of the kernel function.Some frequently used kernels are known to be indefinite.Examples are the sigmoid kernel(Smola etal., 2000; Camps-Valls etal., 2004) or time-warp kernels for time series(Marteau and Gibet, 2014).

Between these three approaches, there is a lack of straightforward empirical procedures,without resorting to complex theoretical reasoning.How can the definiteness of a function be determined? And what impact doesa lack of definiteness have on a model?

Hence, we propose the two related empirical approaches ${\text{A}_{1}}$ and ${\text{A}_{2}}$ introduced inSec.1 to fill the gap.An empirical approach may help to overcome the difficulty of theoretical considerationsor designed kernels.Empirical results may also be a starting point for a more formal approach.Furthermore, it may give a quick answer on whether or not the algorithm will have to be adaptedfor non-PSD matrices(which, if applied by default, would require additional computational effort and may limit the accuracy of derived models).

3.1 Estimating Definiteness with Random Sampling

To estimate the definiteness of a distance or kernel function we propose a simple random sampling approach (${\text{A}_{1}}$).This approach randomly generates $t\in\mathbb{N}$ sets $X_{1},\ldots,X_{t}$.Each set has size $n$, that is, it contains $n\in\mathbb{N}$ candidate samples $X=\{x_{1},\ldots,x_{n}\}$.

For each set, the distance matrix $D$ is computed, containing distances between all candidates in the set.Based on this, $\hat{D}$ is derived from equation(3).Then, the largest eigenvalue $\lambda_{n}$ of $\hat{D}$ is computed.This eigenvalue determines whether $\hat{D}$ is NSD, and hence whether $D$ is CNSD and $K$ PSD.This is repeated for all $t$ sets.The number of times that the largest eigenvalue is positive ($\lambda_{n}>0$) is retained as$n_{\lambda+}$.Accordingly, the proportion of non-CNSD matrices is determined with$p=\frac{n_{\lambda+}}{t}.$

Obviously, all distance measures that yield $p>0$ are proven to be non-CNSD.Hence, an exponential kernel based on these measures is also proven to be indefinite.If $p=0$, CNSDness is not proven or disproven.

In general, the proposed method can be categorized as a randomized algorithm of the complexity class RP (Motwani and Raghavan, 1995, p.21f.).That is, it stops after polynomially many steps, and if the output is “no” then the distance measure is non-CNSDwith probability 1, and if the output is “yes” then the distance measure is CNSD with some probabilitystrictly bounded from zero.

The parameter $p$ is an estimator of how likely a non-CNSD matrix is to occur, for the specified set size $n$.To determine definiteness, the calculation of $\lambda_{n}$ of $\hat{D}$ is not mandatory, but it may be usefulto see how close to zero $\lambda_{n}$ is, to distinguish between a pathological case andcases where the matrix is just barely non-CNSD.We will show in Section5.3 that $\lambda_{n}$ of $\hat{D}$ can be linked to model quality.

Note, that inaccuracies of the numerical algorithmused to compute the eigenvalues might lead to an erroneous sign of the largest eigenvalue.To deal with that, one could try to use exact or symbolic methods or else usea tolerance when checking whether the largest eigenvalue is larger than zero. In thelatter case, a matrix $\hat{D}$ is assumed to be non-NSD if $\lambda_{n}>\epsilon$, where $\epsilon$is a small positive number.

3.2 Estimating Definiteness with Directed Search

If very few sets $X$ yield indefinite matrices, ${\text{A}_{1}}$ may failto find indefinite matrices by pure chance.In such cases, it may be more efficient to replace random sampling with a directed search (${\text{A}_{2}}$).In detail, a set $X$ can itself be viewed as a candidate solution of an optimization problem.The largest eigenvalue $\lambda_{n}$ of the transformed distance matrix $\hat{D}$ is the objective to be maximized.By maximizing the largest eigenvalue, a positive $\lambda_{n}$ may be foundmore quickly (and more reliably).

This optimization problem is strongly dependent onthe kind of solution representation used. Evolutionary Algorithms (EAs)are a good choice to solve this problem, because they are applicableto a wide range of solution representations (e.g., real values, strings, permutations, graphs, and mixed search spaces).EAs use principles derived from natural evolution for the purpose of optimization.That is, EAs are optimization algorithms based on a cyclic repetition ofparent selection, recombination, mutation, and fitness selection (see e.g.,Eiben and Smith, (2003)).These operations can be adapted to a large variety of representations, mainly dependingon suitable mutation and recombination operators. The EA in this study will operate as follows:

4.1 Test Case: Permutations

As a proof-of-concept, we selected 16 distance measuresfor permutations; see Table1 for a complete list.The implementation of these distance measures is taken from the R-package CEGO¹¹1The package CEGO is available on CRAN at http://cran.r-project.org/package=CEGO..

Similarly to Schiavinotto and Stützle, (2007) we define$\Pi^{m}$ as the set of all permutations of the numbers $\{1,2,\ldots,m\}$.A permutation has exactly $m$ elements.We denote a single permutation with $\pi\in\Pi^{m}$ and $\pi=\{\pi_{1},\pi_{2},\ldots,\pi_{m}\}$ where $\pi_{i}$ is aspecific element of the permutation at position $i$.For example, a permutation in this notation is $\pi=\{3,2,1,4,5\}\in\Pi^{5}$.Explanations and formulas (where applicable) for the distance measures are given in appendix A.

4.2 Random Sampling

In a single experiment,$t=10,000$ sets of permutations are randomly generated.Each set contains $n$ permutations and each permutation has $m$ elements.For each set, the largest eigenvalue $\lambda_{n}$ of $\hat{D}$ is computed based on equation(3).To summarize all $t$ sets, the largest $\lambda_{n}$ as well as the ratio $p=\frac{n_{\lambda+}}{t}$ are recorded.The tolerance value used to check whether the largest eigenvalue is positive is $\epsilon$=1e-10.This process is repeated 10 times, to achieve a reliable estimate of the recorded values.

Two batches of experiments are performed. In the first, all 16 distance measures areexamined, with $n=\{4,\ldots,20\}$ and $m=\{4,\ldots,15\}$.In the second batch, larger sizes $n=\{21,\ldots,40,45,50,60,70,80,90,100\}$ are tested,but the permutations are restricted to $m=\{5,\ldots,15\}$ and the distance measures are only LCStr, Insert, Chebyshev, Levenshtein and Interchange.

4.3 Directed Search

To be comparable to the random sampling approach, the budget for each EA run is$10,000$ fitness function evaluations. A run will stop if the budget isexhausted or if $\lambda_{n}>\epsilon=10^{-10}$.The population size of the EA is set to 100. The recombination rate is 0.5, the mutationrate is $1/m$ and truncation selection is used.

To identify bias introducedby the choice of the submutation operator (which may have a strong interaction withthe respective distance measures), each EA run is performed repeatedly withthree different submutation operators:

All 16 distance measures aretested, with
$n=\{4,\ldots,20\}$ and $m=\{4,\ldots,15\}$. With ten repeats, and the three differentsubmutation operators, this results into $97,920$ EA runs, each with $10,000$ fitness function evaluations.The employed EA implementation is part of the R-package CEGO.

4.4 Tests For Other Search Domains

To show that the proposed approach is not limited to the presented permutation distance example, wealso briefly explore other search domains and their respective distances. However, these are not analyzed in further detail. Instead,we provide a list of minimal examples in Appendix B: the smallest (w.r.t.dimension) indefinite distance matrix for eachtested distance measure. The examples include distance measures for permutations, signed permutations,trees, and strings.

5.1 Sampling Results

The proportions of sets with positive eigenvalues ($p$) are summarized in Fig.1.The largest eigenvalues are depicted in Fig.2.Only the five indefinite distance measures, which achieved positive eigenvalues are shown:

Longest Common Substring, Insert, Chebyshev, Levenshtein, Interchange.

No counter-examples are found for the remaining eleven measures. That does notprove that they are CNSD (although some are CNSD, e.g., Euclidean distance, Swapdistance(Jiao and Vert, 2015), Hamming distance(Hutter etal., 2011)),but indicates that it may be unproblematic to use them in practice.Some of the five non-CNSD distance measures were reported to work wellin a surrogate-model optimization framework, e.g., Levenshtein distanceseemed to work well for modeling of scheduling problems(Zaefferer etal., 2014a ).This is mainly due to the fact, that even a non-CNSD distance may yieldPSD kernel matrix $K$, depending on the specific data set and kernel parameters used.We do not suggest that non-CNSD distance measures should be avoided, but thattheir application should be handled with care.

Regarding values of $p$ (Fig.1),some trends can be observed.For indefinite distance measures, increasing theset size($n$) will in general lead to larger values of $p$. Obviously, a larger set is more likely to contain combinations of samples that yieldnegative eigenvalues. In addition, the lower bound for eigenvalues decreases with increasingmatrix size(Constantine, 1985).

In contrast to the set size, increasing the number of permutation elements ($m$) decreases the proportion of positive eigenvalues $p$ in all five cases. This can be attributed to a larger and hence more difficult search space.Overall, none of the distance measures shows exactly the same behavior. LCStr distancehas the least problematic behavior. Only very few sets (of comparatively large size) yielded positive eigenvalues with LCStr distance.Interchange, Levenshtein and Insert distance all have relatively large $p$ for small set sizes $n$. Chebyshev on the other hand,starts to have non-zero $p$ for relatively large $n$. However, the number of permutation elements has only a weak influencein case of Chebyshev distance. Hence, curves for different $m$ are much closer to each other, compared to the other distance measures.Somewhat analogous to $p$, the largest eigenvalues of $\hat{D}$ plotted in Fig.2 are generally increasing for larger $n$,and decreasing with larger $m$.

Our findings are confirmed by some results from literature.Cortes etal., (2004) have shownthat an exponential kernel function based on the Levenshtein distancesis indefinite for strings of more than one symbol. Our experiments show that this result can be easilyrediscovered empirically, for the case of permutations.At the same time, these findings also confirm (and provide reasons for)problems observed with these kinds of kernel functions in a previous study(Zaefferer etal., 2014a ).

As a consistency check, Table2 compares the samplingresults to a brute force approach for $n=\{4,\ldots,8\}$ and $m=4$.It shows that the sampling indeed approximatesthe number of non-CNSD matrices quite well. In this small set, the sampling identified all combinations of distance measure and $n$ that may yield non-CNSD matrices.

5.2 Optimization Results

The average number of fitness function evaluations (i.e., number of times $\lambda_{n}$ of $\hat{D}$ is computed)required to find a matrix $\hat{D}$ with positive $\lambda_{n}$ is depicted in Fig.3.

The optimization results show similar behavior with respect to $n$ and $m$ as the sampling approach. Increasing $m$leads to an increased number of fitness function evaluations. That means, finding positive eigenvalues becomes more difficult with increasingvalues of $m$.Increasing $n$ reduces the number of required fitness function evaluations. That is, finding positive eigenvalues becomes easier.In some cases, this effect disappears for large values of $m$, e.g., for LCStr distance, where the average is more or less constant over$n$, if $m$ is large enough.

Importantly, the comparison to the sampling results clearly shows that the EA has some success in optimizing the largest eigenvalueof the transformed distance matrix $\hat{D}$.In several cases, positive eigenvalues are found by the EA while sampling with the same budgetfailed to find any.Hence, it can be assumed that the fitness landscape based on $\lambda_{n}$ of $\hat{D}$ is sufficiently smooth to allow for optimization.The eigenvalue $\lambda_{n}$ seems to be a good indicator of how close a solution set is to yielding an indefinite kernel matrix.

Furthermore, the expected bias of the used submutation operator becomes visible.For Insert and LCStr distance, the EA with swap mutation works considerably better than the other two variants.Hence, comparisons of these values across different distance measures should be handled with caution.Clearly, other aspects of the optimization algorithm (e.g., selection criteria or recombination operators) might have similar effects.While this bias is troubling, results may still offer interesting insights. The eigenvalue optimizationcould be interpreted as a worst-case scenario that occurs if an iterative learning process strongly correlates with theeigenvalues of the employed distance matrix.

5.3 Verification: Impact on Model Quality

Earlier, we discussed two valuesthat may express the effect of the lack of definitenessin practice, i.e., $p$ and $\lambda_{n}$ of $\hat{D}$.But what do these values imply?

The value $p$ can be seen as a probability of generating indefinite matrices.If we assume that a model is unable to deal with indefinite data,the fraction $p$ is an estimate of how likely a modeling failure is.For other cases, it is very hard to link it to a model performance measuresuch as accuracy without making too many assumptions.We still argue that $p$ provides useful information, especiallywhen the kernel is designed and probed before sampling any data (e.g.,when planning an experiment).We suggest to use $p$ to support an initial decision (e.g., whetherto spend additional time on fixing or otherwise dealing with the indefinite kernel).One advantage is that it is rather easy to interpret.

In contrast, the parameter $\lambda_{n}$ of $\hat{D}$ is more difficult to interpret.But it has the advantage that it may be estimated for a single matrixas well as its average for a set of matrices. Hence, we want to determinewhether the magnitude of this eigenvalue affects model performance.We expect an influence that depends on the choice of model.Consider, e.g., a Gaussian process regression model, as e.g., described by Forrester etal., (2008).The model may be able to mitigate the problematic eigenvalueby assigning larger $\theta$ values to the kernel $k(x,x^{\prime})=\exp(-\theta{d}(x,x^{\prime}))$.For very large $\lambda_{n}$ and thus very large $\theta$, this will lead to kernel matricesthat approximate the unit matrix, which is positive definite.A model with a unit kernel matrix would be able to reproducethe training data, but would predict the process mean for most other data points.Hence, we examine Gaussian process regression models, since they providea transparent and interpretable test case (but similar experimentscould easily be made with support vector regression).

An experimental test has to consider the potential bias of the used data set.We need to be able to reasonably assume that differences in performanceare actually due to the properties of employed distance or kernel(i.e., a kernel performs poorly because the corresponding $\lambda_{n}$ of $\hat{D}$ is high)rather than properties of the data set(i.e., a kernel performs poorly because it does not fit well to the groundtruth of the data set).To that end, we suggest that observations in a test data set are derived fromthe same distances that are used in the model.

Hence, we randomly created data sets $X$ of size $n$ with permutations of dimension $m$,similarly to the random sampling performed earlier.Then, we created training observations by evaluating the distance of each permutation in $X$to a reference permutation $x_{ref}={1,...,m}$, i.e., $y={d}(x,x_{ref})$.A Gaussian process model was then trained with this data, largely followingthe descriptions inForrester etal., (2008).The model was trained with the kernel $k(x,x^{\prime})=\exp(-\theta{d}(x,x^{\prime}))$,and the model parameters (e.g., $\theta$) were determined by maximum likelihood estimation, via the locally biasedversion of the DIviding RECTangles (DIRECT) algorithm (Gablonsky and Kelley, 2001)with 1,000 likelihood evaluations.For each test, the distance chosen to produce the observations $y$ and the distance chosen in the kernel were identical.

The Root MeanSquare Error (RMSE) of the model was evaluated on 1000 randomly chosen permutations.The resulting RMSE values for each training set $X$, as well as the correspondingeigenvalue $\lambda_{n}$ of $\hat{D}$ are shown in Fig.4.The graphic shows a trend that confirms our expectation. It seems that distancesassociated to larger $\lambda_{n}$ tend to produce larger errors.