Permanents, Transport Polytopes and Positive Definite Kernels on Histograms

1 Introduction

Defining meaningful kernels on histograms and clouds of points - and more generally positive measures on an arbitrary space X - is an important topic in the field of kernel methods, as it is directly related to the definition of kernels for structured objects seen as bags-of-components. Since the latter representations are frequently used by practitioners in applications, notably images seen as histograms of colors, texts as bags-of-words or sequences as groups of subsequences, research has been active in this field recently.

In the early applications of kernel methods to complex data structures, histograms were often treated as vectors, and used as such with the standard Gaussian or polynomial kernels [Joachims, 2002]. More adequate positive definite kernels which exploit their specificity have been proposed since. Namely, kernels which take into account the fact that histograms are vectors with nonnegative coordinates [Hein and Bousquet, 2005], and whose sum may be normalized to one, that is cast as discrete probability measures and treated under the light of information geometry [Lafferty and Lebanon, 2005; Lebanon, 2006]. Since such histograms are usually defined on bins which are not equally dissimilar, as is for instance the case with color or amino-acid histograms, further kernels which may take into account an a priori inter-bin similarity where subsequently proposed [Kondor and Jebara, 2003; Cuturi et al., 2005; Hein and Bousquet, 2005].

In this context, a well-known distance for probability measures on a space X which takes explicitly into account the geometry of X is the optimal transportation distance [Villani, 2001], which is usually known as the Monge-Kantorovich (MK) or Wasserstein distance. This distance is also popular in the computer vision community [Rubner et al., 2000] under the name of the earth movers' distance. However, preliminary findings [Naor and Schechtman, 2005] suggest that the MK distance is not negative definite, and cannot thus be used directly to define positive definite kernels, through the Schoenberg theorem 1 and negative exponentiation for instance. Although some approximations of the distance have been used so far to define positive definite kernels in vision applications [Grauman and Darrell, 2004], we propose in this paper to consider not only the optimal transport plan, but the whole of the transportation polytope to characterize the similarity of two histograms r and c.

This idea is rooted in the approach of [Vert et al., 2004] to define a positive definite kernel for strings derived from a set of string manipulations which may map a string m1 to another string m2, namely sequences of deletions, substitutions and insertions of tokens, known as alignments. Vert et al. [2004] consider all possible alignments π between m1 and m2 and associate to each of these alignments π a score S (π) which quantifies how efficiently the sequence π aligns successive tokens of m1 and m2. Rather than considering the score S (π⋆ ) of the optimal alignment π⋆ , known as the Smith-Waterman score in the context of biological sequences, Vert et al. [2004] propose to define a positive definite kernel betweenm1 andm2 through the sum ∑π eβS(π) , β > 0, which in their experimental setting provides a much better performance. Intuitively, the latter sum can be interpreted as the generating function of the set of all alignments π, which may give it more discriminative power than the simple use of the extremum S (π⋆ ).

1 The Schoenberg theorem[Berg et al., 1984, Theorem 3.2.2], states that if a function ψ is a negative definite kernel, equivalently that -ψ is conditionally positive definite [Sch¨olkopf and Smola, 2002], then for all t > 0, exp(-tψ) is a positive definite kernel. The Gaussian kernel is for instance based on the fact that for two vectors x, y, ‖x - y‖2 is negative definite.

In the context of this paper, the set of all alignments π is played by the set of all d × d transportation matrices F ∈ U (r, c) between two discrete histograms r = (r1, . . . , rd) and c = (c1, . . . , cd) of equal sumN; the analog of the BLOSUM distance matrices between amino-acids used in Vert et al. [2004] is an arbitrary d × d negative definite distance matrix T between the bins of the histograms, and finally, the cost π(σ) becomes simply the Frobenius dot-product < F, T >. As is also the case in [Vert et al., 2004], the family of convolution kernels introduced by Haussler [1999] plays an important role in our proofs, notably for the permanent kernel for clouds of points introduced in Section 2. We then propose in Section 3 a kernel between histograms by only taking into account the volume of U (r, c), and we show that this result is a natural consequence of the Robinson-Schensted-Knuth correspondence between transportation matrices and generalized Young tableaux. Finally, inspired by a recent construction obtained by Barvinok [2005] and through a lemma by Berg and Duran [2004], we show in Section 4 that a weighted version of the generating function of U (r, c) can be used to define positive definite kernels. We close the paper with Section 5 by discussing implementation and computational issues brought forward by these kernels, as well as preliminary experimental results led on a subset of the MNIST database of handwritten digits.

2 Permanent kernel for clouds of points

We define in this section a kernel for two clouds of points x = {x1, . . . , xn} and y = {y1, . . . , yn} in a space X endowed with a kernel κ, through a kernel on arbitrary sequence representations x = (x1, . . . , xn) and y = (y1, . . . , yn), which is by definition invariant under reordering of these terms. Recall that for a n × n matrix M = [mij], the permanent of M, per M is defined as the quantity

where σ spans the symmetric group Sn, that is the set of all permutations of {1, . . . , n}. Note that the definition of the permanent of a matrix differs from that of its determinant in that the signatures of the permutations are not taken into account. The permanent of a matrix is also invariant under any permutation of columns or rows of this matrix, a fact that we use in the proof of Proposition 1 below.

Proposition 1 Let X be a set endowed with a kernel κ and Xn the set of clouds of points of X of cardinal n, that is {X = {x1, . . . , xn}, xi ∈ X}. Let X = {x1, . . . , xn} and Y = {y1, . . . , yn} ∈ Xn. Then

is a positive definite kernel on Xn × Xn.

Proof. We first prove the result for two sequences x = (x1, . . . , xn), y = (y1, . . . , yn). Given a permutation σ we write xσ for the sequence (xσ(1), . . . , xσ(n)). The proof now follows from the work of Haussler [1999] on convolution kernels. Namely, we consider the equivalence relation R between two sequences x, y where xRy if and only if there exists a permutation σ ∈ Sn such that xσ = y. Consider now the kernel k for two sequences

which we use to define the convolution kernel K, which is positive definite by definition,

hence the positive definiteness of kper used on clouds of points X, Y represented through any arbitrary pair of sequences x, y.

Suppose for interpretation purposes that the kernel k can be written as k(x, y) = e-d(x,y) where d is an Hilbertian metric [Hein and Bousquet, 2005] on X, that is there exists a mapping φ from X to an arbitrary Hilbert space H such that d(x, y) = ‖φ(x) - φ(y)‖. In that case the permanent can be interpreted as the sum over all possible matchings σ of the weight of each total transport scheme e-dσ where dσ = ∑i d(xi, yσ(i)). The quantity dσ stands for the total transport cost between the two clouds-of-points given the transport plan σ is selected, taken in a feature space H, as illustrated in Figure 1.

Note finally that a possible way to define kernels for two clouds X and Y of sizes n and m respectively is to consider the sum of the pairwise kernels of all their respective subsets of size d ≤ min(n,m), that is

3 The volume of the transport polytope as a kernel for marginals

We write N = {0, 1, . . .} for the set of nonnegative integers, and consider now integral histograms (or margins as in the corresponding statistical literature [Diaconis and Gangolli, 1995]) with identical overall sum and dimension, that is elements of the simplex lattice

We consider the polytope of transport matrices between r and c, restricted to integral matrices, that is

where 1d is the d-dimensional vector of ones. We recall that the optimal transportation cost from r to c is a symmetric function in r and c which is defined as the result of the optimization

where T ∈ Rd,d is an arbitrary distance matrix between bins, and for two square matrices U and V we use the Frobenius dot-product < U, V >= tr(UV ⊤). Note that the optimal plan

can be computed through standard linear-programming methods in polynomial time in d and it is known that F ⋆ is a vertex of the polytope. We will reconsider the cost parameter T in the next section, and focus for the rest of this section on the volume |U (r, c)| of the polytope U (r, c), that is the total number of integral transportation plans. We introduce first the concept of semi-standard Young tableaux.

For two partitions u, v of an integer n, a semi-standard Young tableau of shape u and weight v is a diagram of shape u containing v1 ones, v2 twos, etc., arranged to be weakly right increasing in rows and strictly increasing down columns. Consider for instance for n = 15, the semi-standard tableau

6

of shape u = (6, 4, 4, 1) and weights v = (3, 2, 3, 2, 1, 3, 1). Given two partitions u and v, the Kostka number Ku,v is equal to the number of semi-standard Young tableaux of shape u and weight v.

Proposition 2 The kernel kvol on Σd,N defined as

is symmetric positive definite.

Proof. The proof can be derived either from the theory of symmetric functions or from the Robinson-Schensted-Knuth (RSK) correspondence [Knuth, 1970], with both approaches mentionned in [Diaconis and Gangolli, 1995]. We recall briefly the second proof. The RSK bijective correspondence states that to every matrix M ∈ U (r, c) corresponds one and only pair of semi-standard Young tableaux of identical shape and weights r and c respectively. We hence have that, summing over all possible partitions η of N used as shapes for the Young tableaux,

which is sufficient to prove that the volume |U (r, c)| satisfies Mercer's condition.

4 Weighted generating functions

The volume of the transportation polytope is a special case of the evaluation of the generating function f of U (r, c) [Barvinok, 2006] on a given cost matrix T , in that case the null matrix, where f is more generally defined as

The computation of f for general polytopes, notably the transportation one, is a subject of extensive research, with significative developments carried out in recent years and summarized in [Loera et al., 2004]. The generating function can also be expressed as the total weight of U (r, c) if we use the terminology of [Barvinok, 2005] by setting wij = e-tij . Note that for binary histograms, that is marginals which may either take 1 or 0 values, U (r, c) is known as the Birkhoff polytope and f (T ) corresponds in this case to the permanent kernel defined in Proposition 1 with κ set to e-T . In the general case where r and c may not have binary counts, the volume |U (r, c)| can thus be regarded as a similarity between r and c based exclusively on combinatorial properties, regardless of any prior knowledge T on the distance between the d bins. On the other hand, the optimal plan F ⋆ corresponding to the Monge-Kantorovich distance dMK =< F ⋆ , T > takes into account such an information but does not reflect the information carried out by the distribution of the costs found in the whole polytope. Hence, having in mind Figure 2 we propose to define valid kernels kϕ which consider both criterions, that is we consider the distribution of the cost < F, T > over the whole polytope U (r, c). We do so by introducing weighted versions of the generating function through a function ϕ : U (r, c) →R, defining

Note that the generating functionis recovered when ϕ = 1. Although the generating function might be a good candidate for a kernel between r and c, we do not know at this moment whether its evaluations on arbitrary matrices T are positive definite functions of r and c. We provide instead a family of functions ϕ which ensures this condition:

Proposition 3 Given that T is such that [e-tij ] is positive semidefinite, and having defined for 0 ≤ a < 2 the weight function ϕa : Nd×d → R as

we have that kϕa is a positive definite kernel on Σd,N.

The symmetry of kϕa is ensured by the symmetry of T and ϕa since U (c, r) = U (r, c)⊤ and we have that

We prove the positive-definiteness of kϕa using the following 3 lemmas, which are motivated by a recent characterization carried out by Barvinok [2005] of the generating function of the transport polytope in terms of random permanents.

Lemma 4 Let T be a d × d cost matrix such that [e-tij ]1≤i,j≤d is positive semidefinite and γ = (γ1, . . . , γd) a sequence of nonnegative real numbers. For r, c ∈ Σd,N define the N × N block matrix A as

where each block Ai,j is the ri × cj rectangular matrix with all coefficients set to the constant γiγje-ti,j . Then

is a positive definite kernel on Σd,N × Σd,N .

Proof. we first map each marginal r and c to the corresponding sequences

and ˜ c, and define the positive definite kernel

for the kernel indexed on {1, . . . , d} × {1, . . . , d}. We then have using Proposition 1 that:

Since (r, c) ↦→ 1 r1!···rd! × 1 c1!···cd! is trivially positive definite, so is kγ,T as the product of two positive definite kernels.

We use the following lemma to turn the randomized setting proposed in [Barvinok, 2005] into a sum of positive definite kernels:

Lemma 5 (Berg, Duran) For each 0 < α ≤ 2, the sequence (n!)α , n ∈ N is a determinate Stieltjes moment sequence, that is there exists a unique nonnegative measure μα on [0,∞[ such that ∫∞ 0 xnμα(x) = (n!)α for n ∈ N.

We refer to [Berg and Duran, 2004] for a proof of this result, and more generally to the reference [Berg et al., 1984] for the exposition of the moment problems and their relationship with harmonic analysis on semigroups. Note that in the case where α = 1, μ1 is the standard exponential density, and for α = 0 the measure μ0 can be simply defined as the dirac mass on 1.

Lemma 6 (Barvinok) Let 0 ≤ a < 2 and suppose γ = (γ1, . . . , γd) is distributed as a sequence of independent random variables with identical law μa. Through the identity

we have that kϕa is positive definite.

Proof. E[kγ,T(r, c)] is trivially positive definite as a sum of positive definite kernels. We follow Barvinok's proof to prove the equality, with a slight modification: Barvinok considers standard exponential variable γij arranged in a d × d matrix, while we consider here a sequence of independent random variables γ = (γ1, . . . , γd) which all follow law μa.

Let us consider matrix A defined in Lemma 4. For every permutation σ of SN let

be the corresponding term in per A. Hence

Following Barvinok, with every permutation σ we associate a transport plan D = D(σ) of U (r, c) called the pattern of σ, as follows. Namely D = [dij ]1≤i,j≤d where

that is dij is the number of indices k ∈ {1, . . . ,N} such that (k, σ(k)) is in the (i, j ) block of A. Note that D : σ ↦→ D(σ) ∈ U (r, c) is surjective, but not bijective as we see below. For hσ, we thus have, through Lemma 5 that

At this point we follow exactly Barvinok's proof. Barvinok proves that the number of permutations σ of SN which admit D as a pattern is

yielding

and hence

which concludes the proof.

Note that when a = 0 we obtain the Fisher-Yates distribution on transportation matrices [Diaconis and Gangolli, 1995], that is

whereas the case a = 12 yields a weight which only depends on the diagonal elements of the transport plan,

5 Discussion and Experiments

We discuss in this section complexity issues which may arise when trying to compute the kernels presented above, and we present preliminary results on a pattern recognition task which involves comparing clouds of points.

5.1 Computational issues

To handle clouds of points through Proposition 1 requires the computation of the permanent of a n × n matrix, which is a notoriously difficult problem in combinatorics. Millions of computations of such kernel evaluations, which are usually required to fill in Gram matrices, may not be tractable at the moment when the number n of points exceeds twenty to thirty points. However, and in the case where the kernel κ is bounded between 0 and 1, recent advances 2 in the computation of approximations of the permanent through Sequential Monte Carlo (SMC) techniques [Jerrum et al., 2004] yield a complexity of the order of n7 log4 n. This is still problematic for large n, but we believe that for clouds of points of small size the permanent might be a useful kernel, with the ability of quantifying complex relationships through the power of combinatorics. We propose below in our experiments to compare 2000 images of handwritten digits by sampling artificially 20 black pixels among each image, and compare these clouds of points through the permanent of the pairwise kernels for the points in each cloud. Another issue in that case arises from the numerical stability of the computation of the permanent when the values for κ might be too small, and we do not have an adequate answer to this problem other than simple cross-validation to obtain reasonable entries.

In the more general case of histograms, both the computation of the volume |U (r, c)| and the integration of kϕ0 , which corresponds to the Fisher-Yates distribution, as well as that of kϕ1/2 , may be computed through SMC sampling methodologies presented in recent works [Chen et al., 2006]. For the volume only, exact calculations through toolboxes such as LattE [Loera et al., 2004] are possible, but only tractable for very low dimensions. Diaconis and Gangolli [1995] propose ad-hoc numerical approximations when d is small and N is large,

where

Although these expressions might be symmetrized by averaging 12 (|U (r, c)| + |U (c, r)|), their positive definiteness may not be guaranteed and has yet to be tested on datasets.

5.2 Experiments

Following the previous work of Kondor and Jebara [2003], we have conducted experiments on the first 2000 images (28 × 28 pixels) of the MNIST database of handwritten digits, with approximately 200 images for each digit. For each image X i , we randomly sample a set {xi1 , . . . , xi20 } of 20 distinct black points in the image, that is pixels with an intensity superior to 190 represented as points of the square

2 which we have not used in our experiments.

[0, 1] × [0, 1], and perform a multiclass classification to classify any new image as one of the 10 digits. We do so by applying a simple one-vs-all strategy on 10 classifiers, namely support vector machines, using the Spider toolbox. This setting makes it particularly difficult for most common kernels to compare the images and we consider here two different approaches: first the permanent kernel of Equation (1), second a Gaussian kernel taken between the two images seen as 28 × 28 dimensional vectors, preliminarily smoothed through a smoothing-kernel κ on the pixels, which yields actually a simple summation over Gram matrices as described for instance in [Borgwardt et al., 2006]. The kernel κ used to compute both the permanent as in Equation (1), and to smooth the image in the second case was set to be the Gaussian kernel between pixels κ(x, y) = exp(-‖x - y‖2/σ2 ) with a width σ spanning values 0.1, 0.2 and 0.3. The considered kernels are thus

and we use their normalized counterpart, that is using ˜ k(x, y) = k(x, y)/√k(x, x)k(y, y) instead of k in our experiments. We report the cross-validation errors for these settings for σ taken over 5-fold cross validations in Table 1, which show that the permanent kernels appear as a robust although costly alternative to the smoothed kernel in this preliminary experiment.