On-line Evolutionary Exponential Family Mixture

1 Introduction

Clustering is a fundamental problem in machine learning and data mining. Conventional clustering algorithms, such as k-means [Hartigan and Wong, 1979] and spectral clustering [Ng et al., 2002], focus on static data and assume all the data are I.I.D. (Independent and Identically-Distributed) samples from one underlying distribution. However, in lots of dynamic applications, data come from different time epochs. Due to concept drifting or noise varying, the distribution of epoch data often drifts along time. For example, contents under the topic "life style" in a Bulletin Board System (BBS) often differ from those of one year ago while not deviating too much. The clustering task on this kind of data raised the problem of evolutionary clustering [Chakrabarti et al., 2006]. In this case, the final target is to provide a set of partitions, one for each time epoch. In addition, as the data distributions of adjacent epochs are close to each other, the clustering results of epochs should be smooth along time.

We should distinguish evolutionary clustering from incremental clustering [Charikar et al., 1997]. Incremental clustering gives a single partition for all the data, although the data enter into the algorithm sequentially. Two properties are emphasized in incremental clustering, the first is the one-pass manner of the access to data, and the second is the equivalence between the original non-incremental algorithm and the corresponding incremental one.

The necessity of evolutionary clustering lies in two aspects. First, when distribution drifts, applying a conventional clustering algorithm to overall data may not be appropriate. Second, if we apply a conventional clustering algorithm independently to each epoch data, the smoothness of clustering results along time can not be preserved. The second aspect can be realized from two facts. (1) For non-deterministic clustering algorithms relying on initialization, such as k-means, Gussian Mixture Model (GMM), etc., the clustering results of adjacent epochs may be quite different from each other due to local optima, even when the two distributions are almost the same. (2) For deterministic clustering algorithms, such as spectral clustering and agglomerative hierarchical clustering, data noise may lead to different clustering results between adjacent epochs.

Evolutionary clustering can be off-line or on-line 1 . Two off-line methods have been proposed by [Wang et al., 2007] and [Ahmed and Xing, 2008].

The first on-line method is proposed by [Chakrabarti et al., 2006]. In their approach, the smoothness property is ensured by adding a temporal loss to the original loss of static clustering. The temporal loss penalizes the deviation of current clustering result from the historical. Using the approach, they proposed an evolutionary k-means algorithm: each center at epoch i should be matched to the nearest center at i - 1 as a pair, and distances between all pairs of centers were summed as the temporal loss. As pointed out in [Chi et al., 2007], this heuristic approach could be unstable, i.e., sensitive to small perturbation on the centers. Using the same idea, [Chi et al., 2007] extended spectral clustering to evolutionary setting. Moreover, [Tang et al., 2008] extended evolutionary spectral clustering further to multi-relational clustering. However, in [Chi et al., 2007] and [Tang et al., 2008], data to be clustered at different time epochs should be identical, i.e., data of epochs are "snapshots" of the same set of objects at different time. These kind of methods have difficulties to deal with the scenario when data of different epochs are arbitrary I.I.D.

1 When doing clustering at epoch i, in off-line setting, the overall data of all epochs are available, while in on-line setting, only the data before epoch i are available.

samples from different underlying distributions. However, in many cases, we desire a good solution which is able to deal with the variation of data size and cluster number.

In this paper, we focus on the on-line setting when data of different epochs need not be identical. Starting from a density estimation view to clustering, we propose two general frameworks. In the first framework, i.e., historical data dependent (HDD), current model distribution is designed to approximate both current and historical data distributions. In the second framework, i.e., historical model dependent (HMD), current model distribution is design to approximate both current data distribution and historical model distribution. Both frameworks are based on the general exponential family mixture (EFM) model. As a result, all conventional clustering algorithms based on EFMs can be extended to evolutionary setting under the two frameworks. Experiments on both synthetic and real data sets demonstrate the validation of the two frameworks.

2 Notations and Preliminaries

X = (x1, · · · , xn) denotes observed data, which are i.i.d samples from an unknown underlying distribution F (x) (with density f(x)). About the superscript, x(i) denotes an item at time epoch i, while x[t] denotes an item at the t'th step in an iterating algorithm. Ef [·] is the expectation under distribution f .

2.1 Exponential FamilyMixture (EFM)

An exponential family is a probability distribution set FΨ, from which each density function can be expressed in the

formpΨ(x; θ) = exp {〈θ, T (x)〉 -

where θ, T(x), and Ψ(θ) are called natural parameter, natural statistic, and cumulant function, respectively.

[Banerjee et al., 2005] stated that each exponential family distribution can be uniquely expressed using Bregman divergence

where φ and bφ are functions uniquely determined by Ψ, dφ is the Bregman divergence derived from φ, and μ is the expectation parameter μ(θ) = EpΨ(x,θ) [T(x)]. Parameters μ and θ are linked by

For some widely used exponential families, the specific forms of above parameters can be found in [Banerjee et al., 2005].

A mixture model refers to a parametric distribution model with the following form:

where C is the component number, z ∈ C = {1, · · · ,C} is the component indicator variable, and Ξ = {αz, θz}Cz=1 are model parameters. When the components are taken in an exponential family FΨ, we get the general exponential family mixture (EFM) model. Typical examples of EFMs are GMM, multinomial mixture model (MMM), etc., with different definitions on Ψ or dφ.

2.2 Clustering as Density Estimation

From the view of statistical learning theory, density estimation is to find a model distribution p(x; Ξ) minimizing an expected loss (risk) (Fisher-Wald setting) [Vapnik, 2000]:

on the unknown true distribution F(x). Notice that L(Ξ)) + ∫ f(x) log f(x)dx = KL(f||p), where KL(·||·) denotes Kullback-Leibler (KL) divergence between two distributions. So density estimation is equivalent to minimizing the KL divergence between f(x) and p(x; Ξ).

A mixture model as Eq. (4) can be adopted to estimate f(x). Based on this mixture model, it's well known that L is difficult to minimize, and what will be minimized actually is a variational convex upper bound [Beal, 2003]:

where qx(·) is a distribution of z determined by x. The "≤" is derived from Jensen's inequality, with the "=" holding iff qx(·) = p(·|x; Ξ).

The well known EM procedure is used to minimize the variational bound G:

In E-step, qx(·) actually gives a solution to clustering. If no additional constraints are enforced upon qx(·), the optimal solution is

We call this case soft-clustering, e.g., GMM. If we constrain ∀z ∈ C, qx(z) ∈ {0, 1}, then the optimal solution is

We call this case hard-clustering, e.g., k-means. The superiority of soft-clustering is that in each E-step, the upper bound G is touched by the original loss L, i.e., L(Ξ[t] ) = G(q[t+1] x (·),Ξ[t] ), while in hard-clustering, this property does not hold. However, when p(x; Ξ) is an EFM model, using the Bregman divergence expression (Eq. (2)), hardclustering (Eq. (9)) is efficient to compute.

In M-step, when p(x; Ξ) is an EFM model, simply using Lagrangian method [Beal, 2003], we obtain the closed form of the solution: ∀z ∈ C,α[t+1] z = Ef [qx(z)] (10)

and

Then θ[t+1] z can be obtained by Eq. (3). In fact, using Eq. (2), we do not need θ[t+1] z in the EM iterations.

Typical examples of clustering via EFMs are GMM clustering and k-means.

3 Frameworks

Now consider the setting of on-line evolutionary clustering. At each time epoch i, new data X(i) = (x(i) 1 , . . . , x(i) ni ) arrives, and a partition on X(i) is desired. The underlying distribution is denoted as F (i) (x) (with density f (i) (x)). For each epoch, f (i) is approximated by an EFM model p(i) (x; Ξ(i) ) = ∑Ci z α(i) z p(x; θ(i) z ). The expectation parameter of the component p(x; θ(i) z ) is μ(i) z . The component numbers Ci need not be the same at different epochs.

Following [Chakrabarti et al., 2006; Chi et al., 2007], an first-order Markovian property is assumed for the evolving behavior. Therefore, we only need to consider adjacent epochs i and i+1. From now on, they will be simply denoted as epochs "1" and "2".

From the density estimation view to clustering, the loss of static clustering via an EFM is the divergence between data distribution f and the EFM model distribution p(x; Ξ). In evolutionary setting, the data distributions f (1) and f (2) are assumed close to each other. If model distributions p(1) and p(2) are their good estimates respectively, naturally, current model distribution p(2) should neither deviate much from historical data distribution f (1) nor historical model distribution p(1) . This viewpoint results in our two general frameworks for on-line evolutionary EFM: Historical Data Dependent (HDD) and Historical Model Dependent (HMD).

The general form of loss function for HDD is

Lhdd = (1 - λ) dist(f (2) , p(2) ) + λ dist(f (1) , p(2) ) (12) where the temporal loss dist(f (1) , p(2)) ensures current model distribution p(2) dose not deviate much from historical data distribution f (1) .

The general form of loss function for HMD is

Lhmd = (1 - λ) dist(f (2) , p(2)) + λ dist(p(1) , p(2) ) (13) where the temporal loss dist(p(1) , p(2) ) ensures current model distribution p(2) dose not deviate much from historical model distribution p(1) .

In both frameworks, parameter λ reflects the preference to historical data/model. The dynamic evaluation of λ will be discussed in Sec. 3.3.

3.1 Historical Data Dependent (HDD)

Since the loss of static clustering via an EFM is the KL divergence between true distribution f and the EFM model distribution p(x; Ξ), we also define the temporal loss as dist(f (1) , p(2) ) = KL(f (1) ||p(2) ), then we get the specific form of loss for HDD:

With constant item ignored, it can be easily written as

Notice that (1 - λ)f (2) (x) + λf (1) (x) induces another distribution, denoted by ˜ fλ(x). Then we have

Comparing Eq. (14) with Eq. (5), we can see that HDD is essentially to estimate the density of the deduced distribution ˜ fλ(x) using an EFM. The same EM procedure as Eq. (8, 9, 10, 11) can be used:

E-step: ∀z ∈ C, for soft-clustering

and for hard-clustering

M-step: ∀z ∈ C,α(2),[t+1]

where E ˜ fλ [·] = (1 -λ)Ef (2) [·]+λEf (1) [·]. Notice that, for i = 1, 2, Ef (i) [q[t+1] x (z)] and Ef (i) [q[t+1] x (z)T(x)] are the estimators of α(i) z and μ(i) z on f (i) , respectively. Above result means that, in each M-step, the estimator of parameters Ξ(2) on f (2) is adjusted by the same estimator on f (1) , to ensure that the estimated model distribution p(x; Ξ(2) ) approximates both f (1) and f (2) well.

3.2 Historical Model Dependent (HMD)

In the framework of Eq. (13), an intuitive choice of dist(p(1) , p(2) ) is also the KL divergence. However, KL divergence between two EFMs can not be exactly calculated, and approximate sampling methods are needed, which are time exhausted. We seek other divergence measures.

Rather than KL divergence, Earth Mover Distance (EMD) [Rubner et al., 1998] is another divergence measure between two distributions, which is frequently used to measure divergence between mixture models. EMD between two mixture models is defined as:

where d(p(x; θ(1) l ), p(x; θ(2) z )) is a predefined divergence measure between two components. In this paper, KL divergence is adopted, as KL divergence KL(p(x; θ(1) l )||p(x; θ(2) z )) between the two components from a same exponential family has a closed form

The loss function of HMD will be written as

Lhmd = (1 - λ) KL(f (2) , p(2) ) + λdEMD(p(1) , p(2) ). According to Eq. (19), minimizing Lhmd is equivalent to min Ξ(2) ,w L′hmd(Ξ(2) , w) = (1 - λ) KL(f (2) (x), p(2) (x; Ξ(2) ))

Similar to Eq. (6), L′hmd has a variational upper bound:

Alternative optimization is used to minimize G:

which is just the computation of EMD and can be efficiently solved by linear programming.

The result is identical to Eq. (15,16). The property for softclustering still holds here: in each q-step, the upper bound is touched, i.e., L′hmd (w,Ξ(2) ) = G(q[t+1] x (·),Ξ(2) ,w).

Using Lagrangian method, we can obtain the closed form of the optimal solution for this step: ∀z ∈ C

The result means that, in each Ξ-step, the estimators of expectation parameters μ(2) z on current data distribution f (2) are directly adjusted by the estimators μ(1) l of last epoch.

In fact, the evolutionary k-means in [Chakrabarti et al., 2006] is a special case of HMD with approximately computing of dEMD in w-step. The EFM used in kmeans is the mixture of spherical Gaussians with identical constant variance σ2 and prior 1C : p(x; Ξ) = 1C ∑z N(x;μz, σ2 I). Then the objective function in Eq. (20) is D = 1 2σ2 minw ∑l,z wlz‖μ(1) l - μ(2) z ‖2 , with the same constraints on w as those in Eq.(19, 20). [Chakrabarti et al., 2006] approximate D by ∑z ‖μ(2) z - μ(1) g(z) ‖, where g(z) = argminl ‖μ(2) z -μ(1) l ‖, which means, they assigned each current component to the nearest center at last epoch, then summed the distances between them as the divergence between two mixtures. Based on HMD, we can extend the approach of [Chakrabarti et al., 2006] to all the EFMs, resulting in the approximate HMD algorithm, which is different from HMD in w-step (Eq. (20)):

However, as pointed out in [Chi et al., 2007], this approach could be unstable, i.e., sensitive to small perturbation on the centers.

In both frameworks, the assumption is that epoch data are arbitrary I.I.D. samples from the corresponding epoch distribution, accordingly, the data sizes of different epochs need not be the same. Additionally, in both frameworks, cluster numbers Ci of different epochs are not assumed to be the same, consequently, both frameworks are able to deal with the variation of cluster number. Moreover, using different specific exponential families, both frameworks can produce a large family of evolutionary clustering algorithms.

3.3 Dynamic evaluation of λ

Parameter λ reflects the preference to historical data/model, which should be determined by the dependency between the adjacent data distributions. If current data distribution deviates much from the historical, then the impact of historical data/model should be suppressed. A mechanism for dynamic evaluation of λ is required. However, this problem has not been studied in previous works [Chakrabarti et al., 2006; Chi et al., 2007; Tang et al., 2008].

[Gretton et al., 2007] proposed a non-parametric test statistic to check the dependency between two distributions based on two sets of i.i.d samples. The empirical estimation of the test statistic is τ(X(1) ,X(2) ) = 1 n21 ∑n1 i,j k(x(1) i , x(1) j ) - 2 n1n2 ∑n1,n2 i,j k(x(1) i , x(2) j ) + 1 n22 ∑n2 i,j k(x(2) i , x(2) j ) , where k(·, ·) is an universal kernel. In this paper, the RBF kernel k(x, y) = exp{- ‖x-y‖2 σ2 } is used. In fact, τ measures the discrepancy between the two distributions.

where λ0 ∈ [0, 1] reflects a basic preference to historical data/model, and β reflects the sensitive to variation of the test statistic.

3.4 Comparisons between HDD and HMD

Now we give a comparison analysis to the two frameworks. HDD is efficient in computing and allows users to change the component family FΨ in the EFM model, e.g., from GMM to MMM, while HMD does not allow that. What's more, HMD needs more time to compute EMD, especially when cluster numbers are large. However, in general, if the basic assumption holds that f (1) and f (2) are close, HMD will perform better than HDD in preserving the smoothness of clustering results, which will be explained below.

When statically clustering via an EFM, estimation of p(2) on f (2) produces a solution p0 stuck in a local minimum of loss dist(f (2) , p(2)), as illustrated in Fig. 1 (a.2, b.2).

In HDD, as f (1) and f (2) are assumed close to each other, static loss dist(f (2) , p(2) ) and temporal loss dist(f (1) , p(2) ) are also close to each other (Fig.1 (a.1, a.2)). Then the overall loss Lhdd, the weighted sum of the two losses, is also close to them (Fig.1 (a.3)). Then the candidate solution p0 in static setting (Fig.1 (a.2)) will still be stuck in a nearby local minimum (p∗) in evolutionary setting (Fig.1 (a.3)).

In HMD (Fig.1 (b)), the candidate solution p0 in static setting can be heavily penalized by the large loss resulting from a large deviation from p(1) (the large loss dist(p(1) , p0)) in evolutionary setting. The penalty can drag out the solution from p0 and push it toward another minimum more close to p(1) .

4 Experiments

We demonstrate the validation of the two frameworks by experiments on three typical clustering algorithms based on EFMs, i.e., GMM, k-means, and multinomial mixture model (MMM). Evolutionary GMM is tested on a synthetic data designed by ourselves. Evolutionary k-means and evolutionary MMM are tested on a real text data set.

4.1 Data sets

The GMM data set are samples from an evolving 2D GMM model with noise. In 20 epochs, all the parameters of GMM model are slowly evolving. Data size is also varying. Five epochs and the overall data are illustrated in Fig. 2. This data set will be used to provide experiential evidences to the comparison analysis in Sec. 3.4. We will also demonstrate the necessity of dynamic evaluation of λ on this data set.

The real data set is "NSF Research Awards Abstracts" 2 , which consists of the abstracts describing NSF awards for basic research, covering 14 years from 1990 to 2003. We extract the field "NSF program" indicating the research area as the class label. A subset containing the top 10 classes covering 13 years (1990-2002) is selected as our experimental data. This subset consists of 19,728 documents with 15,412 identical words. There are 10 classes in the first 9 years, and 9 classes in the following 4 years. For evolutionary k-means, the tf-idf feature is used, while for evolutionary MMM, the word count feature is used.

2 http://kdd.ics.uci.edu/databases/nsfabs/nsfawards.html

4.2 Algorithms

Besides HDD and HMD, another two algorithms are considered: first, the static baseline, i.e. clustering via an EFM independently at each epoch, denoted as IND; second, the approximate HMD as Eq. (22), denoted as APP-HMD. In evolutionary k-means, the APP-HMD reduces to the algorithm of [Chakrabarti et al., 2006].

We cannot compare with PCQ and PCM of [Chi et al., 2007], as they cannot deal with the case when epoch data are arbitrary I.I.D. samples from epoch distributions, as pointed out in Sec. 1.

Besides the four algorithms, to illustrate the necessity of dynamic evaluating λ, we also run HDD and HMD with static λ on GMM data, which will be denoted as HDD-S and HMDS, respectively.

4.3 Criterions

The clustering quality at each epoch will be measured by Normalized Mutual Information (NMI), which is a widely used criterion for clustering. High value on NMI reflects good consistency with the true class label. The temporal smoothness of clustering results will be measured by the two types of temporal loss: KL(f (1) ||p(2) ) and dEMD(p(1) , p(2) ). They will be called data measured temporal loss (DTL) and model measured temporal loss (MTL), respectively. LowDTL andMTL reflect good smoothness of clustering results along time.

4.4 Methodology

For all algorithms, the data epochs are traversed through by N (=50) times. In one traverse, at each epoch, an identical initialization generated randomly is imposed to all algorithms. Then the criterions, e.g., NMI, DTL, MTL are calculated at each epoch. The mean and standard deviation of the criterions at each epoch are calculated across the N runs.

On GMM data, we demonstrate the necessity of dynamic evaluation of λ by another experiment: we replace the 16th epoch with the 4th epoch, which makes the epoch distribution change abruptly at the 16th and 17th epochs. We run HDD-S and HMD-S and compare the performance of them with that of HDD and HMD.

-0.2

-0.4

-0.6

-0.8

4.5 Results

The results are illustrated in Fig. (3). In general, compared to IND, HDD and HMD enhance the clustering quality at each epoch, i.e., higher mean value and much lower deviation value on NMI (Fig.3(a), 3(d), 3(g)), meanwhile, HDD and HMD better preserve the smoothness of clustering results along time, i.e., much lower mean and deviation value on DTL (Fig. 3(e), 3(b), 3(h)) and MTL (Fig. 3(f), 3(c), 3(i)). APP-HMD does not perform well, the results on NMI are even worse than those of IND on GMM data.

Results on GMM data (Fig.3(g),3(h),3(i)) provide evidences to our analysis in Sec.3.4. Due to local optima, IND gives results with large deviation (large value on deviation of NMI, mean of DTL and mean of MTL). HDD is easy to be stuck by the same local optima, resulting in the almost same result as that of IND. HMD gives the best performance.

Fig. 3(j) demonstrates the necessity of dynamic evaluation of λ: at 16th and 17th epochs, due to the abrupt change of data distribution, with static λ(= λ0), the historical data/model harms current clustering. However, dynamic evaluation of λ as Eq. (23) gives rather low value of λ at the two epochs, suppressing the impact of historical data/model. For clarity, only mean values on NMI are plot in Fig. 3(j).

5 Conclusion

We deal with the problem of evolutionary clustering where distribution of data evolves along time. We propose two general density estimation based online frameworks. They give uniform evolutionary solutions to all the conventional clustering algorithms based on EFMs.

Acknowledgments

This research was supported by National Science Foundation of China (Grant No. 60835002 and 60721003).