数理解析研究所講究録1093 巻 1999 年 21-26

21Computational Geometry on Statistical Manifolds for Clustering(Extended Abstract)

Mary Inaba*, Hiroshi Imai* and Kunihiko Sadakane*
Abstract
This paper investigates computational-geometric aspects of clustering problems on the statistical manifolds in information geometry. This extends the traditional Euclidean computational
geometry to a new one in the space of information geometry. Computational-geometric concepts
such as Voronoi diagrams in the space of information geometry are described, and then someresults on the geometric clustering problem based on the concepts are outlined.

1

Introduction

In information geometry [1, 2, 7], differential geometric properties of probabilistic distributions
and other stochastic systems have been studied. A set of parametrized probability distributions
form a Riemannian manifold by their parameters, and, the exponential family of probabilistic
distribution is the most typical and-behaved family. This family contains normal distribution,
Poisson, finite discrete, and exponential distribution as special cases. Furthermore, this has a nicedifferential-geometric property, called “
flat,” and a divergence in general form, which is adistance-like function between two probability distributions and contains both
-Leiblerdivergence and squared Euclidean distance as its special cases $[1, 2]$ . The Voronoi diagram in theflat space is studied in [8, 9, 10].
In this paper, we extend our results for the Euclidean clustering problem [6] to clustering
problems by divergence on statistical manifolds. Clustering problem is to group similar objects
under some criteria, and, in general it is-hard. Geometric -clustering problem is to find agood partition,a -clustering, of the given set of points $p_{i}=(x_{i})(i=1, \ldots, n)$ in the-dimensional space into
nonempty subsets. We first introduce the weighted
Voronoi diagram by divergence, and analyze its complexity. This generalized Voronoi diagram
share nice properties with the Euclidean diagrams. Our Voronoi diagram by divergence is thenapplied to the general mixture clustering case. With this unified approach via the divergence, we
could characterize an optimal clustering for the pure variance criterion in the Euclidean case which
was left open in [6]. We propose a randomized 2-clustering approximation algorithm using random
sampling technique, together with analyzing its approximation performance.
$\mathrm{w}\mathrm{e}\mathrm{U}$$\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{U}\mathrm{y}$$\mathrm{K}\mathrm{u}\mathrm{U}\mathrm{b}\mathrm{a}\mathrm{C}\mathrm{k}$$\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{U}\mathrm{y}$$k$$\mathrm{N}\mathrm{P}$$\mathrm{c}\mathrm{a}\mathrm{U}\mathrm{e}\mathrm{d}$

$d$

$k$

$S$

$k\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{i}^{\text{・}}\mathrm{n}\mathrm{t}$$S_{1},$$n$$\ldots,$$S_{k}$$\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{o}\mathrm{r}\mathrm{e}\mathrm{t}\mathrm{i}_{\mathrm{C}}\mathrm{a}\mathrm{U}\mathrm{y}$2Statistical manifolds of probability distributions

For the statistical estimation, in a traditional form, first we assume an underlying distribution
such as normal distribution or Poisson distribution, then, from a set of observed data, we estimate
the parameters of the distribution, such as mean or deviation in the normal distribution case.
In this sense, once distribution has been assumed, statistical estimation can be regarded as theestimation of parameter of the distribution. Here, we regard a statistical distribution characterizedby parameters, as a point in the -dimensional parametric space, geome.trically structures by theproperties of the distributions.
$d$

$d$

*Department of Information Science, University of Tokyo, Hongo, Bunkyo-ku, Tokyo, 113-0033 Japan22
by their paA set of parametrized probabihty distributions form a Riemanmian manifold
rameters. For example, a class of one-dimensional normai distribution with mean and standard
deviation form a manifold $\lambda 4=\{[\mu, \sigma]|\sigma>0\}$ , the upper half plane. This section describes
fundamental properties of this manifold for a wide and well-behaved class of probabihty distributions, called the exponential family. Since we will use two dual coordinates, -coordinate anduse the tensor-coordinate, which generalizes the polarity with respect to a paraboloid, wenotation.
$\mathcal{M}$$\mu$$\sigma$$\theta$$\mathrm{w}\mathrm{i}\mathrm{U}$$\eta$

2.1

Exponential family

belongs to the exponential family if its probaA probability distribution parametrized bybihty density function $f(x;\theta)$ with probability variable (vector) is expressed as$\theta=[\theta^{i}]$$x$$f(x; \theta)=\exp[C(x)+\sum_{i}.\theta^{i}F_{1(}.x)-\psi(\theta)]$Since$\int f(x;\theta)\mathrm{d}_{X}=1,$$\psi$.

is given by$\psi(\theta)=\log\int\exp[C(X)+\sum\theta^{i}F_{1}.(x)]\mathrm{d}_{X}i$

For this$\theta=[\theta^{\}]$, we define$\eta=[\eta_{i}]$by$\eta_{i}=\int F_{i}(x)f(_{X};\theta)\mathrm{d}X$and are two coordinate systems on the manifold
exponential family. is also given by$\mathcal{M}$$\theta$$\eta$.

of parameters of the distributions in the$\eta$$\eta_{i}=\frac{\partial\psi(\theta)}{\partial\theta^{i}}$In the case of the exponential family, the dual potential function
system by$\varphi(\eta)$

is defined in the -coordinate$\varphi(\eta)=\int f(x;\theta)(\log f(X;\theta)-^{o())}X\mathrm{d}_{X}$in the right-hand side is that corresponding to in the left-hand side. Note that when$C(x)\equiv 0$ , this potential function
becomes the minus of entropy of distribution. 9 is then given
by

where

$\theta$$\eta$$\varphi$$\theta^{i}=\frac{\partial\varphi}{\eta_{i}}$.

of points .
and $\eta=\eta(p)$ give two coordinate systems on the manifold
In fact,
Examples of the exponential family are normal distributions, finite discrete distribution, exponential distribution, etc., see $[1, 2]$ .
$\mathcal{M}$$\theta=\theta(p)$

2.2

$p$

Properties of the divergence

In the sequel, we adopt the Einstein’s notation.
$\theta^{i}\eta_{i}\equiv\sum_{i}\theta i\eta i$

The statistical manifold of the exponential family has nice differential-geometric structure,
flat. In order to investigate the discrete proxinity structure of this manifold, in this paper,for the exponential
it suffices to consider the -coordinate and the -coordinate of the manifold
denote the and -coordinate values
andfamily as two dual affine coordinate systems.
, and $\eta(p)=[\eta_{1}(p), \ldots,\eta_{d}(p)]$ . In the duallyfor a point on , that is,.
flat space, we can define a distance-like function divergence between two points and on$\mathrm{c}\mathrm{a}\mathrm{U}\mathrm{e}\mathrm{d}$$\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{U}\mathrm{y}$$\Lambda 4$$\theta$$\eta$$\theta(p)$$p$$\eta(p)$$\theta-$$\eta$$\theta(p)=[\theta^{1}(p), \ldots , \theta^{d}(p)]$$\mathcal{M}$$p$

Definition 1 (Divergence) Consider the two potential functions, define the divergence $D(p||q)$ by
tial family. For two points$\psi,$ $\varphi:\mathcal{M}arrow \mathrm{R}$$p,$$q\in \mathcal{M}$$D(p||q)=\psi(p)+\varphi(q)-\theta^{i}(p)\eta_{i}(q)$$q$$\mathcal{M}$for the exponen-

23
The pair of potential functions are connected via the Legendre transformation, that is,$\mathit{9}^{i}=\frac{\partial\varphi}{\partial\eta_{i}}$$\psi,$$\varphi$,$\eta_{\dot{2}}=\frac{\partial\psi}{\partial\theta^{i}}$are strictly convex, and$\varphi(q)=\max\{p\in s\theta^{i}(p)\eta i(q)-\psi(p)\}$,$\psi(p)=\mathrm{m}\mathrm{a}\mathrm{x}q\in s^{\{\mathit{9}^{i}(p})\eta i(q)-\varphi(q)\}$Hence, $D(p||q)\geq 0$ , and $D(p||q)=0$ iff $p=q$ . But, unhike the distance, $D(p||q)\neq D(q||p)$ ingeneral.
Next, we consider the relation of $D(p||q)$ with the potential function and a tangent hyperplane.
Add a new coordinate , corresponding to the height, to the -coordinate system, and consider the, lift it up to the graph-space. For,,graph $z=\varphi$ in the, the height difference of a pointand consider the tangent hyperplane. Then, for a point
is given bylifted to the graph$\varphi$$z$$\eta$$[\eta, z]$$(\eta_{1}(p),\eta_{2}(p),$ $\ldots$$p\in \mathcal{M}$$\eta_{d}(p),$ $\varphi(p))$$q\in \mathcal{M}$$z=\varphi(\eta)$$\varphi(q)-\mathit{9}i(p)\eta_{i}(q)+\theta i(p)\eta_{i}(p)-\varphi(p)=\psi(p)+\varphi(q)-\theta^{i}(p)\eta i(q)=D(p||q)$

By the symmetric duality, this holds also in 9-coordinate system, i.e., the divergence $D(p||q)$is also the difference of the height at the point between the potential function and tangent
hyperplane on on the point .
The divergence has such a nice and natural meaning, which was used to analyze the -Voronoi
diagram as stated and cited in Theorem 1.
In the Euclidean case, which corresponds to a self-dual case, $\psi=\varphi=\sum_{i=1}^{d}x_{i}2/2$ and. The divergence is a half of the square of the Euclidean distance. For the exponential family,
the divergence coincides with the Kullback-Leibler divergence $D_{K}(q||p)$ , also known as the relative
flat structure is an extension of the ordinary Euclidean case, and theentropy. Thus, thisdivergence is an extension of the squared Euclidean distance.
Furthermore, the maximum likelihood method in statistical inference can be interpreted in anatural way in the -coordinate system, say, taking the average or the orthogonal projection to
obtain maximum likelihood estimators. For detaik, see $[1, 2]$ .
$\psi$$p$$\psi$$q$$\nabla^{*}$$\theta^{i}=\eta_{i}=$$x_{i}$$\mathrm{d}\mathrm{u}\mathrm{a}\mathrm{U}\mathrm{y}$$\eta$

3

Weighted Voronoi diagrams by divergence

The Voronoi diagram by the divergence is investigated in $[8, 9]$ . In extending the -pair sumof squared Euclidean distances to the divergence case, multiplicatively and additively weighted
beVoronoi diagrams are useful. Hence, this section investigates such weighted diagrams. As
seen, the weighted diagram has similar structures as the weighted Euclidean diagram [3], and thisresult may be viewed as an extension of [3].
We begin with a non-weighted case.
$\mathrm{a}\mathrm{U}$$\mathrm{w}\mathrm{i}\mathrm{U}$Definition 2 ( -Voronoi diagram) For generator points $r^{(j)}(j=1, \ldots, k)$ , thediagram consists of Voronoi regions $V(r^{(j)})$ defined as follows in .
$k$$\nabla^{*}$$\nabla^{*}$-Voronoi$f\mathit{9}J$$V(r^{(j)})=,\cap\{p|D(p^{(}|j)|p)<D(p^{(j’)}||p)\}j\neq j$$V(r^{(j)})(j=1, \ldots, k)$

For the$\nabla^{*}$

partition the manifold, which is called the$\nabla^{*}$-Voronoi diagram.

-Voronoi diagram, the following holds.

-Voronoi diagram can be obtained as the projection to
Theorem 1 (Onishi, Imai [9]) The
hyperplanes
envelope
which are tangent hyperplanes in the
upper
the
the manifoldofof$z=\varphi(p)$.
atcoordinate of the graph$\nabla^{*}$$[\eta, z]-$$\mathcal{M}$$[\eta(p), \varphi(p)]$-Voronoi diagram can be bounded by
By this theorem, the combinatorial complexity of thepolytopes.
the upper bound theorem for convex.
The weighted Voronoi diagram by the divergence is defined as$\nabla^{*}$$\mathrm{f}\mathrm{o}\mathrm{U}_{0}\mathrm{w}\mathrm{S}$24, a multiDefinition 3 (Weighted -Voronoi diagram) Suppose that, for each points
is givenpoint
regiongiven.
VoronoiTheweightareadditive
anofplicative weight$r^{(j)}$$k$$\nabla^{*}$$r^{(j)}$$\tilde{w}^{(j)}$$w^{(j]}\cdot\cdot and$$by$$V(r^{(j)})=j’ \bigcap_{\neq j}\{p|w^{(j}(r|(j)|p))_{D}+\tilde{w}^{(j})<w^{(j’)}D(r^{(j’)}||p)+\tilde{w}\}(j’)$

The collection of $V(r^{(j)})(j=1, \ldots, k)$ is called the weighted
weighted Voronoi diagram.

$\nabla^{*}$-Voronoi diagram, or simply theConcerning the the combinatorial complexity of this weighted diagram, we can show that each, and hence obtain the following.
Voronoi region has complexity$o(n^{\lfloor\frac{d+1}{2}\rfloor})$

Theorem 2 The combinatorial complexity of the weighted divergence Voronoi diagram is$o(n^{\lfloor\frac{d+3}{2}\rfloor})$.

Clustering by divergence4

For a given set of points $p^{(l)}(l=1, \ldots,n)$ on the manifold
whose uniori is .
into nonempty disjoint subsets
$S$

$S$

$n$

$k$

$S_{1},$$\ldots,$$\mathcal{M}$, a -clustering is a partition of
$k$

$S$

$S_{k}$

Problem 1 (Divergence-sum clustering)$S_{\mathrm{j}}( \min_{tj=1,\ldots,k)}\sum j=1(kp\sum D(r|(j)|\iota)\in Sjp\mathrm{t}^{1)})$$(\dot{g}),$Here,is a representative point for , and, since the sum of divergence is minimized at thein the -coordinate. This clustering criterion corcentroid,
is simply set to the centroid of
responds to maximizing the Classification Maximum Likelihood (CLM) for the exponential family[4].
Next, generalizing the weighted Euclidean case, we consider the following.
$r^{(j)}$$S_{j}$$r^{(j)}$$S_{j}$$\eta$

Problem 2 (Multiplicatively weighted divergence-sum clustering)$\dot{\mathrm{m}}\mathrm{n}\sum_{j=1}^{k}w(|sj|)p^{()}\sum_{\iota\in sj}D(r^{(j)}||p)(1)$

where $w(s)=s^{\alpha}$ , and the following cases have their own meanings:$\bullet$$\bullet$$\bullet$$w(s)=1$ (this case corresponds to the simple sum case)-pair sum of squared Euclidean distances in the
Euclidean case)$w(s)=1/s$ (this case corresponds to the variance)$w(s)=s$ (this case corresponds to the$\mathrm{a}\mathrm{U}$for cluster
is simply set to the centroid of
As in the first problem,
We further extend the problem to a mixture case.
$r^{(j)}$$S_{j}$$S_{j}$

in the -coordinate.
$\eta$

bdongs to
Problem 3 (Mixed divergence clustering) In this clustering model, each point
every duster in some sense. Specifically, $\zeta(j,l)(\geq 0)$ denotes how much point belongs to cluster, where$p^{(l)}$$p^{l}$$S_{j}$$\sum_{j=1}^{k}\zeta(j, l)=1$

for each .
$l$Then, the dustering problem is to find an optimd -clustering for$k$$\min\sum_{j=1}^{k}w(\sum_{l=1}^{n}\zeta(j,l))\mathrm{t}=\sum^{n}\zeta(j,l)D(\tilde{r}^{(j)}||p^{(})\iota \mathrm{t})$25
In this problem,$\tilde{r}^{(j)}$is the weighted centroid defined by$\eta_{i}(\tilde{r}^{(j)})=\frac{1}{\sum_{l=1}^{n}\zeta(j,l)}\sum_{=l1}^{n}\zeta(j, l)\eta i(p^{(l}))$Now, we describe our results for these problems. First, we have the following.
Theorem 3 An optimal dustering for the divergence-sum clustering problem is identical with apartition by the-Voronoi diagram generated by the centroids of clusters.
$\nabla^{*}$The number of partitions of points induced by the nabla*-Voronoi diagram is shown to be$O(n^{\min\{d}-d+k-2\}k,dk)$ and furthermore such partitions can be enumerated in time proportional tothe bound in [5]. This number is polynomial when and are regarded as constants, thus revealing
the geometry helps in solving the clustering problem for moderate and .
Concerning the Problem 2, at this point, it seems hard to generalize a result for the Euclideancase in [6] to this general case, but its mixture version of Problem 3 can be again characterized by. Furthermore, this characterization solves unsolved
the weighted-Voronoi diagram as
pure variance problem in the Euclidean case with $w(s)=1/s$ .
To analyze optimal clusterings in the mixed divergence clustering problem, we adopt infinitesimal analysis.
Now, define
by
$n$

$d$

$k$

$d$

$\nabla^{*}$$k$$\mathrm{f}\mathrm{o}\mathrm{U}\mathrm{o}\mathrm{w}\mathrm{s}$$W_{j}$$W_{j}=w( \sum_{l=1}^{n}\zeta(j,l))\iota=\sum\zeta(j,l)D(_{\tilde{\Gamma}}(j)||p^{(l})1n)$also depends on
where we recall
lemma.
We now have the$\tilde{r}^{(j)}$$\zeta(j, l)$.

$\mathrm{f}\mathrm{o}\mathrm{u}_{0}\mathrm{W}\mathrm{i}\mathrm{n}\mathrm{g}$Lemma 1$\frac{\partial W_{j}}{\partial\zeta(j,l)},=w(\iota=\sum_{1}^{n}\zeta(j, l))D(\tilde{r}^{(j)}||p^{(})\iota’)+\frac{\partial w(\sum_{\iota=}^{n}1\zeta(j,l))}{\partial\zeta(j,l)},\sum_{l=1}^{n}\zeta(j, l)D(\tilde{r}^{(j)}||p)(l)$Theorem 4 Suppose $w(s)=s^{\alpha}$ . For an opti.$mal$ solution of the mixed divergence clustering probwith multiplicative weight
lem, consider the weighted -Voronoi diagram for weighted centroids
and additive weight$\tilde{r}^{(j)}$$\nabla^{*}$$w( \sum_{l=1}^{n}\zeta(j, l))$$\alpha(\sum_{l=1}\zeta(jn, l).)^{\alpha}-1\sum_{\iota=1}^{n}\zeta(j,l)D(\tilde{r}^{(}|j)|p)(l)$$\zeta(j, l)=1,$
in the Voronoi region
Then, for point.
it is on some boundary face on the Voronoi region of$of\tilde{r}^{(j)},$$p^{(l)}$.

, for points$an\dot{d}$$p^{(l)}$with $0<\zeta(j, l)<1$ ,$\tilde{r}^{(j)}$As mentioned above, this general characterization solves an unsolved pure variance-clustering
problem in the Euclidean case with $w(s)=1/s$ to which no characterization by the Voronoi diagram
has not yet been known.

5

Random sampling algorithm

technique to divergenceWe extend approximate algorithm for 2-clustering using random
sum problem, based on the algorithm in the Euclidean case [6]. The idea of this algorithm is that,sampled data might not reflect the cost function of whole data, but, sampled data can reflect thepossible partitioning on sampled data, and, compute thecentroid with high probability, so, try
centroid using sampled data, then, compute cost function using whole data and get minimum one.
However, in the Euclidean case, the divergence is directly connected with the variance, theclustering cost function, while in general cases it is not. Hence, as for analysis of approximation
ratio we restrict ourselves to the case of finite discrete distribution. Using the Kullback-Leibler
divergence, we apply the tail distribution analysis using this divergence.
$\mathrm{S}\mathrm{a}\mathrm{m}_{\mathrm{P}^{1\dot{\mathrm{m}}\mathrm{g}}}$$\mathrm{a}\mathrm{U}$26[Randomized 2-clustering algorithm with divergence]by $m$ independent draws at random;
1. Sample a subset of $m$ points
of in the -coordinate system, execute the
2. For every linearly separable 2-clustering
following:$T$$\mathrm{h}\mathrm{o}\mathrm{m}S$$T$$(T_{1},T_{2})$$\eta$

in the -coordinate system, respectively;
Compute the centroids and of andof by dividing by the hyperplane with the same diverFind a 2-clustering-coordinate system,
inthegence between andamong theseand maintain theCompute the value of Cost
values;$t_{2}$$t_{1}$$T_{2}$$S$$(s_{1}, s_{2})$$t_{2}$$t_{1}$$T_{1}$$\eta$$S$$\eta$$(s_{1})+\mathrm{C}\mathrm{o}\mathrm{s}\mathrm{t}(s_{2})$$\mathrm{m}\ddot{\mathrm{m}}\mathrm{m}\mathrm{u}\mathrm{m}$This randomized algorithm is an approximation algorithm, and its approximation ratio may beevaluated, and we have the following.
Theorem 5 Suppose that there is an optimal 2-clustering such that the sizes of each cluster
are within some constant factor to each other. Let $D$ be the minimum among the averages offor each duster in the optimal dustering. Then, for some constant a’ with $\alpha>d>0$ ,in $O(nm^{d})$ time, whose sum of divergences is within
the randomized algonthm finds a$D(\overline{q}(s_{j})||p)(\mathrm{t})$$\mathit{2}-dusie\dot{n}ng$

a factor of $1+c$ with probability at least $1-4d \exp(-2\alpha’(1-\exp(-\frac{\mathrm{c}D}{n})2m)$ .

References[1] S. Amari: Differential Geometrical Methods in Statistics. Lecture Notes in Statistics, Vol.28,
Springer-Verlag, New York, 1985.
[2] S. Amari and H. Nagaoka: Method of Information Geometry (in Japanese). Iwanami Shoten,
Tokyo, 1993.
[3] F. Aurenhammer and H. Edelsbrumer: An Optimal Algorithm for Constructing the Weighted
Voronoi Diagram in the Plane. Pattem Recognition, Vol.17 (1984), 251-257.
[4] G. Celeux and G. Govaert: A Classification EM Algorithm for Clustering and Two Stochastic
Versions. Computational Statistics b Data Analysis, Vol.14 (1992), pp.315-332.
[5] M. Inaba and H. Imai: The Number of Partitions of n Points Induced by the Voronoi Diagram via the Conjugacy Generated by k Points. Japanese-Hungarian Symposium on Discrete
Mathematics and Its Applications, March 1999.

[6] M. Inaba, N. Katoh and H. Imai: Applications of Weighted Voronoi Diagrams and Randomization to Variance-Based -Clustering. Proceedings of the 10th ACM Symposium on Computational Geometry, 1994, pp.332-339.
[7] M. K. Murray and J. W. Rice: Differential Geometry and Statistics. Chapman&HaU, 1993.
$k$[8] K. Onishi and H. Imai: Voronoi Diagrams for an Exponential Family of Probabihity DistribuComputation,
tions in Information Geometry. Japan-Korea Joint Workshop on, 1997, pp.1-8.
[9] K. Onishi and H. Imai: Riemamian Computational Geometry –Voronoi Diagram andDelaunay-type Triangulation in Dually Flat Space, submitted for publication.
$Alg_{\mathit{0}\gamma}.ithmS.and$$\dot{\mathrm{F}}\mathrm{u}\mathrm{k}\mathrm{u}\mathrm{o}\mathrm{k}\mathrm{a}$M. Inaba, F. Takeuchi and K. Imai: Voronoi Diagrams[10] K. Sadakane, H. Imai, K.
by Divergences with Additive Weights. Proceedings of the 14th Annual ACM Symposium on
Computational Geometry, 1998, pp.403-404.
$\mathrm{O}\mathrm{n}\mathrm{i}\mathrm{s}\mathrm{h}\mathrm{i}_{s}$