Distributed Clustering Based on Sampling Local Density Estimates

1 Introduction

Knowledge discovery is a process aiming at the extraction of previously unknown and implicit knowledge out of large databases, which may potentially be of added value for some given application [Fayyad et al., 1996].

Data mining, which is devoted to the automated extrac­ tion of unknown patterns from given data, is a central element among the steps of the overall knowledge discovery process; the steps include preparation of the data to be analyzed as well as evaluation and visualization of the discovered knowledge. The large variety of data mining techniques which have been developed over the past decade include: methods for pattern­ based similarity search, cluster analysis, decision-tree based classification, generalization taking the data cube or attribute­ oriented induction approach, and mining of association rules [Chen et a/., 1996]

The increasing demand to scale up to massive data sets which are inherently distributed over networks with limited bandwidth and computational resources has led to methods for parallel and distributed knowledge discovery [Kargupta eta/., 2000]. The related pattern extraction problem in dis­ tributed knowledge discovery is referred to as distributed data mining. Distributed data mining is expected to perform par­ tial analysis of data at individual sites and then to send the outcome as partial result to other sites where it is sometimes aggregated to the global result.

One of the most common approaches of business applica­ tions to perform distributed data mining is to centralize dis- Department of Electronics, Computer Science and Systems University ofBologna Via Rasi e Spinelli, 176 I-47023 Cesena FC, Italy gmoro@deis.unibo.it tributed data into a data warehouse on which to apply the usual data mining techniques. Data warehousing is a popular technology which integrates data from multiple data sources into a single repository in order to efficiently execute complex analysis queries [Moro and Sartori, 200 1]. However, despite its commercial success, this approach may be impractical or even impossible for some business settings, for instance:

• when huge amounts of data are (frequently) produced at different sites and the cost for their centralization cannot scale in terms of communication, storage and computa­ tion;
• whenever data owners cannot or do not want to release information, for instance to protect privacy or because disclosing such information may result in a competitive advantage or a considerable commercial added value.

One of the most studied data mining techniques in central­ ized environments is data clustering. The goal of this tech­ nique is to decompose or partition a data set into groups such that both intra-group similarity and inter-group dissimilarity are maximized. Despite the success of data clustering in cen­ tralized environments, only a few approaches to the problem in a distributed environment are available to date.

In this work we present KDEC, a novel approach to dis­ tributed data clustering based on sampling density estimates. In KDEC each data source transmits an estimate of the prob­ ability density function of its local data to a helper site, and then executes a density based clustering algorithm that is driven by the overall density estimate, which is built by the helper from the samples of the local densities.

The paper is organized as follows. In Section 2 we de­ scribe related work and highlight differences with respect to our approach. Section 3 and 4 present the KDEC scheme to distributed data clustering. Finally, Section 5 concludes the paper and outlines ongoing and future research work.

2 Related work

In [Johnson and Kargupta, 1999] a tree clustering approach is taken to build a global dendrogram from individual dendro­ grams that are computed at local data sites subject to a given set of requirements. In contrast to the approach presented in this paper, the distributed data sets are assumed to be hetero­ geneous, therefore every site has access only to a subset of the features of an object. The proposed solution implements a distributed version of the single-link clustering algorithm which generates clusters that are substantially different from the ones generated by density-based methods. In particular, it suffers from the so-called chaining effect, by which any of two well separated and internally homogeneous groups of objects connected only by a dense sequence of objects are re­ garded as a single cluster. [Kargupta et al., 2001] proposes a technique for distributed principal component analysis, Col­ lective PCA. It is shown that the technique satisfies efficiency and data security requirements and can be integrated with ex­ isting clustering methods in order to cluster distributed, high­ dimensional heterogeneous data. Since the dimensionality of the data is reduced prior to clustering by applying PCA, the approach is orthogonal to ours. Another related research direction deals with incremental clustering algorithms. The BIRCH [Zhang et al., 1996] and related BUBBLE method [Ganti et al., 1999], compute the most accurate clustering, given the amount of memory available, while minimizing the number of 1/0 operations. It uses a dynamic index struc­ ture of nodes that store synthetic, constant-time maintainable summaries of sets of data objects. The method is sufficiently scalable requiring O(NlogN) time and linear 1/0. However, since it uses the centroid to incrementally aggregate objects, the method exhibits a strong bias towards globular clusters. IncrementalDBSCAN [Ester et al., 1998] is a dynamic clus­ tering method supporting both insertions and deletions, which is shown to be equivalent to the well-known static DBSCAN algorithm. Since in turn DBSCAN can be shown to be equiv­ alent to a method based on density estimation when the kernel function is the square pulse and the clusters are density-based, IncrementalDBSCAN is less general than methods based on kernel density estimates. Its time complexity is O(NlogN).

3 Data Clustering

3.1 The cluster analysis problem

Cluster analysis is a a descriptive data mining task which aims at decomposing or partitioning a usually multivariate data set into groups such that the data objects in one group are similar to each other and are different as possible from those in other groups. Therefore, a clustering algorithm 51(·) is a mapping from any data set S of objects to a clustering of S, that is, a collection of pairwise disjoint subsets of S. Clus­ tering techniques inherently hinge on the notion of distance between data objects to be grouped, and all we need to know is the set of interobject distances but not the values of any of the data object variables. Several techniques for data cluster­ ing are available but must be matched by the developer to the objectives of the considered clustering task [Grabmeier and Rudolph, 2002]. In partition-based clustering, for example, the task is to partition a given data set into multiple disjoint sets of data objects such that the objects within each set are as homogeneous as possible. Homogeneity here is captured by an appropriate cluster scoring function. Another option is based on the intuition that homogeneity is expected to be high in densely populated regions of the given data set. Con­ sequently, searching for clusters may be reduced to searching for dense regions of the data space which are more likely to be populated by data objects. That leads us to the approach of density estimation based clustering.

3.2 Density estimation based clustering

In density estimation (DE) based clustering the search for densely populated regions is accomplished by estimating a so-called probability density function from which the given data set is assumed to have arisen. Many techniques for DE-based clustering are available from the vast KDD liter­ ature [Ankerst et al., 1999; Ester et al., 1996; Schikuta, 1996; Hinneburg and Keirn, 1998] and statistics [Silverman, 1986]. In both areas, the proposed clustering methods require the computation of a non-parametric estimation of the density function from the data. One important family of non­ parametric estimates is known as kernel estimators. The idea is to estimate a density function by defining the density at any data object as being proportional to a weighted sum of all objects in the data set, where the weights are defined by an appropriately chosen kernel function. In the following we introduce kernel-based density estimation [Parzen, 1962; Silverman, 1986] and our approach to density estimation based clustering.

Let us assume a set S = {x; J i = 1, .. . ,N} t;;; JRn of data points or objects. Kernel estimators originate from the in­ tuition that the higher the number of neighbouring data ob­ jects x; of some given object x E JRn, the higher the density at this object x. However, there can be many ways of cap­ turing and weighting the influence of data objects. When given the distance between one data object x and another x; as an argument, the influence of x; may be quantified by us­ ing a so called kernel function. A kernel function K (x) is a real-valued, non-negative function on lR which has finite integral over R When computing a kernel-based density es­ timation of the data set S, any element x; in S is regarded as to exert more influence on some x E JRn than elements which are farther from x than the element. Accordingly, kernel functions are often non-increasing with Jxj. Promi­ nent examples of kernel functions are the square pulse func­ tion i(sign(x+ 1)- sign(x- 1)), and the Gaussian function vk exp( -t.xZ). 1

A kernel-based density estimate <PK,h [S] ( ·): JRn --4 JR+ is defined, modulo a normalization factor, as the sum over all data objects x; inS of the distances between x; and x, scaled by a factor h, called window width, and weighted by the kernel function K:

The influence of data objects and the smoothness of the estimate is controlled by both the window width h and the shape of kernel K: h controls the smoothness of the estimate, whereas K determines the decay of the influence of a data object according to the distance. Even if the number N of data objects is very large, in practice it is not necessary to compute N distances for calculating the kernel density estimate at a given object x. In fact, the value of commonly used kernel functions is negligible for distances larger than a few h units;

1 Where Vx E lR \ {0} : sign(x) = x/lxl. and sign(O) = 0.

it may even be zero if the kernel has bounded support, as it is the case, for example, for the square pulse. Using kernel­ based density estimation, it is straightforward to decompose the clustering problem into three phases as follows.

1. Choose a window width h and a kernel function K.
2. Compute the kernel-based density estimate {pK,h [S] (.X) from the given data set.
3. Detect regions of the data space where the value of the estimate is high and group all data objects of space re­ gions into corresponding clusters.

In the literature, many different definitions of cluster have been proposed formalizing the clusters referred to in step 3 above. A density-based cluster [Ester et a!., 1996] collects all data objects included in a region where density exceeds a threshold. Center-defined clusters [Hinneburg and Keirn, 1998] are based on the idea that every local maximum of {p corresponds to a cluster including all data objects which can be connected to the maximum by a continuous, uphill path in the graph of (p. Finally, an arbitrary-shape cluster [Hinneb­ urg and Keirn, 1998] is the union of center-defined clusters having their maxima connected by a continuous path whose density exceeds a threshold.

Algorithm 1 (DE-cluster) implements the computation of center-defined clusters by a climbing procedure driven by the density estimate. The main procedure is DECluster, taking as inputs an instanceS of the class of data objects, the kernel function K, the window width h, and returning a clustering represented by C, which stores a mapping from each .X; to the unique integer label of x; 's cluster. It is assumed that S is an instance of a class which provides the following meth­ ods: get(i) to access object .X; given index i, NQ (k,x) and Radius(k,x) to retrieve, given x E JR, the indexes and maxi­ mum distance of X's k nearest neighbours. Uphill computes the steepest direction on the graph of the estimated density as the versor of its gradient, computed by function DEGradient (cf. [Hinneburg and Keirn, 1998]). Uphill then moves in that direction a fraction o < 1 of the distance S.Radius(k,x) of the kth nearest neighbour ofx, and finally returns the index of the nearest neighbour in S of the reached position. Every nested call to FixedPoint marks the current object x; as visited and calls Uphill to get the index of the next data object x 1 . If such object has already been visited, the proximity of a local maximum has been reached and j is taken as new cluster la­ bel. Otherwise x 1 is inductively assumed to lie at the bottom end of a path leading to the proximity of a local maximum, and to be already labeled accordingly. (If x 1 is not marked as clustered, a recursive call ensures that the assumption holds.)

The complexity of the DE-cluster algorithm is that of call­ ing N = S.coimt times FixedPoint. At the beginning of every iteration in DECJuster, the sets of clustered and visited ob­ jects are equal. FixedPoint is never called with a clustered object as argument, and visits unclustered objects at most once. Therefore, even if the number of visited data objects in one call of FixedPoint is bounded only by N, the number of visited data objects in all calls is still only N. For each visited object a single k-nearest neighbour query suffices to compute the gradient and the next uphill object. The methods NQ ( ·, ·) and Radius ( ·, ·) can be efficiently implemented by Algorithm 1 DE-cluster: Clustering based on density estima­ tion

functDEGradient(x,S ,K() ,h) = foreach i E S.NQ (k,x) do g := g+ (x; -.X) *K(d(x,x;)/h) od g. funct Uphill ( i, S, K () , h) = .X:= S.get(i) v := DEGradient(x,S ,K ,h) I IIDEGradient(x,S,K,h)ii S.NQ (1,x+o *S.Radius(k,x) * v). proc FixedPoint(i ,S ,K() ,h,C) C.setVisited(i) j := Uphill(i,S,K,h) !f C.clustered (j) then C.setLabeJ (i, C.getLabeJ (j)) else !f C. visited(j) then C.setLabeJ (i , j) else FixedPoint(j,S,K,h,C) C.setLabeJ (i , C.getLabel (j)) ft C.setCJustered (i). proc DECluster(S,K() ,h ,k ,C) fori:= 1 to S.count do !f-, C.clustered (i) then FixedPoint ( i, S ,K , h, C) ft od.

equipping the class of S with a spatial access method like the KD-, or MVP-, or M-tree. Therefore, the time complexity of DECluster is O(Nq(N) ), where q(N) is the cost of a k near­ est neighbour query in any such access method. Note that in many practical cases, q(N) is very close to logN.

4 Distributed Data Clustering

The body of work on applications of data clustering in dis­ tributed environments, the problem of so called distributed data clustering (DDC), is comparatively small. In this sec­ tion we adopt the kernel density estimation based clustering approach presented above for the distributed case assuming homogeneous data, which means that a data object cannot be split across two sites.

4.1 The DDC Problem

We define the problem of homogeneous distributed data clus­ tering for a clustering algorithm Jt as follows. Let S = {.X; I i = 1, ... ,N} s;;; JRn be a data set of objects. Let L1,j = I, ... ,M, be a finite set of sites. Each site L 1 stores one data set DJ, and it will be assumed that S = UX= 1 DJ. The DDC problem is to find a site clustering c 1 residing in the data space of L1, for j = 1, ... ,M, such that

• (i). C1 = {CnD1 : C E Jt(S)} (correctness requirement)
• (ii). Time and communications costs are minimized (efficiency requirement)
• (iii). At the end of the computation, the size of the subset of S which has been transferred out of the data space of any site L J is minimized (privacy requirement).

The traditional solution to the homogeneous distributed data clustering problem is to simply collect all the distributed data sets Dj into one centralized repository where their union S is computed, and the clustering C of the union S is com­ puted and transmitted to the sites. Such approach, however, does not satisfy our problem's requirements both in terms of privacy and efficiency. We therefore propose a different ap­ proach yielding a kernel density estimation based clustering scheme, called KDEC.

4.2 The KDEC Scheme for DDC

The key idea of the KDEC scheme is based on the follow­ ing observation: Although the density estimate computed on each local data set gives information on the distribution of the objects in the data set, it conceals the objects themselves. Moreover, the local density estimate can be coded to provide a more compact representation of the data set for the purpose of transmission. In the sequel, we tacitly assume that all sites L 1 agree on using a global kernel function K and a global window width h. We will therefore omit K and h from our notation, and write <P[S](x) for <PK,h[S](x). Density estimates in the form of Equation (1) are additive, i.e. the global den­ sity estimate <P[S](x) can be decomposed into the sum of the site density estimates, one estimate for every data set Dr

Thus, the local density estimates can be transmitted to and summed up at a distinguished helper site yielding the global estimate which can be returned to all sites. Each site Lj may then apply, in its local data space, the hill-climbing technique of Algorithm I (DE-cluster) to assign clusters to the local data objects. There is nevertheless a weakness in such a plan: the definition of a density estimate explicitly refers to all the data objects x;. Hence, knowing how to manipulate the estimate entails knowing the data objects, which contradicts the pri­ vacy requirement. However, only an intensional, algebraic definition of the estimate includes knowledge of the data ob­ jects. Multidimensional sampling theory provides the basis for an alternative extensional representation of the estimate which makes no explicit reference to the data objects.

The theoretical idea of sampling is to represent a function f by a sampling series, that is, a summation of suitable ex­ pansion functions weighted by the values of f at a discrete subset of its domain [Higgins, 1996]. In the following, let the i th coordinate of x E !Rn be denoted by X(i), and let Diag[u], u = [u(l l, ... , u(n)jT E !Rn, denote the n x n diagonal matrix having diagonal u, i.e., defined by Diag[u];J = 0 if i =f. j, Diag[u];; = uUl. Further let 1: = [·(1), ... , -.(n)jT E !Rn a vector of sampling periods. The sampled form of <!>[S](x) at intervals 1: is the sequence { <i>z[S]}, z E zn, defined by

where · is the inner product between vectors. Therefore, <Pz[S] is the sequence of the values of <!>[S] (x) at all the real, n­ dimensional vectors whose i th coordinates are spaced by a multiple of the i th sampling period -.(i), i = I, ... , n. The sam­ pled forms of the local density estimates are defined in a sim­ ilar way by <Pz[D1] = <P[D1](Diag[z]·1:), for j = 1, ... ,M. It is immediate to see by (2) that addivity holds for the sampled forms:

Therefore, after receiving the sampled forms <Pz[Dj] of the M density estimates, the helper site can compute by (4) the sampled form of the overall estimate and transmit it to the sites LJ. Sites LJ can then cluster local data with respect to the overall density estimate, using the gradient of the sampling series

where Sinc(x) = ITI'= l sinc(x(il), and

as needed in the hill-climbing function.

We briefly discuss the extent to which the series (5) can be used to represent <P [ S] (x). It is well known that, under mild conditions, sampling a function g(x) is an invertible transformation if, for every coordinate i = 1, ... , n, there is a frequency f(i) such that the Fourier transform of g differs from zero only in the interval [-f(il,t(il], and the samples are computed with a period not greater than -.(i) = 1 I (2/(i)) (cf. [Higgins, 1996]). Under these assumptions, the value of the sampling series computed at x equals g(x). Unfortunately, most popular kernel functions (hence summations of kernel functions) do not satisfy these hypotheses since the support of their Fourier transform is unbounded. Consequently sam­ pling density estimates yields an information loss. However, it can be shown that the Fourier transform of a kernel density estimate is negligible everywhere except for l!(i) I not greater than 1 I h. Therefore, the global density estimate can be re­ constructed from its samples by (5) introducing only a small error if 't(i) < hl2.

It is worth noting that the infinite series (5) need not be ap­ proximated, if it has finitely many nonzero terms. The latter case holds if the used kernel function has a bounded support, since the density estimate will also have bounded support. If, however, the kernel function has unbounded support, like the Gaussian kernel, then the density estimate can be approxi­ mated by regarding its value to be zero everywhere except inside an appropriately chosen bounded region.

According to this approach, we propose the following al­ gorithmic KDEC scheme for computing the kernel density es­ timation based clusters for local data spaces at M distributed data sites L1 (see Algorithm 2). Every local site runs the pro­ cedure SiteDECluster whereas the helper site runs Helper. SiteDECluster is passed a reference H to the helper and the local data set D, and returns a clustering in a class instance C. Helper is passed a list of references L to the local sites. The procedure SiteNegotiate carries out a negotiation with the other local sites through HelperNegotiate at the helper site to determine the sampling periods 1:, the bounding cor­ ners of the sampling rectangle z1, z2 E zn, the kernel K, and Algorithm 2 KDEC: distributed clustering based on density estimation

functSampJe(D,'t,z1,hK(),h) = s.setComers(z1 ,zz) foreach z E RectangJe(s.getComers()) do s.Add(z,DE(Diag[Z] ·'t,D,K,h)) od s. funct Uphill (x, 't, GlobaJSam) = if II SeriesGradient(x, 't, GlobalSam) II > E then x + 8 * SeriesGradient(x, 't, GlobaJSam) else x ft. funct FixedPoint(x,D, GlobalSam,C) ji := Uphill(x,'t, GlobalSam) q := D.RangeQ (&,ji) if :liE q :C. clustered (i) then Id := C.getLabel (i) elsifji = x then Id := D.Add(x) else Id : = FixedPoint (ji, D, Global Sam, C) !! foreach i E D.RangeQ (&,x) do C.setLabel(i,Id); C.setClustered(i) od Id. proc SiteDECJuster(D,H,C) = SiteNegotiate(H, 't,z1 ,z2 ,K,h) Send(H, SampJe(D,'t,z1 ,hK,h)) GlobalSam := Receive(H) fori:= 1 to D.count do if · C. clustered (i) then FixedPoint(xi,D, GlobalSam,C)!! od. proc HeJper(L) = HeJperNegotiate (L) for j := 1 to L.count do Sam:= Receive(L.get(j)) foreach z E Rectangle (Sam .getComers ()) do GlobalSam.Sum(z, Sam.get(z)) od od for j := 1 to L.count do Send(GlobaJSam,L.get(j)) od.

the window width h. The negotiation procedures contain ap­ propriate handshaking primitives to ensure that all sites par­ ticipate and exit negotiations only if an agreement has been reached. Each local site computes the sampled form of the estimate of D, by calling function Sample, and sends it to the helper. (Function DE computes the density estimate and is omitted for brevity.) The helper receives the sampled es­ timates, sums them by sampling indexes into a global sam­ ple, and retuins it to all sites. Procedures Send and Receive implement appropriate blocking and handshaking to ensure the transmission takes place. Each local site uses the global sample in functions FixedPoint and Uphill to compute the values of the gradient of the global density estimate. Func­ tion SeriesGradient can be easily derived from (5). Local sites perform a DE-cluster algorithm to compute the corre­ sponding local data clusters. The details of the hill-climbing strategy are however different from Algorithm 1 because the sites are allowed access to local data objects only. Uphill advances a fraction 8 of the gradient in its direction, if the gradi­ ent's norm exceeds a threshold E. If a &-neighbourhood of the object returmed by Uphill contains an already clustered data object, the current cluster label Jd is set from that object's la­ bel. Otherwise, checking whether Uphill returned the same space object signals to FixedPoint that the proximity of the local maximum has been reached. The maximum is marked by adding the current space objectx as a dummy object to the local data set; this ensures that subsequent paths converging to the same local maximum will use the same cluster label as the current path. Method D.Add(·) returns the identifier of the added object, which is used as current cluster label. If neither case holds, the label Jd is obtained by a recursive call. Finally, all objects in a small neighbourhood of the current object are labeled by Id. Note that adding dummy objects has effect only on the range queries, and does not modify the density estimate.

4.3 Complexity of the KDEC scheme

In terms of the complexity in computation and communica­ tion one crucial point of the KDEC scheme is how many sam­ ples have to be computed and transferred among the sites. In most cases, to obtain good density estimates, h must not be less than a small multiple of the smallest object distance. As 't(i) ~ h/2, the number of samples should rarely exceed the number of objects, if only space regions where the density estimate is not negligible are sampled. Since the size of a sample is usually much smaller than the size of a data ob­ ject, the overall communication costs of our DDC approach will be in most cases significantly lower than in a centralized approach. Of course, the precise number of samples depends on the bounding region that is being sampled by every site. In Algorithm 2 the site L 1 determines autonomosly the rectangle that contains the computed samples.

The computational costs of the KDEC scheme in terms of used CPU cycles and 1/0 do not exceed the one in the cen­ tralized approach where clustering is performed on data col­ lected in a single repository. The computational complexity is linear in the number of samples. The precise cost of compu­ tation of any KDEC-based DDC algorithm as an instance of the proposed scheme largely depends also on the used kernel function and local clustering algorithm. The DE-cluster algo­ rithm we developed for the KDEC scheme in Section 3.2 is of complexity O(Nq(N)), where q(N) is the cost of a nearest neighbour query (which in practical cases is close to logN). Algorithm 2 implements a slightly different approach in the hill-climbing function than Algorithm 1, since the function does not use data objects to direct the uphill path. How­ ever, preliminary results of experiments conducted on a pro­ totype implementation show good scalability of the approach, in terms of number of executed range queries.

5 Conclusion

Due to the explosion in the number of autonomous data sources there is a growing need for effective approaches to distributed knowledge discovery and data mining. In this pa­ per we have presented KDEC, a novel scheme for distributed data clustering which computes the density estimation, to perform the clustering, from sampled forms of local densities at each data source site.

The approach exploits statistical density estimation and in­ formation theoretic sampling to minimize communications between sites. Moreover, the privacy of data is preserved to a large extent by never transmitting data values but ker­ nel based density estimation samples outside the site of ori­ gin. The approach does not require CPU and 110 costs sig­ nificantly higher than a similar centralized approach and its communication costs may be lower. Ongoing research fo­ cuses in particular on implementations of a multiagent system for KDEC-based distributed data clustering in a peer-to-peer network, and investigation on methods to mitigate the risk of security and privacy violations in distributed data mining en­ vironments.