Thème Sym

Nielsen -Jean-Daniel Boissonnat

Richard Nock

Frank Nielsen

Jean-Daniel Boissonnat

Bregman Voronoi Diagrams: Properties, Algorithms and Applications

Keywords: Computational Information Geometry, Voronoi diagram, Delaunay triangulation, Bregman divergence, Quantification, Sampling, Clustering Diagrammes de Bregman-Voronoi : Propriétés, Algorithmes et Applications Géométrie de l'information, géométrie algorithmique, diagrammes de Voronoï, divergence de Bregman, quantification, échantillonnage, classification

The Voronoi diagram of a finite set of objects is a fundamental geometric structure that subdivides the embedding space into regions, each region consisting of the points that are closer to a given object than to the others. We may define many variants of Voronoi diagrams depending on the class of objects, the distance functions and the embedding space. In this paper, we investigate a framework for defining and building Voronoi diagrams for a broad class of distance functions called Bregman divergences. Bregman divergences include not only the traditional (squared) Euclidean distance but also various divergence measures based on entropic functions. Accordingly, Bregman Voronoi diagrams allow to define information-theoretic Voronoi diagrams in statistical parametric spaces based on the relative entropy of distributions. We define several types of Bregman diagrams, establish correspondences between those diagrams (using the Legendre transformation), and show how to compute them efficiently. We also introduce extensions of these diagrams, e.g. k-order and k-bag Bregman Voronoi diagrams, and introduce Bregman triangulations of a set of points and their connexion with Bregman Voronoi diagrams. We show that these triangulations capture many of the properties of the celebrated Delaunay triangulation. Finally, we give some applications of Bregman Voronoi diagrams which are of interest in the context of computational geometry and machine learning.

Introduction and prior work

The Voronoi diagram vor(S) of a set of n points S = {p 1 , ..., p n } of the d-dimensional Euclidean space R d is defined as the cell complex whose d-cells are the Voronoi regions {vor(p i )} i∈{1,..,n} where vor(p i ) is the set of points of R d closer to p i than to any other point of S with respect to a distance function δ:

vor(p i ) def = {x ∈ R d | δ(p i , x) ≤ δ(p j , x) ∀ p j ∈ S}.
Points {p i } i are called the Voronoi sites or Voronoi generators. Since its inception in disguise by Descartes in the 17th century [5], Voronoi diagrams have found a broad spectrum of applications in science. Computational geometers have focused at first on Euclidean Voronoi diagrams [5] by considering the case where δ(x, y) is the Euclidean distance ||x -y|| = d i=1 (x i -y i ) 2 . Voronoi diagrams have been later on defined and studied for other distance functions, most notably the L 1 distance ||x -y|| 1 = d i=1 |x i -y i | (Manhattan distance) and the L ∞ distance ||x -y|| ∞ = max i∈{1,...,d} |x i -y i | [START_REF] Boissonnat | Algorithmic Geometry[END_REF]5]. Klein further presented an abstract framework for describing and computing the fundamental structures of abstract Voronoi diagrams [START_REF] Klein | Concrete and Abstract Voronoi Diagrams[END_REF][START_REF] Boissonnat | Curved Voronoi diagrams[END_REF].

In artificial intelligence, machine learning techniques also rely on geometric concepts for building classifiers in supervised problems (e.g., linear separators, oblique decision trees, etc.) or clustering data in unsupervised settings (e.g., k-means, support vector clustering [2], etc.). However, the considered data sets S and their underlying spaces X are usually not metric spaces. The notion of distance between two elements of X needs to be replaced by a pseudo-distance that is not necessarily symmetric and may not satisfy the triangle inequality. Such a pseudo-distance is also referred to as distortion, (dis)similarity or divergence in the literature. For example, in parametric statistical spaces X , a vector point represent a distribution and its coordinates store the parameters of the associated distribution. A notion of "distance" between two such points is then needed to represent the divergence between the corresponding distributions.

Very few works have tackled an in-depth study of Voronoi diagrams and their applications for such a kind of statistical spaces. This is all the more important even for ordinary Voronoi diagrams as Euclidean point location of sites are usually observed in noisy environments (e.g., imprecise point measures in computer vision experiments), and "noise" is often modeled by means of Normal distributions (so-called "Gaussian noise"). To the best of our knowledge, statistical Voronoi diagrams have only been considered in a 4-page short paper of Onishi and Imai [START_REF] Onishi | Voronoi diagram in statistical parametric space by Kullback-Leibler divergence[END_REF] which relies on Kullback-Leibler divergence of dD multivariate normal distributions to study combinatorics of their Voronoi diagrams, and subsequently in a 2-page video paper of Sadakane et al. [START_REF] Sadakane | Voronoi diagrams by divergences with additive weights[END_REF] which defines the divergence implied by a convex function and its conjugate, and present the Voronoi diagram with flavors of information geometry [START_REF] Amari | Methods of Information Geometry[END_REF] (see also [START_REF] Onishi | Voronoi diagrams for an exponential family of probability distributions in information geometry[END_REF] and related short communications [START_REF] Inaba | Geometric clustering for multiplicative mixtures of distributions in exponential families[END_REF][START_REF] Inaba | Geometric clustering models for multimedia databases[END_REF]). Our study of Bregman Voronoi diagrams generalizes and subsumes these preliminary studies using an easier concept of divergence: Bregman divergences [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF]6] that do not rely explicitly on Vor(p 6 ) Figure 1: Ordinary Euclidean Voronoi diagram of a given set S of seven sites. In the bounded Voronoi cell vor(p 6 ), every point p ∈ vor(p 6 ) is closer to p 6 than to any other site of S (with respect to the Euclidean distance). Dashed segments denote infinite edges delimiting unbounded cells.

INRIA convex conjugates. Bregman divergences encapsulate the squared Euclidean distance and many widely used divergences, e.g. the Kullback-Leibler divergence. It should be noticed however that other divergences have been defined and studied in the context of Riemannian geometry [START_REF] Amari | Methods of Information Geometry[END_REF]. Sacrifying for some generality, while not very restrictive in practice, allows a much simpler treatment and our study of Bregman divergences is elementary and does not rely on Riemannian geometry.

In this paper, we give a thorough treatment of Bregman Voronoi diagrams which elegantly unifies the ordinary Euclidean Voronoi diagram and statistical Voronoi diagrams. Our contributions are summarized as follows: Since Bregman divergences are not symmetric, we define two types of Bregman Voronoi diagrams. One is an affine diagram with convex polyhedral cells while the other one is curved. The cells of those two diagrams are in 1-1 correspondence through the Legendre transformation. We also introduce a third-type symmetrized Bregman Voronoi diagram.

We present a simple way to compute the Bregman Voronoi diagram of a set of points by lifting the points in a higher dimensional space using an extra dimension. This mapping leads also to combinatorial bounds on the size of these diagrams. We also define weighted Bregman Voronoi diagrams and show that the class of these diagrams is identical to the class of affine (or power) diagrams. Special cases of weighted Bregman Voronoi diagrams are the k-order and k-bag Bregman Voronoi diagrams.

We define two triangulations of a set of points. The first one captures some of the most important properties of the well-known Delaunay triangulation. The second triangulation is called a geodesic Bregman triangulation since its edges are geodesic arcs. Differently from the first triangulation, this triangulation is the geometric dual of the first-type Bregman Voronoi diagram of its vertices.

We give a few applications of Bregman Voronoi diagrams which are of interest in the context of computational geometry and machine learning.

The outline of the paper is as follows: In Section 2, we define Bregman divergences and recall some of their basic properties. In Section 3, we study the geometry of Bregman spaces and characterize bisectors, balls and geodesics. Section 4 is devoted to Bregman Voronoi diagrams and Section 5 to Bregman triangulations. In Section 6, we select of few applications of interest in computational geometry and machine learning. Finally, Section 7 concludes the paper and mention further ongoing investigations.

Notations. In the whole paper, X denotes an open convex domain of R d and F : X → R a strictly convex and differentiable function. F denotes the graph of F , i.e. the set of points (x, z) ∈ X × R where z = F (x). We write x for the point (x, F (x)) ∈ F. ∇F , ∇ 2 F and ∇ -1 F denote respectively the gradient, the Hessian and the inverse gradient of F .

Bregman divergences

In this section, we recall the definition of Bregman 1 divergences and some of their main properties ( §2.1). We show that the notion of Bregman divergence encapsulates the squared Euclidean distance as well as several well-known information-theoretic divergences. We introduce the notion of dual divergences ( §2.2) and show how this comes in handy for symmetrizing Bregman divergences ( §2.3). Finally, we prove that the Kullback-Leibler divergence of distributions that belong to the exponential family of distributions can be viewed as a Bregman divergence ( §2.4).

Definition and basic properties

For any two points p and q of X ⊆ R d , the Bregman divergence 2 D F (•||•) : X → R of p to q associated to a strictly convex and differentiable function F (called the generator function of the divergence) is defined as

D F (p||q) def = F (p) -F (q) -∇F (q), p -q , (1)
where ∇F = [ ∂F ∂x1 ... ∂F ∂x d ] T denotes the gradient operator, and p, q the inner (or dot) product:

d i=1 p i q i .
Informally speaking, Bregman divergence D F is the tail of the Taylor expansion of F . See [START_REF] Csiszár | Why least squares and maximum entropy? An axiomatic approach to inference for linear inverse problems[END_REF] for an axiomatic characterization of Bregman divergences as "permissible" divergences.

Lemma 1 The Bregman divergence D F (p||q) is geometrically measured as the vertical distance between p and the hyperplane H q tangent to F at point q: D F (p||q) = F (p) -H q (p).

Proof: The tangent hyperplane to hypersurface F : z = F (x) at point q is H q : z = F (q) + ∇F (q), xq . It follows that D F (p||q) = F (p) -H q (p) (see Figure 2).

We now give some basic properties of Bregman divergences. The first property seems to be new. The others are well known. First, observe that, for most functions F , the associated Bregman divergence is not symmetric, i.e. D F (p||q) = D F (q||p) (the symbol || is put to emphasize this point, as is standard in information theory). The following lemma proves this claim.

Lemma 2 Let F be properly defined for D F to exist. Then D F is symmetric if and only if the Hessian ∇ 2 F is constant on X .

Proof: (⇒) From Eq. 1, the symmetry D F (p||q) = D F (q||p) yields:

F (p) = F (q) + 1 2 p -q, ∇F (q) + ∇F (p) . (2)
1 Lev M. Bregman historically pioneered this notion in the seminal work [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF] on minimization of a convex objective function under linear constraints. See http://www.math.bgu.ac.il/serv/segel/bregman.html. We gratefully acknowledge him for sending us this historical paper.

2 See Java applet at http://www.csl.sony.co.jp/person/nielsen/BregmanDivergence/ A Taylor expansion of F around q using the Lagrange form of the remainder also yields:

INRIA F X p q p q H q D F (p||q)
F (p) = F (q) + p -q, ∇F (q) + 1 2 (p -q) T ∇ 2 F (q)(p -q) + 1 6 p -q, ∇F 3 (r pq ) , (3)
with r pq on the line segment pq. Equations ( 2) and (3) yield the following constraint:

p -q, ∇F (p) = p -q, ∇F (q) + (p -q) T ∇ 2 F (q)(p -q) + 1 3 p -q, ∇F 3 (r pq ) . (4)
On the other hand, if we make the Taylor expansion of ∇F around q and then multiply both sides by pq, we separately obtain:

p -q, ∇F (p) = p -q, ∇F (q) + (p -q) T ∇ 2 F (q)(p -q) + 1 2 p -q, ∇F 3 (s pq ) ,
with s pq on the line segment pq. However, for this to equal Eq. ( 4), we must have pq, ∇F 3 (r pq ) = (3/2) p-q, ∇F 3 (s pq ) for each p and q in X . If we pick p and q very close to each other, this equality cannot be true, except when the third differentials are all zero on r pq and s pq . Repeating this argument over each subset of X having non zero measure, we obtain that the third differentials of F must be zero everywhere but on subsets of X with zero measure, which implies that the second differentials (the Hessian of F , ∇ 2 F ) are constant everywhere on X . (⇐) Assume the hessian ∇ 2 F is constant on X . In this case, because F is strictly convex, the Hessian ∇ 2 F is positive definite, and we can factor it as ∇ 2 F = P -1 DP where D is a diagonal matrix and P a unitary rotation matrix. Reasoning in the basis of X formed by P, each element x is mapped to Px, and we have

F (x) = i d i x 2 i
, where the d i 's are the diagonal coefficients of D. The symmetry of D F is then immediate (i.e., D F is a generalized quadratic distance).

Property 1 (Non-negativity) The strict convexity of generator function F implies that, for any p and q in X , D F (p||q) ≥ 0, with D F (p||q) = 0 if and only if p = q. Property 2 (Convexity) Function D F (p||q) is convex in its first argument p but not necessarily in its second argument q.

Bregman divergences can easily be constructed from simpler ones. For instance, multivariate Bregman divergences D F can be created from univariate generator functions coordinate-wise as

F (x) = d i=1 f i (x i ) with ∇F = [ df1 dx1 ... df d dx d ] T .
Because positive linear combinations of strictly convex and differentiable functions are strictly convex and differentiable functions, new generator functions (and corresponding Bregman divergences) can also be built as positive linear combinations of elementary generator functions. This is an important property as it allows to handle mixed data sets of heterogenous types in a unified framework.

Property 3 (Linearity) Bregman divergence is a linear operator, i.e., for any two strictly convex and differentiable functions F 1 and F 2 defined on X and for any λ ≥ 0:

D F1+λF2 (p||q) = D F1 (p||q) + λD F2 (p||q).
Property 4 (Invariance under linear transforms) G(x) = F (x) + a, x + b, with a ∈ R d and b ∈ R, is a strictly convex and differentiable function on X , and

D G (p||q) = D F (p||q).
Examples of Bregman divergences are the squared Euclidean distance (obtained for F (x) = x 2 and the generalized quadratic distance function F (x) = x T Qx where Q is a positive definite matrix. When Q is taken to be the inverse of the variance-covariance matrix, D F is the Mahalanobis distance, extensively used in computer vision. More importantly, the notion of Bregman divergence encapsulates various information measures based on entropic functions such as the Kullback-Leibler divergence based on the (unnormalized) Shannon entropy, or the Itakura-Saito divergence based on Burg entropy (commonly used in sound processing). Table 1 lists the main univariate Bregman divergences.

Legendre duality

We now turn to an essential notion of convex analysis: Legendre transform that will allow us to associate to any Bregman divergence a dual Bregman divergence.

Let F be a strictly convex and differentiable real-valued function on X . The Legendre transformation makes use of the duality relationship between points and lines to associate to F a convex conjugate function F * : R d → R given by [START_REF] Rockafellar | Convex Analysis[END_REF]:

F * (y) = sup x∈X { y, x -F (x)}.
The supremum is reached at the unique point where the gradient of G(x) = y, x -F (x) vanishes or, equivalently, when y = ∇F (x).

INRIA

Dom. X

Function F Gradient Inv. grad. Divergence DF (p||q)

R

Squared function

Squared loss (norm)

x 2 2x x 2 (p -q) 2 R+, α ∈ N Norm-like Norm-like α > 1 x α αx α-1 ( x α ) 1 α-1 p α + (α -1)q α -αpq α-1 R + Unnorm. Shannon entropy Kullback-Leibler div. (I-div.) x log x -x log x exp(x) p log p q -p + q Exponential Exponential loss R exp x exp x log x exp(p) -(p -q + 1) exp(q) R + * Burg entropy Itakura-Saito divergence -log x -1 x -1 x p q -log p q -1 [0, 1] Bit entropy Logistic loss x log x + (1 -x) log(1 -x) log x 1-x exp x 1+exp x p log p q + (1 -p) log 1-p 1-q Dual bit entropy Dual logistic loss R log(1 + exp x) exp x 1+exp x log x 1-x log 1+exp p 1+exp q -(p -q) exp q 1+exp q [-1, 1] Hellinger-like Hellinger-like - √ 1 -x 2 x √ 1-x 2 x √ 1+x 2 1-pq √ 1-q 2 - p 1 -p 2 Table 1: Some common univariate Bregman divergences D F .
As is well-known, F * is strictly convex. To see this, consider the epigraph epi(F * ), i.e. the set of points (y, z) such that F * (y) ≤ z. Clearly, (y, z) ∈ epi(F * ) iff G x (y) = y, x -F (x) ≤ z for all x ∈ X . Therefore, epi(F * ) = ∩ x∈X epi(G x ). Since G x (y) is an affine function, epi(G x ) is a half-space and epi(F * ) being the intersection of half-spaces is a convex set, which proves that F * is convex. The strict convexity follows from the fact that otherwise, F would not be differentiable in at least one point z ∈ X : at this point, y α , z -F (z) ≥ y α , x -F (x), ∀x ∈ X , and y α = αy 1 + (1 -α)y 2 , ∀α ∈ [0, 1], y 1 y 2 being a segment on which F * is not strictly convex. Thus, y 1 y 2 would be a subdifferential of F in z contradicting the fact that F is differentiable.

For convenience, we write x = ∇F (x) (omitting the F in the x notation as it should be clear from the context). Figure 3 gives a geometric interpretation of the Legendre transformation. Using this notation, Eq. 1 can be rewritten as

D F (p||q) = F (p) -F (q) -q , p -q . ( 5
)
Since F is a strictly convex and differentiable real-valued function on X , its gradient ∇F is well defined as well as its inverse ∇ -1 F . Writing X for the gradient space

{∇F (x) = x |x ∈ X }, the convex conjugate F * of F is the function: X ⊂ R d → R defined by F * (x ) = x, x -F (x). (6)
F : z = F (x) (0, -F * (y )) Deriving this expression, we get ∇F * (x ), dx = x, dx + x , dx -∇F (x), dx = x, dx = ∇ -1 F (x ), dx , from which we deduce that ∇F * = ∇ -1 F . From Eq. 6, we also deduce (F * ) * = F . From the above discussion, it follows that D F * is a Bregman divergence, which we call the Legendre dual divergence of D F . We have :

y z =< y , x > -F * (y )
Lemma 3 D F (p||q) = F (p) + F * (q ) -p, q = D F * (q ||p )
Proof: By Eq. 5, D F (p||q) = F (p) -F (q) -pq, q , and, according to Eq. 6, we have F (p) = p , p -F * (p ) and F (q) = q , q -F * (q ). Hence, D F (p||q) = p , p -F * (p )p, q + F * (q ) = D F * (q ||p ) since p = ∇F -1 ∇F (p) = ∇F * (p ).

Observe that, when D F is symmetric, D F * is also symmetric. The Legendre transform of the quadratic form

F (x) = 1 2 x T Qx, where Q is a symmetric invertible matrix, is F * (y) = 1 2 y T Q -1 y (corresponding divergences D F and D F * are both generalized quadratic distances).
To compute F * , we use the fact that ∇F * = ∇ -1 F and obtain F * as F * = ∇ -1 F . For example, the Hellinger-like measure is obtained by setting F (x) = -√ 1 -x 2 (see Table 1). The inverse gradient is x √ 1+x 2 and the dual convex conjugate is

xdx √ 1+x 2 = √ 1 + x 2 .
Integrating functions symbolically may be difficult or even not possible, and, in some cases, it will be required to approximate numerically the inverse gradient ∇ -1 F (x).

INRIA

Let us consider the univariate generator functions defining the divergences of Table 1. Both the squared function F (x) = x 2 and Burg entropy F (x) = -log x are self-dual, i.e. F = F * . This is easily seen by noticing that the gradient and inverse gradient are identical (up to some constant factor).

For the exponential function F (x) = exp x, we have F * (y) = y log y-y (the unnormalized Shannon entropy) and for the dual bit entropy F (x) = log(1 + exp x), we have F * (y) = y log y 1-y + log(1 -y), the bit entropy. Note that the bit entropy function is a particular Bregman generator satisfying F (x) = F (1 -x).

Symmetrized Bregman divergences

For non-symmetric d-variate Bregman divergences D F , we define the symmetrized divergence

S F (p, q) = S F (q, p) = 1 2 (D F (p||q) + D F (q||p)) = 1 2 p -q, p -q .
An example of such a symmetrized divergence is the symmetric Kullback-Leibler divergence (SKL) widely used in computer vision and sound processing (see for example [START_REF] Le | Ent-Boost: Boosting Using Entropy Measure for Robust Object Detection[END_REF]).

A key observation is to note that the divergence S F between two points of X can be measured as a divergence in X ×X ⊂ R 2d . More precisely, let x = [x x ] T be the 2d-dimensional vector obtained by stacking the coordinates of x on top of those of x , the gradient of F at x. We have :

Theorem 1 S F (p, q) = 1 2 D F (p||q) where F (x) = F (x) + F * (x )
and D F is the Bregman divergence defined over X × X ⊂ R 2d for the generator function F .

Proof: Using Lemma 3, we have

S F (p, q) = 1 2 (D F (p||q) + D F (q||p)) = 1 2 (D F (p||q) + D F * (p ||q )) = 1 2 D F (p||q)
It should be noted that x lies on the d-manifold

X = {[x x ] T | x ∈ R d } of R 2d
. Note also that S F (p, q) is symmetric but not a Bregman divergence in general while D F is a non symmetric Bregman divergence in X × X .

Exponential families

Parametric statistical spaces and exponential families

A statistical space X is an abstract space where coordinates of vector points θ ∈ X encode the parameters of statistical distributions. The dimension d = dim X of the statistical space coincides with the finite number of free parameters of the distribution laws. For example, the space

X = {[µ σ] T | (µ, σ) ∈ R × R + * } of univariate normal distributions N (µ, σ
) is a 2D parametric statistical space, extensively studied in information geometry [START_REF] Amari | Methods of Information Geometry[END_REF] under the auspices of differential geometry. A prominent class of distribution families called the exponential families E F [START_REF] Amari | Methods of Information Geometry[END_REF] admits the same canonical probability distribution function

p(x|θ) def = exp{ θ, f (x) -F (θ) + C(x)}, (7)
where f (x) denotes the sufficient statistics and θ ∈ X represents the natural parameters. Space X is thus called the natural parameter space and, since log x p(x|θ)dx = log 1 = 0, we have 2 summarizes the various relevant parts of the canonical decompositions of some of these usual statistical distributions. Observe that the product of any two distributions of the same exponential family is another exponential family distribution that may not have anymore a nice parametric form (except for products of normal distribution pdfs that yield again normal distribution pdfs). Thus exponential families provide a unified treatment framework of common distributions. Note, however, that the uniform distribution does not belong to the exponential families.

F (θ) = log x exp{ θ, f (x) + C(x)

Kullback-Leibler divergence of exponential families

In such statistical spaces X , a basic primitive is to measure the distortion between any two distributions. The Kullback-Leibler divergence (also called relative entropy or information divergence, I-divergence) is a standard information-theoretic measure between two statistical distributions d 1 and d 2 defined as KL(d

1 ||d 2 ) def = x d 1 (x) log d1(x)
d2(x) dx. This statistical measure is not symmetric nor does the triangle inequality holds.

The link with Bregman divergences comes from the remarkable property that the Kullback-Leibler divergence of any two distributions of the same exponential family with respective natural parameters θ p and θ q is obtained from the Bregman divergence induced by the cumulant function of that family by swapping arguments. By a slight abuse of notations, we denote by KL(θ p ||θ q ) the oriented Kullback-Leibler divergence between the probability density functions defined by the respective natural parameters, i.e. KL(θ

p ||θ q ) def = x p(x|θ p ) log p(x|θp) p(x|θq) dx.
The following theorem is the extension to the continuous case of a result mentioned in [6].

INRIA

Exponential family

Canonical probability density function

: exp{ θ, f (x) -F (θ) + C(x)} Natural Sufficient Cumulant function F (θ) Dens. Norm. parameters θ statistics f (x) C(x)
Bernouilli B(q) (Tossing coin with Pr(heads) = q and Pr(tails

) = 1 -q) log q 1-q x log(1 + exp θ) 0
Multinomial M(q 1 , ..., q d+1 ) (Extend Bernouilli with Pr(x i ) = q i and i q i = 1)

θ i = log qi 1- P d j=1 qi f i (x) = x i log(1 + d i=1 exp θ i ) 0 Beta β(θ 1 , θ 2 ) (Bernouilli conjugate prior) [θ 1 θ 2 ] T [log x log(1 -x)] T log B(θ 1 + 1, θ 2 + 1) 0
F (θ) = log Γ(θ1+1)Γ(θ2+1) Γ(θ1+θ2-2) (with Γ(x) = ∞ 0 t x-1 exp(-t)dt) = (x -1)Γ(x -1)) Univariate Normal N (µ, σ 2 ) [ µ σ 2 -1 2σ 2 ] T [x x 2 ] T - θ 2 1 4θ2 + 1 2 log(-π θ2 ) 0 Multivariate Normal N (µ, Σ) [Σ -1 µ -1 2 Σ -1 ] [x xx T ] 1 2 µ T Σ -1 µ + 1 2 log det(2πΣ) 0 Rayleigh R(σ 2 ) (used in ultrasound imageries) -1 2σ 2 x 2 log -1 2θ log x Laplacian L(θ) (used in radioactivity decay) θ -x -log θ 0 Poisson P(λ) (counting process) log λ x exp θ -log x! Gamma γ(θ 1 , θ 2 ) (waiting times in Poisson processes) [θ 1 θ 2 ] T [log x x] T log Γ(θ 1 + 1) + (θ 2 + 1) log(-θ 2 ) 0 Dirichlet D(α) (varying proportion model ||x|| = 1, conjugate prior of Multinomial) θ i = α i -1 f i (x) = log x i log Γ( i θ i + d) -i Γ(θ i + 1) 0
Table 2: Canonical decompositions of usual exponential families.

Theorem 2 The Kullback-Leibler divergence of any two distributions of the same exponential family with natural parameters θ p and θ q is obtained from the Bregman divergence induced by the cumulant function F as: KL(θ p ||θ q ) = D F (θ q ||θ p ).

Before proving the theorem, we note that

∇F (θ) = x f (x) exp{ θ, f (x) -F (θ) + C(x)}dx . (8)
The coordinates of

µ def = ∇F (θ) = [ x f (x)p(x|θ)dx] = E θ (f (x)
) are called the expectation parameters. As an example, consider the univariate normal distribution N (µ, σ) with sufficient statistics [x x 2 ] T (see Table 2). The expectation parameters are

µ = ∇ F (θ) = [µ µ 2 + σ 2 ] T , where µ = x x p(x|θ)dx and µ 2 + σ 2 = x x 2 p(x|θ)dx.
We now prove the theorem. Proof:

KL(θ p ||θ q ) = x p(x|θ p ) log p(x|θ p ) p(x|θ q ) dx = x p(x|θ p )(F (θ q ) -F (θ p ) + θ p -θ q , f (x) )dx = x p(x|θ p ) (D F (θ q ||θ p ) + θ q -θ p , ∇F (θ p ) + θ p -θ q , f (x) ) dx = D F (θ q ||θ p ) + x p(x|θ p ) θ q -θ p , ∇F (θ p ) -f (x) )dx = D F (θ q ||θ p ) - x p(x|θ p ) θ q -θ p , f (x) dx + θ q -θ p , ∇F (θ p ) (Eq. 8) = D F (θ q ||θ p )

Dual parameterizations and dual divergences

The notion of dual Bregman divergences introduced earlier and dual parameterizations extend naturally to statistical spaces. Since, µ = ∇F (θ) (Eq. 8), the convex conjugate of F (θ) is F * (µ) = θ, µ -F (θ) (Eq. 6). From Lemma 3, we then deduce the following theorem.

Theorem 3 D F (θ p ||θ q ) = D F * (µ q ||µ p ) where F * denote the convex conjugate of F .
Table 3 presents some examples of dual parameterizations of exponential families (i.e., the natural θ-parameters and expectation µ-parameters and dual Legendre cumulant functions), and describe the corresponding Bregman divergences induced by the Kullback-Leibler divergences.

Finally, we would like to point out that Banerjee et al. [6] have shown that there is a bijection between the regular exponential families and a subset of the Bregman divergences called regular Bregman divergences.

INRIA

Bernouilli dual divergences: Logistic loss/binary relative entropy

F (θ) = log(1 + exp θ) D F (θ||θ ) = log 1+exp θ 1+exp θ -(θ -θ ) exp θ 1+exp θ f (θ) = exp θ 1+exp θ = µ F * (µ) = µ log µ + (1 -µ) log(1 -µ) D F * (µ ||µ) = µ log µ µ + (1 -µ) log 1-µ 1-µ f * (µ) = log µ 1-µ = θ
Poisson dual divergences: Exponential loss/Unnormalized Shannon entropy

F (θ) = exp θ D F (θ||θ ) = exp θ -exp θ -(θ -θ ) exp θ f (θ) = exp θ = µ F * (µ) = µ log µ -µ D F * (µ ||µ) = µ log µ µ + µ -µ f * (µ) = log µ = θ
Table 3: Examples of dual parameterizations of exponential families and their corresponding Kullback-Leibler (Bregman) divergences for the Bernoulli and Poisson distributions.

Elements of Bregman geometry

In this section, we discuss several basic geometric properties that will be useful when studying Bregman Voronoi diagrams. Specifically, we characterize Bregman bisectors, Bregman balls and Bregman geodesics. Since Bregman divergences are not symmetric, we describe several types of Bregman bisectors in §3.1. We subsequently characterize Bregman balls by using a lifting transform that extends a construction well-known in the Euclidean case ( §3.2). Finally, we characterize geodesics and show an orthogonality property between bisectors and geodesics in §3.3.

Bregman bisectors

Since Bregman divergences are not symmetric, we can define several types of bisectors. The Bregman bisector of the first type is defined as

H F (p, q) = {x ∈ X | D F (x||p) = D F (x||q)}.
Similarly, we define the Bregman bisector of the second type as

H F (p, q) = {x ∈ X | D F (p||x) = D F (q||x)}.
These bisectors are identical when the divergence is symmetric. However, in general, they are distinct, the bisectors of the first type being linear while the bisectors of the second type are potentially curved (but always linear in the gradient space, hence the notation). More precisely, we have the following lemma Lemma 4 The Bregman bisector of the first type H F (p, q) is the hyperplane of equation:

H F (p, q) : x, p -q + F (p) -p, p -F (q) + q, q = 0
The Bregman bisector of the second type H F (p, q) is the hypersurface of equation

H F (p, q) : x , q -p + F (p) -F (q) = 0
(a hyperplane in the gradient space X ).

It should be noted that p and q lie necessarily on different sides of

H F (p, q) since H F (p, q)(p) = -D F (p||q) < 0 and H F (p, q)(q) = D F (q||p) > 0.
From Lemma 3, we know that D F (x||y) = D F * (y ||x ) where F * is the convex conjugate of F . We therefore have The bisector H F (p, q) for the symmetrized Bregman divergence S F is given by H F (p, q) : x, qp + x , qp + p, pq, q = 0. Such a bisector is not linear in x nor in x . However, we can observe that the expression is linear in x = [x x ] T . Indeed, proceeding as we did in §2.3, we can rewrite the above equation as

H F (p, q) = ∇ -1 F (H F * (q , p )), H F (p, q) = ∇ -1 F (H F * (q , p )).
H F (p, q) : x x , q -p q -p + p, p -q, q = 0.
which shows that H F (p, q) is the projection on X of the intersection of the hyperplane

H(p, q) of R 2d with the d-dimensional manifold X = {x = [x x ] T | x ∈ X }.

Bregman spheres and the lifting map

We define the Bregman balls of, respectively, the first and the second types as

B F (c, r) = {x ∈ X | D F (x||c) ≤ r} and B F (c, r) = {x ∈ X | D F (c||x) ≤ r}
The Bregman balls of the first type are convex while this is not necessarily true for the balls of the second type as shown in Fig. 5 for the Itakura-Saito divergence (defined in Table 1). The associated bounding Bregman spheres are obtained by replacing the inequalities by equalities.

From Lemma 3, we deduce that

B F (c, r) = ∇ -1 F (B F * (c , r)). (9)
Let us now examine a few properties of Bregman spheres using a lifting transformation that generalizes a similar construct for Euclidean spheres (see [START_REF] Boissonnat | Algorithmic Geometry[END_REF][START_REF] Teillaud | The space of spheres, a geometric tool to unify duality results on voronoi diagrams[END_REF]).

Let us embed the domain X in X = X × R ⊂ R d+1 using an extra dimension denoted by the Z-axis. For a point x ∈ X , recall that x = (x, F (x)) denotes the point obtained by lifting x onto F (see Figure 1). In addition, write Proj X (x, z) = x for the projection of a point of X onto X .

INRIA

Source space X

Gradient space X (a) Gradient Space: Itakura-Saito dual p'(-1.88760873,-1.38808518) q'(-1.16516903,-3.43833618) D*(p',q')=0.44835617 D*(q',p')=0.66969016 Let p ∈ X and H p be the hyperplane tangent to F at point p of equation

z = H p (x) = x -p, p + F (p),
and let H ↑ p denote the halfspace above H p consisting of the points x = [x z] T ∈ X such that z > H p (x). Let σ(c, r) denote either the first-type or second-type Bregman sphere centered at c with radius r (i.e., ∂B F (c, r) or ∂B F (c, r)).

The lifted image σ of a Bregman sphere σ is σ = {(x, F (x)), x ∈ σ}. We associate to a Bregman sphere σ = σ(c, r) of X the hyperplane

H σ : z = x -c, c + F (c) + r, (10)
parallel to H c and at vertical distance r from H c (see Figure 6). Observe that H σ coincides with H c when r = 0, i.e. when sphere σ is reduced to a single point.

Lemma 5 σ is the intersection of F with H σ . Conversely, the intersection of any hyperplane H with F projects onto X as a Bregman sphere. More precisely, if the equation of H is z = x, a + b, the sphere is centered at c = ∇ -1 F (a) and its radius is a, c -F (c) + b.

Proof: The first part of the lemma is a direct consequence of the fact that D F (x||y) is measured by the vertical distance from x to H y (see Lemma 1). For the second part, we consider the hyperplane H parallel to H and tangent to F. From Eq. 10, we deduce

a = c . The equation of H is thus z = x -∇ -1 F (a), a + F (∇ -1 F (a)
). It follows that the divergence from any point of σ to c, which is equal to the vertical distance between H and H , is

∇ -1 F (a), a -F (∇ -1 F (a)) + b = a, c -F (c) + b. INRIA (a) Squared Euclidean distance (b) Itakura-Saito divergence
F G = F ∩ G and H G = H ∩ G .
Note that F G is the graph of the restriction of F to G and that H G is a hyperplane of G . We can therefore apply Lemma 5 in G , which proves the lemma for Bregman spheres of the first type. The case of Bregman spheres of the second type follows from the duality of Eq. 9.

Union and intersection of Bregman balls

Theorem 4 The union of n Bregman balls has combinatorial complexity Θ(n

) and can be computed in time Θ(n log n + n d+1 2

).

Proof: To each ball, we can associate its bounding Bregman sphere σ i which, by Lemma 5, is the projection by Proj X of the intersection of F with a hyperplane H σi . The points of F that are below H σi projects onto points that are inside the Bregman ball bounded by σ i . Hence, the union of the balls is the projection by Proj X of the complement of F ∩ H ↑ where

H ↑ = ∩ n i=1 H ↑ σi .
H ↑ is a convex polytope defined as the intersection of n half-spaces. The theorem follows from McMullen's theorem that bounds the number of faces of a polytope [START_REF] Mcmullen | The maximum numbers of faces of a convex polytope[END_REF], and Chazelle's optimal convex hull/half-space intersection algorithm [START_REF] Chazelle | An optimal convex hull algorithm in any fixed dimension[END_REF]. The result for the balls of the second type is deduced from the result for the balls of the first type and the duality of Eq. 9.

Very similar arguments prove the following theorem (just replace H ↑ σi by the complementary halfspace H ↓ σi ).

Theorem 5 The intersection of n Bregman balls has combinatorial complexity Θ(n

) and can be computed in time Θ(n log n + n d+1 2

).

Circumscribing Bregman spheres. There exists, in general, a unique Bregman sphere passing through d + 1 points of R d . This is easily shown using the lifting map since, in general, there exists a unique hyperplanes of R d+1 passing through d + 1 points. The claim then follows from Lemma 5. Deciding whether a point x falls inside, on or outside a Bregman sphere σ ∈ R d passing through d + 1 points of p 0 , ..., p d will be crucial for computing Bregman Voronoi diagrams and associated triangulations. The lifting map immediately implies that such a decision task INRIA reduces to determining the orientation of the simplex (p 0 , ..., pd , x) of R d+1 , which in turn reduces to evaluating the sign of the determinant of the (d + 2) × (d + 2) matrix (see [START_REF] Nielsen | Visual Computing: Geometry, Graphics, and Vision[END_REF])

InSphere(x; p 0 , ..., p d ) = 1 ... 1 1 p 0 ... p d x F (p 0 ) ... F (p d ) F (x)
If one assumes that the determinant 1 ... 1 p 0 ... p d is non-zero, InSphere(x; p 0 , ..., p d ) is negative, null or positive depending on whether x lies inside, on, or outside σ.

Projection, orthogonality and geodesics

We start with an easy property of Bregman divergences.

Property 5 (Three-point property) For any triple p, q and r of points of X , we have:

D F (p||q) + D F (q||r) = D F (p||r) + p -q, r -q .
The following lemma characterizes the Bregman projection of a point onto a closed convex set W. We now introduce the notion of Bregman orthogonality. We say that pq is Bregman orthogonal to qr iff D F (p||q) + D F (q||r) = D F (p||r) or equivalently (by the Three-point property), if and only if pq, rq = 0. Observe the analogy with Pythagoras' theorem in Euclidean space (see Figure 7). Note also that the orthogonality relation is not symmetric: the fact that pq is Bregman orthogonal to qr does not necessarily imply that qr is Bregman orthogonal to pq. More generally, we say that I ⊆ X is Bregman orthogonal to J ⊆ X (I ∩ J = ∅) iff for any p ∈ I and r ∈ J, there exists a q ∈ I ∩ J such that pq is Bregman orthogonal to qr.

Notice that orthogonality is preserved in the gradient space. Indeed, since pq, rq = rq , pq , pq is Bregman orthogonal to qr iff r q is Bregman orthogonal to q p . Let Γ F (p, q) be the image by ∇ -1 F of the line segment p q , i.e. By analogy, we rename the line segment pq as

Γ F (p, q) = {x ∈ X : x = (1 -λ)p + λq , λ ∈ [0, 1]}.
Λ(p, q) = {x ∈ X : x = (1 -λ)p + λq, λ ∈ [0, 1]} In the Euclidean case (F (x) = 1 2 x 2 ), Γ F (p, q) = Λ(p, q)
is the unique geodesic path joining p to q and it is orthogonal to the bisector H F (p, q). For general Bregman divergences, we have similar properties as shown next.

Lemma 8 Γ F (p, q) is Bregman orthogonal to the Bregman bisector H F (p, q) while Λ(p, q) is Bregman orthogonal to H F * (p, q). Proof: Since p and q lie on different sides of H F (p, q), Γ F (p, q) must intersect H F (p, q). Fix any distinct x ∈ Γ(p, q) and y ∈ H F (p, q), and let t ∈ Γ(p, q) ∩ H F (p, q). To prove the first part of the lemma, we need to show that yt, xt = 0.

Since t and x both belong to ∈ Γ F (p, q), we have tx = λ(p -q ), for some λ ∈ R, and, since y and t belong to H F (p, q), we deduce from the equation of H F (p, q) that yt, pq = 0. We conclude that yt, xt = 0, which proves that Γ F (p, q) is indeed Bregman orthogonal to H F (p, q).

The second part of the lemma is easily proved by using the fact that orthogonality is preserved in the gradient space as noted above.

Figure 8 shows Bregman bisectors and their relationships with respect to Λ(p, q) and Γ F (p, q).

We now focus on characterizing Bregman geodesics. First, recall that a parameterized curve C between two points p 0 and p 1 is defined as a set C = {p λ } 1 λ=0 , which is continuous. In Riemannian geometry, geodesics are the curves that minimize the arc length with respect to the Riemannian metric [START_REF] Amari | Methods of Information Geometry[END_REF][START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF]. Since embedding X with a Bregman divergence does not yield a metric space, we define the following curve lengths:

Γ (C) = 1 λ=0 D F (p 0 ||p λ )dλ , (11)
INRIA

Λ (C) = 1 λ=0 D F (p λ ||p 0 )dλ. ( 12
)
We now characterize the dual pair of geodesics and their lengths as follows:

Lemma 9 Curve Γ F (p 0 , p 1 ) (respectively straight line segment Λ(p 0 , p 1 )) minimizes

H - p : y, h < b. Now, H p in- tersects Γ(p 0 , p 1 ) since H p separates p 0 from p 1 . Indeed, H p (p 0 ) = p 0 -p, p 0 -p = D F (p 0 ||p) + D F (p||p 0 ) > 0 and H p (p 1 ) = p 1 -p, p 0 -p = λ-1 λ p 1 -p, p 1 -p < 0 where p = λp 0 + (1 -λ)p 1 (with λ ∈]0, 1[).
Therefore any connected path C joining p 0 to p 1 has to intersect H p .

To finish up, consider function f

: [0, 1] → C with f (0) = p 0 , f (1) = p 1 , and f (λ) ∈ C ∩ H p λ otherwise, where it is understood that p λ is hereafter a point of Γ F (p 0 , p 1 ). Since f (λ) ∈ H p(λ) , we have D F (f (λ)||p 0 ) = D F (f (λ)||p λ ) + D F (p λ ||p 0 ) ≥ D F (p λ ||p 0 ), with equality if and only if f (λ) = p λ . Thus we have Γ (Γ F (p 0 , p 1 )) = 1 λ=0 D F (p λ ||p 0 )dλ ≤ 1 λ=0 D F (f (λ)||p 0 )dλ ≤ Γ (C) .
The case of Λ(p 0 , p 1 ) follows similarly from Legendre convex duality.

Corollary 1 Since Γ F (p 0 , p 1 ) = Γ F (p 1 , p 0 ) (respectively, since Λ(p 0 , p 1 ) = Λ(p 1 , p 0 )) we deduce that Γ F (p 0 , p 1 ) minimizes also Observe also that Γ F (p, q) is the unique geodesic path joining p to q in X for the metric image by ∇ -1 F of the Euclidean metric.

Finally, we give a characterization of these geodesics in information-theoretic spaces. Recall that Banerjee et al. [6] showed that Bregman divergences are in bijection with exponential families. This was emphasized by Theorem 2 that proved that the Kullback-Leibler divergence of probability density functions of the same exponential family E F is a Bregman divergence D F for the cumulant function F . From this standpoint, Λ(p, q) and Γ F (p, q) minimize the total Kullback-Leibler divergence, a characteristic that we choose to call the information length of a curve. Since the Kullback-Leibler divergence is not symmetric, this justifies for the existence of two geodesics, one which appears to be linear when parameterized with the natural affine coordinate system (θ), and the other that is linear in the expectation affine coordinate system (µ). See also [START_REF] Amari | Methods of Information Geometry[END_REF].

INRIA

Bregman Voronoi diagrams

Let S = {p 1 , ..., p n } be a finite point set of X ⊂ R d . To each point p i is attached a d-variate continuous function D i defined over X . We define the lower envelope of the functions as the graph of min 1≤i≤n D i and their minimization diagram as the subdivision of X into cells such that, in each cell, arg min i f i is fixed.

The Euclidean Voronoi diagram is the minimization diagram for D i (x) = xp i 2 . In this section, we introduce Bregman Voronoi diagrams as minimization diagrams of Bregman divergences (see Figure 10).

We define three types of Bregman Voronoi diagrams in §4.1. We establish a correspondence between Bregman Voronoi diagrams and polytopes in §4.2 and with power diagrams in §4.3. These correspondences lead to tight combinatorial bounds and efficient algorithms. Finally, in §4.4, we give two generalizations of Bregman Voronoi diagrams; k-order and k-bag diagrams.

We note S = {∇ F (p i ), i = 1, . . . , n} the gradient point set associated to S.

Three types of diagrams

F (p i ) def = {x ∈ X | D F (x||p i ) ≤ D F (x||p j ) ∀p j ∈ S.}
Since the Bregman bisectors of the first-type are hyperplanes, the cells of any diagram of the first-type are convex polyhedra. Therefore, first-type Bregman Voronoi diagrams are affine diagrams [4,5].

Similarly, the minimization diagram of the D i , i = 1, . . . , n, is called the second-type Bregman Voronoi diagram of S, which we denote by vor F (S). A cell in vor F (S) is associated to each site p i and is defined as above with permuted divergence arguments:

vor F (p i ) def = {x ∈ X | D F (p i ||x) ≤ D F (p j ||x) ∀p j ∈ S.}
In contrast with the diagrams of the first-type, the diagrams of the second type have, in general, curved faces. For asymmetric Bregman divergences D F , we can further consider the symmetrized Bregman divergence S F = D F and define a third-type Bregman Voronoi diagram vor F (S). The cell of vor F (S) associated to site p i is defined as:

vor F (p i ) def = {x ∈ X | S F (x, p i ) ≤ S F (x, p j ) ∀p j ∈ S.}
From the Legendre duality between divergences, we deduce correspondences between the diagrams of the first and the second types. As usual, F * is the convex conjugate of F .

(p i ) = {x ∈ X | D F * (p i ||x ) ≤ D F * (p j ||x ) ∀p j ∈ S } = ∇ -1 F (vor F * (p i ))
. The proof of the second part follows the same path.

Hence, constructing the second-type curved diagram vor F (S) reduces to constructing an affine diagram in the gradient space X (and map the cells by ∇F -1 ).

INRIA

Let us end this section by considering the case of symmetrized Bregman divergences introduced in §2.

Bregman Voronoi diagrams from polytopes

Let H pi , i = 1, . . . , n, denote the hyperplanes of X defined in §3.2. For any x ∈ X , we have following Lemma 1

D F (x||p i ) ≤ D F (x||p j ) ⇐⇒ H pi (x) ≥ H pj (x).
The first-type Bregman Voronoi diagram of S is therefore the maximization diagram of the n linear functions H pi (x) whose graphs are the hyperplanes H pi (see Figure 10). Equivalently, we have Theorem 6 The first-type Bregman Voronoi diagram vor F (S) is obtained by projecting by Proj X the faces of the (d + 1)-dimensional convex polyhedron H = ∩ i H ↑ pi of X + onto X .

From McMullen's upperbound theorem [START_REF] Mcmullen | The maximum numbers of faces of a convex polytope[END_REF] and Chazelle's optimal half-space intersection algorithm [START_REF] Chazelle | An optimal convex hull algorithm in any fixed dimension[END_REF], we know that the intersection of n halfspaces of R d has complexity Θ(n ) and can be computed in optimal time Θ(n log n + n d+1 2

). The third-type Bregman Voronoi diagram for the symmetrized Bregman divergence of a set of n d-dimensional points has complexity O(n d ) and can be obtained in time O(n d ).

Apart from Chazelle's algorithm, several other algorithms are known for constructing the intersection of a finite number of halfplanes, especially in the 2-and 3-dimensional cases. See [START_REF] Boissonnat | Algorithmic Geometry[END_REF]5] for further references.

D i (x) = ||p i -x|| 2 -r 2 .
The power bisector of any two balls B(p i , r i ) and B(p j , r j ) is the radical hyperplane of equation 2 x, p j -p i +||p i || 2 -||q j || 2 +r 2 j -r 2 i = 0. Thus power diagrams are affine diagrams. In fact, as shown by Aurenhammer [3,[START_REF] Boissonnat | Algorithmic Geometry[END_REF], any affine diagram is identical to the power diagram of a set of corresponding balls. In general, some balls may have an empty cell in their power diagram.

Since Bregman Voronoi diagrams of the first type are affine diagrams, Bregman Voronoi diagrams are power diagrams [3,[START_REF] Boissonnat | Algorithmic Geometry[END_REF] in disguise. The following theorem makes precise the correspondence between Bregman Voronoi diagrams and power diagrams (see Figure 11).

Theorem 8

The first-type Bregman Voronoi diagram of n sites is identical to the power diagram of the n Euclidean spheres of equations

x -p i , x -p i = p i , p i + 2(F (p i ) -p i , p i ), i = 1, . . . , n. Proof: We have D F (x||p i ) ≤ D F (x||p j ) ⇐⇒ -F (p i ) -x -p i , p i ≤ -F (p j ) -x -p j , p j
Multiplying twice the last inequality, and adding x, x to both sides yields

x, x -2 x, p i -2F (p i ) + 2 p i , p i ≤ x, x -2 x, p j -2F (p j ) + 2 p j , p j ⇐⇒ x -p i , x -p i -r 2 i ≤ x -p j , x -p j -r 2 j ,
where r 2 i = p i , p i + 2(F (p i ) -p i , p i ) and r 2 j = p j , p j + 2(F (p j ) -p j , p j ). The last inequality means that the power of x with respect to the Euclidean (possibly imaginary) ball B(p i , r i ) is no more than the power of x with respect to the Euclidean (possibly imaginary) ball B(p j , r j ).

As already noted, for F (x) = 1 2 x 2 , vor F (S) is the Euclidean Voronoi diagram of S. Accordingly, the theorem says that the centers of the spheres are the p i and r 2 i = 0 since p i = p i . Figure 11 displays affine Bregman Voronoi diagrams3 and their equivalent power diagrams for the squared Euclidean, Kullback-Leibler and exponential divergences.

Note that although the affine Bregman Voronoi diagram obtained by scaling the divergence D F by a factor λ > 0 does not change, the equivalent power diagrams are not strictus senso identical since the centers of corresponding Euclidean balls and radii are mapped differently. See the example of the squared Euclidean distance depicted in Figure 11(a). Since Power diagrams are well defined "everywhere", this equivalence relationship provides a natural way to extend the scope of definition of Bregman Voronoi diagrams from X ⊂ R d to the full space R d . (That is, Bregman Voronoi diagrams are power diagrams restricted to X .)

To check that associated balls may be potentially imaginary, consider for example, the Kullback-Leibler divergence. The Bregman generator function is F (x) = i x i log x i and the gradient is

∇F (x) = [log x 1 ... log x d ] T . A point p = [p 1 ... p d ] T ∈ X maps to a Euclidean ball of center p = [log p 1 ... log p d ] T with radius r 2 p = i (log 2 p i -2p i ).
Thus for points p with coordinates p i > 1 2 log p 2 i for i ∈ {1, ..., d}, the squared radius r 2 p is negative, yielding an imaginary ball. See Figure 11(b).

It is also to be observed that not all power diagrams are Bregman Voronoi diagrams. Indeed, in power diagrams, some balls may have empty cells while each site has necessarily a non empty cell in a Bregman Voronoi diagram (See Figure 11 and Section 4.4 for a further discussion at this point).

Since there exist fast algorithms for constructing power diagrams [START_REF] Pion | 3d triangulation data structure[END_REF], Theorem 8 provides an efficient way to construct Bregman Voronoi diagrams.

Generalized Bregman divergences and their Voronoi diagrams

Weighted Bregman Voronoi diagrams

Let us associate to each site p i a weight w i ∈ R. We define the weighted divergence between two weighted points as WD F (p i ||p j ) def = D F (p i ||p j ) + w i -w j . We can define bisectors and weighted Bregman Voronoi diagrams in very much the same way as for non weighted divergences. The Bregman Voronoi region associated to the weighted point (p i , w i ) is defined as vor

F (p i , w i ) = {x ∈ X | D F (x||p i ) + w i ≤ D F (x||p j ) + w j ∀p j ∈ S}.
Observe that the bisectors of the first-type diagrams are still hyperplanes and that the diagram can be obtained as the projection of a convex polyhedron or as the power diagram of a finite set of balls. The only difference with respect to the construction of Section 4.2 is the fact that now the hyperplanes H pi are no longer tangent to F since they are shifted by a z-displacement of length w i . Hence Theorem 7 extends to weighted Bregman Voronoi diagrams.

Theorem 9

The weighted Bregman Voronoi diagrams of type 1 or 2 of a set of n d-dimensional points have complexity Θ(n d+1 2

) and can be computed in optimal time Θ(n log n + n d+1 2

).

k-order Bregman Voronoi diagrams

We define the k-order Bregman Voronoi diagram of n punctual sites of X as the subdivision of X into cells such that each cell is associated to a subset T ⊂ S of k sites and consists of the points of X whose divergence to any site in T is less than the divergence to the sites not in T . Similarly to the case of higher-order Euclidean Voronoi diagrams, we have: Proof: Let S 1 , S 2 , . . . denote the subsets of k points of S and write

D i (x) = 1 k pj ∈Si D F (x||p j ) = F (x) - 1 k pj ∈Si F (p j ) + 1 k pj ∈Si x -p j , p j = F (x) -F (c i ) -x -c i , c i + w i = W D F (x||c i )
where

c i = ∇ -1 F 1 k j∈Si p j and the weight associated to c i is w i = F (c i ) -c i , c i - 1 k j∈Si F (p j ) + p j , p j . Hence, S i is the set of the k nearest neighbors of x iff D i (x) ≤ D j (x)
i ] T inducing a Bregman divergence D Fα i (x||p i ) def = D αi (x||p i ) anchored at that site. Let us consider the first-type of k-bag Bregman Voronoi diagram (kbag BVD for short). The first-type bisector K F (p i , p j ) of two weighted points (p i , α i ) and (p j , α j ) is the locus of points x at equidivergence to p i and p j . That is,

K F (p i , p j ) = {x ∈ X | D αi (x||p i ) = D αj (x||p j )}.
The equation of the bisector is simply obtained using the definition of Bregman divergences (Eq. 1) as

F αi (x) -F αi (p i ) -x -p i , ∇F αi (p i ) = F αj (x) -F αj (p j ) -x -p j , ∇F αi (p j ) .

This yields the equation of the first-type bisector K

F (p i , p j ) k l=1 (α (l) i -α (l) j )F l (x) -x, ∇F αj (p j ) -∇F αi (p i ) + c = 0, (13)

INRIA

where c is a constant depending on weighted sites (p i , α i ) and (p j , α j ). Note that the equation of the first-type k-bag BVD bisector is linear if and only if α i = α j (i.e., the case of standard BVDs).

Let us consider the linearization lifting x → x = [x F 1 (x) ... F k (x)] T that maps a point x ∈ R d into a point x in R d+k . Then Eq. 13 becomes linear, namely x, a + c = 0 with

a = ∇F αj (p j ) -∇F αi (p i ) α i -α j ∈ R d+k .
That is, first-type bisectors of a k-bag BVD are hyperplanes of R d+k . Therefore the complexity of a k-bag Voronoi diagram is at most O(n k+d 2

), since it can be obtained as the intersection of the affine Voronoi diagram in R d+k with the convex d-dimensional ) and can be computed within the same time bound.

submanifold {x = [x F 1 (x) ... F k (x)] T | x ∈ R d }.
Further, using the Legendre transform, we define a second-type (dual) k-bag BVD. We have

∇F α = k l=1 α l ∇F l and F * α = ∇F -1 α . (Observe that F * α = k l=1 α l F * l in general.
) k-bag Bregman Voronoi diagrams are closely related to the anisotropic diagrams of Labelle and Shewchuk [START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF] that associate to each point x ∈ X a metric tensor M x which tells how lengths and angles should be measured from the local perspective of x. Labelle and Shewchuk relies on a deformation tensor (ideally defined everywhere) to compute the distance between any two points p and q from the perspective of x as d x (p, q) = (p -q) T M x (p -q). Let d x (p) = d x (x, p). The anisotropic Voronoi diagram, which approximates the ideal but computationally prohibitive Riemannian Voronoi diagram, is defined as the arrangement of the following anisotropic Voronoi cells:

Vor(p i ) = {x ∈ X | d pi (x) ≤ d pj (x) ∀j ∈ {1, ..., n}}, ∀i ∈ {1, ..., n}.
It follows that all anisotropic Voronoi cells are non-empty as it is the case for k-bag Bregman Voronoi diagrams.

Hence, the site weights of a k-bag Bregman Voronoi diagram sparsely define a tensor divergence that indicates how divergences should be measured locally from the respective bag of divergences. Noteworthy, our study of k-bag Bregman Voronoi diagrams shows that the anisotropic Voronoi diagram also admits a second-type anisotropic Voronoi diagram, induced by the respective dual Legendre functions of the Bregman basis functions of the quadratic distance monomials. The Legendre dual of a quadratic distance function d M (p, q) = (pq) T M(p-q) induced by positive-definite matrix M is the quadratic distance d M -1 . (Matrix M is itself usually obtained as the inverse of a variance-covariance matrix Σ in so-called Mahalanobis distances.)

Bregman triangulations

Consider the Euclidean Voronoi diagram vor(S) of a finite set S of points of R d (called sites). Let f be a face of vor(S) that is the intersection of k d-cells of vor(S). We associate to f a dual face f * , namely the convex hull of the sites associated to the subset of cells. If no subset of d + 2 sites lie on a same sphere, the set of dual faces (of dimensions 0 to d) constitutes a triangulation embedded in R d whose vertices are the sites. This triangulation is called the Delaunay triangulation of S, noted del(S). The correspondence defined above between the faces of vor(S) and those of del(S) is a bijection that satisfies: f ⊂ g ⇒ g * ⊂ f * . We say that del(S) is the geometric dual of vor(S). See Figure 12.

A similar construct is known also for power diagrams. Consider the power diagram of a finite set of balls of R d . In the same way as for Euclidean Voronoi diagrams, we can associate a triangulation dual to the power diagram of the balls. This triangulation is called the regular triangulation of the balls. The vertices of this triangulation are the centers of the balls whose cell is non empty.

We derive two triangulations from Bregman Voronoi diagrams. One has straight edges and captures some important properties of the Delaunay triangulation. However, it is not always the geometric dual of the corresponding Bregman Voronoi diagram. The other one has curved (geodesic) edges and is the geometric dual of the Bregman Voronoi diagram.

Bregman Delaunay triangulations

Let Ŝ be the lifted image of S and let T be the lower convex hull of Ŝ, i.e. the collection of facets of the convex hull of Ŝ whose supporting hyperplanes are below Ŝ. We assume in this section that S is in general position if there is no subset of d + 2 points lying on a same Bregman sphere. Equivalently (see Lemma 5), S is in general position if no subset of d + 2 points pi lying on a same hyperplane.

INRIA

Under the general position assumption, each vertex of H = ∩ i H ↑ pi is the intersection of exactly d + 1 hyperplanes and the faces of T are all simplices. Moreover the vertical projection of T is a triangulation del F (S) = Proj X (T ) of S embedded in X ⊆ R d . Indeed, since the restriction of Proj X to T is bijective, del F (S) is a simplicial complex embedded in X . Moreover, since F is convex, del F (S) covers the (Euclidean) convex hull of S, and the set of vertices of T consists of all the pi . Consequently, the set of vertices of del F (S) is S. We call del F (S) the Bregman Delaunay triangulation of S (see Fig. 13). When F (x) = ||x|| 2 , del F (S) is the Delaunay triangulation dual to the Euclidean Voronoi diagram. This duality property holds for symmetric Bregman divergences (via polarity) but not for general Bregman divergences.

We say that a Bregman sphere σ is empty if the open ball bounded by σ does not contain any point of S. The following theorem extends a similar well-known property for Delaunay triangulations whose proof (see, for example [START_REF] Boissonnat | Algorithmic Geometry[END_REF]) can be extended in a straightforward way to Bregman triangulations using the lifting map introduced in Section 3.2.

Theorem 12

The first-type Bregman sphere circumscribing any simplex of del F (S) is empty. del F (S) is the only triangulation of S with this property when S is in general position.

Several other properties of Delaunay triangulations extend to Bregman triangulations. We list some of them.

Theorem 13 (Empty ball) Let ν be a subset of at most d + 1 indices in {1, . . . , n}. The convex hull of the associated points p i , i ∈ ν, is a simplex of the Bregman triangulation of S iff there exists an empty Bregman sphere σ passing through the p i , i ∈ ν.

The next property exhibits a local characterization of Bregman triangulations. Let T (S) be a triangulation of S. We say that a pair of adjacent facets f 1 = (f, p 1 ) and f 2 = (f, p Let S be a given set of points, del F (S) its Bregman triangulation, and T (S) the set of all triangulations of S. We define the Bregman radius of a d-simplex τ as the radius noted r(τ ) of the smallest Bregman ball containing τ . The following result is an extension of a result due to Rajan for Delaunay triangulations [START_REF] Rajan | Optimality of the Delaunay triangulation in R d[END_REF].

Theorem 15 (Optimality) We have del F (S) = min T ∈T (S) max τ ∈T r(τ ).

The proof mimics Rajan's proof [START_REF] Rajan | Optimality of the Delaunay triangulation in R d[END_REF] for the case of Delaunay triangulations.

Bregman geodesic triangulations

We have seen in Section 4.3 that the Bregman Voronoi of a set of points S is the power diagram of a set of balls B centered at the points of S (Theorem 8). Write reg F (B ) for the dual regular triangulation dual to this power diagram. This triangulation4 is embedded in X and has the points of S as its vertices (see Figure 14). The image of this triangulation by ∇ -1 F is a curved triangulation whose vertices are the points of S. The edges of this curved triangulation are geodesic arcs joining two sites (see Section 3.3). We call it the Bregman geodesic triangulation of S, noted del F (S) (see Figure 15). Observe that del F (S) is, in general, distinct from del F (S), the Bregman Delaunay triangulation introduced in the previous section. However, when the divergence is symmetric, both triangulations are combinatorially equivalent and dual to the Bregman Voronoi diagram of S. Moreover, they coincide exactly when F is the squared Euclidean distance.

Applications

In this section, we give some applications related to computational geometry and machine learning.

Centroidal Bregman Voronoi diagrams and Lloyd quantization

Let D be a domain of X and p(x) be a density function defined over D. We define the Bregman centroid of D as the point c * ∈ D such that c * = argmin c∈D x∈D p(x)D F (x||c) dx.

The following lemma states that the mass Bregman centroid of D is uniquely defined and independent of F .

Lemma 11

The Bregman centroid of D coincides with the mass centroid of D.

Proof:

∇ c x∈D p(x) D F (x||c) dx = ∇ c x∈D p(x) (F (x) -F (c) -x -c, ∇F (c) )dx = - x∈D p(x) ∇ 2 F (c)(x -c)dx = -∇ 2 F (c)( x∈D p(x) xdx -c x∈D p(x) dx). Hence, c * = R x∈D p(x) xdx R
x∈D p(x) dx . When x is a random variable following the probability density p(x), x∈D p(x) D F (x||c) dx is called the distortion rate associated to the representative c, the optimal distortion-rate function x∈D p(x) D F (x||c * ) dx is called the Bregman information, and c * is called the Bregman representative. The above result states that the optimal distortion rate exists and does not depend on the choice of the Bregman divergence, and that the Bregman representative c * is the expectation E(x) of x. This result extends an analogous result in the discrete case (finite point sets) studied in [6].

Computing a centroidal Bregman Voronoi diagram of k points can be done by means of Lloyd's algorithm [START_REF] Lloyd | Least squares quantization in PCM[END_REF]. We select an initial set of k points. Then, we iteratively compute a Bregman Voronoi diagram and move the sites to the Bregman centroids of the corresponding cells in the diagram. Upon convergence, the output of the algorithm is a local minimizer of f

((p i , V i ), i = 1, . . . , k) = k i=1 x∈Vi D F (x||p i ) dx , where {p i } k
i=1 denotes any set of k points of X and {V i } k i=1 denotes any tesselation of X into k regions. See [START_REF] Du | Centroidal Voronoi tesselations: Applications and algorithms[END_REF] for a further discussion and applications of centroidal Voronoi diagrams.

ε-nets

Lloyd's algorithm intends to find a best set of k points for a given k so as to minimize a least-square criterion. Differently, we may want to sample a compact domain D ⊂ X up to a INRIA given precision while minimizing the number of samples. Instead of a least-square criterion, we define the error associated to a sample P as error(P ) = max x∈D min pi∈P D F (x||p i ). A finite set of points P of D is an ε-sample of D iff error(P ) ≤ ε.

An ε-sample P is called an ε-net if it satisfies the sparsity condition: max(D F (p||q), D F (q||p)) > ε for any two points p and q in P .

We will see how to construct an ε-net. For simplicity, we assume in the rest of the section that D is a convex polytope. Extending the results to more general domains is possible.

Let P ⊂ D, vor F (P ) be the Bregman Voronoi diagram of P and vor F |D (P ) be its restriction to D. Write V for the set of vertices of vor F |D (P ). V consists of vertices of vor F (P ) and intersection points between the edges of vor F (P ) and the boundary of D. The following lemma states that error(P ) can be computed by examining only a finite number of points, namely the points of V .

Lemma 12 error(P ) = max v∈V min pi∈P D F (x||p i ).

Proof: Let x ∈ D, p x the point of P closest to x and V x the associated cell of vor F |D (P ) (which contains x). V x is a bounded polytope whose vertices belongs to V . Let w be the vertex of V x most distant from p x . We have D F (x||p x ) ≤ D F (w||p x ). This is a consequence of the convexity of F and of the fact that D F (x||p) is measured by the vertical distance between x and H p (Lemma 1).

An ε-net of D can be constructed by the following greedy algorithm originally proposed by Ruppert in the context of mesh generation [START_REF] Ruppert | A Delaunay refinement algorithm for quality 2-dimensional mesh generation[END_REF]. See also [START_REF] Eldar | The farthest point strategy for progressive image sampling[END_REF]. We initialize the sample set P 0 with d points of D lying at distance greater than ε from one another. Then, at each step, the algorithm looks for the point v i of D that is the furthest (for the considered Bregman divergence) from the current set of samples P i . By Lemma 12, this step reduces to looking at the vertices of vor F |D (P i ). If D F (x||v i ) ≤ ε, the algorithm stops. Otherwise, we take v i as a new sample point, i.e. p i+1 = v i , we update the set of sample points, i.e. P i+1 = P i ∪ {p i+1 }, and insert p i+1 in the Bregman Voronoi diagram of the sample points. Upon termination, the set of sample points P t satisfies the hypothesis of Lemma 12 and therefore P t is an ε-sample of D. Moreover, for any two points p and q of P t , we have D F (p||q) > ε or D F (q||p) > ε, depending on whether p has been inserted after or before q. Indeed, we only insert a point if its divergence to the points of the current sample is greater than ε. Hence, P t is an ε-net of D.

To prove that the algorithm terminates, we need the following lemma. Given a Bregman ball B(c, r), we define the biggest Euclidean ball EB(c, r ) contained in B(c, r) and the smallest Euclidean ball EB(c, r ) containing B(c, r).

Lemma 13 Let F be a strictly convex function of class C 2 , there are constants γ and γ (that do not depend on c nor on r) such that r 2 ≥ γ r and r 2 ≤ γ r.

Proof: According to Taylor's formula, there exists a point t of the open segment xc such that

F (x) = F (c) + x -c, ∇F (c) + 1 2 (x -c) T ∇ 2 F (t)(x -c).
Hence,

D F (x||c) = F (x) -F (c) -x -c, c = 1 2 (x -c) T ∇ 2 F (t)(x -c), (14)
where t is a point of the open segment xc.

Since F is strictly convex, the Hessian matrix is positive definite (i.e., x T ∇ 2 F (t)x > 0 for all x in X ), and the domain D being compact, there exist two constants η and η such that, for any y

∈ D, 0 < η ≤ ||∇ 2 F (y)|| ≤ η . If x -c 2 > 2r
η (Fröbenius matrix norm), we deduce from Equation ( 14) that D F (x||c) > r. Therefore, B(c, r) ⊂ EB(c, 2r η ). If xc 2 ≤ 2r η , we have using again Equation ( 14)

D F (x||c) ≤ η 2 x -c 2 ≤ r.
Therefore, EB(c, 2r η ) ⊂ B(c, r). Let p and q be two points such that D F (p||q) = r. Observing that EB(p, r ) ⊆ EB(p, pq ) ⊆ EB(p, r ), we deduce from the above lemma that γ r ≤ pq ≤ γ r

and

γ γ D F (p||q) ≤ D F (q||p) ≤ γ γ D F (p||q).
Another consequence of the lemma is that the volume of any Bregman ball of radius at least r > 0, is bounded away from 0 (when F is of class C 2 ). Hence, since D is compact, the algorithm cannot insert infinitely many points and therefore terminates. Moreover, the size of the sample output by the algorithm can be bounded, as stated in the next lemma. Write D ≤ε = {x| ∃y ∈ D, xy ≤ ε}.

Lemma 14 If F is of class C 2 , the algorithm terminates. If P t denotes the final set of sample points, we have

|P t | = O vol(D) d/2
. Proof: We have already shown that the algorithm terminates. Let P t be the set of points that have been inserted by the algorithm, excluding the initial set (of constant size). Let τ (x) = inf{r : |EB(x, r) P t | ≥ 2} and B p = EB(p, τ (p) 2 ), p ∈ P t . It is easy to see that τ is 1-Lipschitz and that the Euclidean balls B p , p ∈ P t are disjoint. Let q be a point of P t closest to p : τ (p) = pq and, as noticed above, max(D F (p||q), D F (q||p)) > ε. Eq. 15 then implies that τ (p) = pq ≥ √ γ ε. Consider now the midpoint m of pq and write t for the point of P t that minimizes D F (m||.). Since D is convex, m ∈ D and, according to the definition of q, mp ≤ mt . Eq. 15 and the fact that P t is an ε-sample of D then yield mt ≤ √ γ ε. In summary, we have

γ ε ≤ τ (p) = p -q ≤ 2 γ ε. (16)

INRIA

The right inequality shows that all the balls B p , p ∈ P t , are contained in D ≤η where η = √ γ ε. We can now bound the size of P t .

D ≤η dx τ d (x) ≥ p∈Pt Bp∩D ≤η dx τ d (x) (the balls Bp have disjoint interiors) ≥ p∈P vol(Bp∩D ≤η ) ( 3 2 τ (p)) d (τ (x) ≤ τ (p) + p -x ≤ 3 2 τ (p)) ≥ C 3 d |P t | where C = π p p! if d = 2p and C = 2 2p-1 (p-1)! π p-1 (2p-1)! if d = 2p -1.
Using again the Lipschitz property of τ and Eq 16, we have for all

x ∈ B p τ (x) ≥ τ (p) -x -p ≥ 1 2 τ (p) ≥ 1 2 γ ε
We deduce

|P t | ≤ 6 √ γ d 1 Cε d/2 D ≤η dx.
A geometric object O is said α-fat [7] if the ratio r + r -of the radius r + of the smallest ball enclosing O over the radius r -of the largest ball inscribed in O is bounded by α: Thus, since Bregman balls are fat on compact domains, we can build efficient datastructures for point location with applications to piercing (geometric 0-transversal) and others, as described in [START_REF] Efrat | Dynamic data structures for fat objects and their applications[END_REF].

VC-dimension, classification and learning

Some important classification rules rely on Voronoi diagrams; furthermore, the analysis of classification rules (complexity or statistical generalization) sometimes makes use of concepts closely related to Voronoi diagrams. Extending the rules and analyses to arbitrary Bregman divergences, with important related consequences (such as the eventual lost of convexity) is thus particularly interesting for classification, and we review here some notable consequences.

In supervised classification, we are generally interested in capturing the joint structure of X and a set of classes, {0, 1} in the simplest case. For this objective, we build representations of concepts, i.e. functions that map X to the set of classes. A concept class H is a set of concept representations h : X → {0, 1}; for example, should h be a Bregman ball, it would classify 0 the points outside the ball, and 1 the points inside. Armed with these definitions, our supervised classification problem becomes the following one. A so-called target concept, c, which is unknown, labels the points of X ; we have access to its labeling throughout a sampling process: we retrieve examples (i.e., pairs (x, c(x))), independently at random, according to some unknown but fixed distribution D over the set {(x, c(x)) : x ∈ X }. The question is: what are the conditions on H that guarantee the possibility to build, within reasonable time, some h ∈ H agreeing as best as possible with c, with high probability? While the complexity requirement is usual in computer science, the fact that we require adequacy with high probability better than systematically is also a necessary requirement, as there is always the possibility of an extremely bad sampling that would prevent any efficient learning (e.g. we have drawn the same example all the time). In general, rather than directly sampling the domain, we work with a finite data set S of examples which is supposed to be sampled this way.

From the statistical standpoint, learning requires to find a good balance between the accuracy, i.e. the goodness-of-fit of h as measured on S, and the capacity of H, i.e. its ability to learn (or fit in generalization) the data with the smallest number of errors. Consider for example geometric figures in the plane and the "square" concept. Intuitively, an H with too large capacity is like the person who picks a huge quantity of geometric figures including squares, memorizes each of them, and then rejects every square that would not exactly be in its collection (edge lengths, colors, etc.). An H with too little capacity is like the lazy person who keeps as sole concept the fact that squares have four edges. Both extremal situations mean little generalization capabilities, but for different reasons.

There have been intensive lines of works on the measures of this capacity, and one of the most popular is the VC-dimension [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF]. Informally, the VC-dimension of H is the size of the largest dataset S for which H shatters S, i.e. for which H contains all the classifiers that could perform any of the 2 |S| possible labelings of the data. To be more formal, let Π H (S) = {(h(p 1 ), h(p 2 ), ..., h(p n )) | h ∈ H} denote the set of all distinct tuples of labels on S that can be performed by elements of H. While it always holds that Π H (n) ≤ 2 n , the maximal n for which Π H (n) = 2 n is the VC-dimension of H, VCdim(H). The importance of the VC-dimension comes from the fact that it allows to bound the behavior of the empirical optimal classifier in a distribution-free manner [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF]. In particular, if the VC-dimension is finite, the average error probability of the empirical optimal classifier tends to 0 when the INRIA size of the training data set increases. The following lemma proves that the VC-dimension of Bregman balls is the same as for linear separators, and this does not depend on the choice of F . Proof: We use the lifting map introduced in Section 3.2. Given a set S of points in R d , we lift them onto F, obtaining Ŝ ∈ R d+1 .

Let B F be a Bregman ball and write σ for the Bregman sphere bounding B F . From Lemma 5, we know that, for any p ∈ R d , p ∈ B iff p ∈ H ↓ σ . For a given function F , let B F denote the set of all Bregman balls, and let H F denote the set of all lower halfspaces of R d+1 . It follows from the observation above that B shatters S iff H shatters Ŝ. Hence the VC dimension of B over the sets of points of R d is equal to the VC dimension of H over the sets of points of F ⊂ R d+1 .

Since the points of Ŝ are in convex position, they are shattered by H iff the affine hull of their convex hull is of dimension strictly less than the dimension of the embedding space, i.e. d + 1, which happens iff |S| < d + 2. Indeed otherwise, the subset of vertices of any facet of the upper convex hull of Ŝ cannot be obtained by intersecting Ŝ with a lower halfspace (an upper halfspace would be required). Hence, the VC dimension of Bregman balls is at most d + 1.

It is exactly d + 1 since any set of d + 1 points on F in general position generates a ddimensional affine hull A that cannot be shattered by less than d + 1 hyperplanes of A. The same result plainly holds for hyperplanes of R d+1 since we can associate to each hyperplane h of A a hyperplane H of R d+1 such that h = H ∩ A.

This result does not fall into the general family of VC bounds for concept classes parameterized by polynomial-based predicates [START_REF] Goldberg | Bounding the Vapnik-Chervonenkis dimension of concept classes parameterized by real numbers[END_REF], it is mostly exact, and it happens not to depend on the choice of the Bregman divergence. This has a direct consequence for classification, which is all the more important as Bregman balls are not necessarily convex (see Figure 5). Because the capacity of Bregman balls is not affected by the divergence, if we fit this divergence in order to minimize the empirical risk (risk estimated on S), then there is an efficient minimization of the true risk (risk estimated on the full domain X ), as well. There is thus little impact (if any) on overfitting, one important pitfall for classification, usually caused by over-capacitating the classifiers by tuning too many parameters. Some applications of our results in supervised learning also meet one of the oldest classification rule: the k-Nearest Neighbors (k-NN) rule [START_REF] Fix | Discrimatory analysis, nonparametric discrimination[END_REF], in which a new observation receives the majority class among the set of its k nearest neighbors, using e.g. k-order Voronoi diagrams of S (Section 4.4). Various results establish upperbounds for the k-NN rule that depend on the Bayes risk (the true risk of the best possible rule) [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF]. The choice of the proximity notion between observations (it is often not a metric for complex domains) is crucial: if it is too simple or oversimplified, it degrades the k-NN results and may even degrade Bayes risk as well; if it is too complicated or complexified, it may degrade the test results via the capacity of the rule. Searching for accurate "distance" notions has been an active field of research in machine learning in the past decade [START_REF] Wilson | Improved heterogeneous distance functions[END_REF]. Our results on the linearity of the Bregman Voronoi diagrams essentially show that we can mix arbitrary Bregman divergences for heterogenous data (mixing binary, real, integer values, etc.) without losing anything from the capacity standpoint.

Range spaces of finite VC-dimensions have found numerous applications in Combinatorial and Computational Geometry. We refer to Chazelle's book for an introduction to the subject and references wherein [START_REF] Chazelle | The Discrepancy Method[END_REF]. In particular, Brönnimann and Goodrich [START_REF] Brönnimann | Optimal set covers in finite VC-dimension[END_REF] have proposed an almost optimal solution to the disk cover algorithm, i.e. to find a minimum number of disks in a given family that cover a given set of points. Theorem 17 allows to extend this result to arbitrary Bregman ball cover (see also [START_REF] Even | Hitting sets when the VC-dimension is small[END_REF]).

Conclusion

We have defined the notion of Bregman Voronoi diagrams and showed how these geometric structures are a natural extension of ordinary Voronoi diagrams. Bregman Voronoi diagrams share with their Euclidean analogs surprisingly similar combinatorial and geometric properties. We hope that our results will make Voronoi diagrams and their relatives applicable in new application areas. In particular, Bregman Voronoi diagrams based on various entropic divergences are expected to find applications in information retrieval (IR), data mining, knowledge discovery in databases, image processing (e.g., see [START_REF] Inaba | Geometric clustering models for multimedia databases[END_REF]). The study of Bregman Voronoi diagrams raises the question of revisiting computational geometry problems in this new light. This may also allow one to tackle uncertainty ('noise') in computational geometry for fundamental problems such as surface reconstruction or pattern matching.

A limitation of Bregman Voronoi diagrams is their combinatorial complexity that depends exponentially on the dimension. Since many applications are in high dimensional spaces, building efficient data-structures is a major avenue for further research.

Figure 2 :

 2 Figure 2: Visualizing the Bregman divergence. D F (.||q) is the vertical distance between F and the hyperplane tangent to F at q.

Figure 3 :

 3 Figure 3: Legendre transformation of a strictly convex function F : The z-intercept (0, -F * (y )) of the tangent hyperplane H y : z = y , x -F * (y ) of F at ŷ defines the value of the Legendre transform F * for the dual coordinate y = ∇F (y). Any hyperplane passing through an other point of F and parallel to H y necessarily intersects the z-axis above -F * (y ).

Figure 4

 4 Figure 4 depicts several first-type and second-type bisectors for various pairs of primal/dual Bregman divergences.The bisector H F (p, q) for the symmetrized Bregman divergence S F is given by

Figure 4 :Figure 5 :

 45 Figure 4: Bregman bisectors: first-type linear bisector and second-type curved bisector are displayed for pairs of primal/dual Bregman divergences: (a) exponential loss/unnormalized Shannon entropy, (b) logistic loss/dual logistic loss, and (c) self-dual Itakura-Saito divergence. The grid size of R 2 in X and X is ten ticks per unit. First-type (primal linear/dual curved) and second-type (primal curved/dual linear) bisectors are respectively drawn in red and blue.

Figure 6 :Lemma 6

 66 Figure 6: Two Bregman circles σ and the associated curves σ obtained by lifting σ onto F. The curves σ are obtained as the intersection of the hyperplane H σ with the convex hypersurface F. 3D illustration with (a) the squared Euclidean distance, and (b) the Itakura-Saito divergence.

Lemma 7 (

 7 Bregman projection) For any p, there exists a unique point x ∈ W that minimizes D F (x||p). We call this point the Bregman projection of p onto W and denote it p W . Proof: If it is not the case, then define x and y two minimizers with D F (x||p) = D F (y||p) = l. Since W is convex, (x+y)/2 ∈ W and, since D F is strictly convex in its first argument (see Section 2.1), D F ((x+y)/2||p) < D F (x||p)/2+D F (y||p)/2. But D F (x||p)/2+D F (y||p)/2 = l yielding a contradiction.

Figure 7 :

 7 Figure 7: Generalized Pythagoras' theorem for Bregman divergences: The projection p W of point p to a convex subset W ⊆ X . For convex subset W, we have D F (w||p) ≥ D F (w||p W ) + D F (p W ||p) (with equality for and only for affine sets W).

Figure 8 :

 8 Figure 8: Bregman bisectors and their relationships with respect to Λ(p, q) (straight line segments) and Γ F (p, q) (bold curves), for the Itakura-Saito divergence (left) and I-divergence (right). Bold curves become linear in X ; colors depict the Bregman orthogonality relationships of Lemma 8.

1 λ=0D

 1 F (p 0 ||p λ )dλ (respectively 1 λ=0 D F (p λ ||p 0 )dλ) over all curves C = {p λ } 1 λ=0 . Proof: For any curve C between p 0 and p 1 , we measure the Γ length as Γ (C) = λ D F (p λ ||p 0 )dλ. Fix some inner point p ∈ Γ F (p 0 , p 1 )\{p 0 , p 1 }. From the three-point property (Property 5), the set of points {y ∈ X | D F (y||p 0 ) = D F (y||p) + D F (p||p 0 )} is the hyperplane H p : y, h = b (h is a perpendicular vector to the hyperplane) which splits X into two open half-spaces H + p : y, h > b, and

1 λ=0D

 1 F (p 1 ||p λ )dλ (respectively, minimizes also1 λ=0 D F (p λ ||p 1 )dλ) over all curves C = {p λ } 1 λ=0 .

Corollary 2

 2 Suppose p(.|θ 0 ) and p(.|θ 1 ) are probability density functions of the same exponential family E F . Then Γ F (θ 0 , θ 1 ) (resp. Λ(θ 0 , θ 1 )) minimizes Γ (C) = 1 λ=0 KL(θ 0 ||θ λ )dλ (resp. Λ (C) = 1 λ=0 KL(θ λ ||θ 0 )dλ) over all curves C = {p(.|θ λ )} 1 λ=0 .

Because

  Bregman divergences are not necessarily symmetric, we associate to each site p i two types of distance functions, namely D i (x) = D F (x||p i ) and D i (x) = D F (p i ||x). The minimization diagram of the D i , i = 1, . . . , n, is called the first-type Bregman Voronoi diagram of S, which we denote by vor F (S). The d-dimensional cells of this diagram are in 1-1 correspondence with the sites p i and the d-dimensional cell of p i is defined as vor

Figure 9 :

 9 Figure 9: Three types of Bregman Voronoi diagrams for (a) the Kullback-Leibler and (b) the Itakura-Saito divergences. First-type affine Bregman Voronoi diagram (red), second-type Bregman Voronoi diagram (blue) and symmetrized Bregman Voronoi diagram (green).

Figure 9

 9 Figure 9 illustrates these Bregman Voronoi diagrams for the Kullback-Leibler and the Itakura-Saito divergences. Note that the Euclidean Voronoi diagram is a Bregman Voronoi diagram since vor(S) = vor F (S) = vor F (S) for F (x) = x 2 .For asymmetric Bregman divergences D F , we can further consider the symmetrized Bregman divergence S F = D F and define a third-type Bregman Voronoi diagram vor F (S). The cell of vor F (S) associated to site p i is defined as:

Lemma 10

 10 vor F (S) = ∇ -1 F (vor F * (S )) and vor F (S) = ∇ -1 F (vor F * (S )). Proof: By Lemma 3, we have D F (x||y) = D F * (y ||x ), which gives vor F

3 :

 3 S F (p, q) = D F (p||q) = D F (q||p) where F is a 2d-variate function and x = [x x ] T . As already noted in §2.3, x lies on the d-manifold X = {[x x ] T | x ∈ R d }. It follows that the symmetrized Voronoi diagram vor F (S) is the projection of the restriction to X of the affine diagram vor F ( S) of R 2d where S = {p i , p i ∈ S}. Hence, computing the symmetrized Voronoi diagram of S reduces to: 1. computing the first-type Bregman Voronoi diagram vor F ( S) of R 2d , 2. intersecting the cells of this diagram with the manifold X , and 3. projecting all points of vor F ( S) ∩ X to X by simply dropping the last d coordinates.

d 2 )

 2 and can be computed in optimal-time Θ(n log n + n d 2 ) for any fixed dimension d. From Theorem 6 and Lemma 10, we then deduce the following theorem. Theorem 7 The Bregman Voronoi diagrams of type 1 or 2 of a set of n d-dimensional points have complexity Θ(n d+1 2

Figure 10 :

 10 Figure 10: Voronoi diagrams as minimization diagrams. The first row shows minimization diagrams for the Euclidean distance and the second row shows minimization diagrams for the Kullback-Leibler divergence. In the first column, the functions are the non-linear functions D i (x) and, in the second column, the functions are the linear functions H pi (x), both leading to the same minimization diagrams. Isolines are shown in green.

INRIA 4 . 3

 43 Bregman Voronoi diagrams from power diagramsThe power distance of a point x to a Euclidean ball B = B(p, r) is defined as ||p -x|| 2 -r 2 . Given n balls B i = B(p i , r i ), i = 1, . . . , n, the power diagram (or Laguerre diagram) of the B i is defined as the minimization diagram of the corresponding n functions

INRIAFigure 11 :

 11 Figure 11: Affine Bregman Voronoi diagrams (left column) can be computed as power diagrams (right column). Illustrations for the squared Euclidean distance (a), Kullback-Leibler divergence (b), and exponential divergence (c). Circles are drawn either in grey to denote positive radii, or in red to emphasize imaginary radii. Observe that although some cells of the power diagrams may be empty, all cells of the affine Bregman Voronoi diagram are necessarily non-empty.

Theorem 10

 10 The k-order Bregman Voronoi diagram of n d-dimensional points is a weighted Bregman Voronoi diagram.

  for all j or, equivalently, iff x belongs to the cell of c i in the weighted Bregman Voronoi diagram of the c i . k-bag Bregman Voronoi diagrams Let F 1 , ..., F k be k strictly convex and differentiable functions, and α = [α 1 ... α k ] T ∈ R k + a vector of positive weights. Consider the d-variate function F α = k l=1 α l F l . By virtue of the positive additivity property rule of Bregman basis functions (Property 3), D Fα is a Bregman divergence. Now consider a set S = {p 1 , ..., p n } of n points of R d . To each site p i , we associate a weight vector α i = [α

Theorem 11

 11 The k-bag Voronoi diagram (for k > 1) on a bag of d-variate Bregman divergences of a set of n points of R d has combinatorial complexity O(n k+d 2

Figure 12 :

 12 Figure 12: Ordinary Voronoi diagram (red) and geometric dual Delaunay triangulation (blue).

Figure 13 :

 13 Figure 13: Bregman Delaunay triangulation as the projection of the convex polyhedron T .

  2 ) of T (S) is regular iff p 1 does not belong to the open Bregman ball circumscribing f 2 and p 2 does not belong to the open Bregman ball circumscribing f 1 (the two statements are equivalent for symmetric Bregman divergences). Theorem 14 (Locality) Any triangulation of a given set of points S (in general position) whose pairs of facets are all regular is the Bregman triangulation of S.

Theorem 16

 16 The Bregman geodesic triangulation del F (S) is the geometric dual of the 1sttype Bregman Voronoi diagram of S. Proof: We have, noting * ≡ for the dual mapping, and using Theorem 8 vor F (S) ≡ pow(B ) * ≡ reg(B ) = ∇F (del F (S)).

Figure 14 :Figure 15 :

 1415 Figure 14: First-type Kullback-Leibler Bregman Voronoi diagram (a) obtained from the corresponding power diagram (b), and its associated dual regular triangulation rooted at gradient vertices (blue).

Corollary 3 For C 2 = O( 1 )

 321 r + r -≤ α. Euclidean balls are therefore 1-fat, namely the fattest objects. It has been shown that considering the fatness factor for a set of objects yields in practice efficient tailored data-sensitive algorithms[7] by avoiding bad configurations of sets of skinny objects. A direct consequence of Lemma 13 is that Bregman balls (in fixed dimensions) are fat (i.e., α = O(1)) on any compact domain: Bregman generator functions, Bregman balls on any compact domain are fat.Proof: Indeed, consider any Bregman ball defined on a compact domain for a C 2 strictly convex and differentiable Bregman generator function F . Its fatness α is upper bounded by γ γ , where γ and γ are the two constants (depending on F and D) of Lemma 13. Recall that Lemma 13 considers concentric Euclidean balls ham sandwiching a Bregman ball, all centered at position c. We have α ≤ r + r -≤ r + since r - c ≤ r -and r + c ≥ r + , where r + c (respectively, r - c ) denote the radius of the smallest enclosing (respectively, largest inscribed) Euclidean ball centered at c. The fatness property simply means that we can cover any Bregman ball by a constant number of (convex) Euclidean balls.

Theorem 17

 17 The VC dimension of the class of all Bregman balls B F of R d (for any given strictly convex and differentiable function F ) is d + 1.

  }dx. F is called the cumulant function or the logpartition function. F fully characterizes the exponential family E F while term C(x) ensures density normalization. (That is, p(x|θ) is indeed a probability density function satisfying x p(x|θ)dx = 1.)When the components of the sufficient statistics are affinely independent, this canonical representation is said to be minimal, and the family E F is called a full exponential family of order d = dim X . Moreover, we consider regular exponential families E F that have their support domains topologically open. Regular exponential families include many famous distribution laws such as Bernoulli (multinomial), Normal (univariate, multivariate and rectified), Poisson, Laplacian, negative binomial, Rayleigh, Wishart, Dirichlet, and Gamma distributions. Table

RR n°6154

See Java applet at http://www.csl.sony.co.jp/person/nielsen/BVDapplet/ RR n°6154

Applet at http://www.csl.sony.co.jp/person/nielsen/BVDapplet/ INRIA

http://www.inria.fr ISSN 0249-6399

Acknowledgements

Frédéric Chazal, David Cohen-Steiner and Mariette Yvinec are gratefully acknowledged for their comments on this paper. The work by the second author has been partially supported by the project GeoTopAl (1555) of the Agence Nationale de la Recherche (ANR).