Clustering For Point Pattern Data

TL;DR

This paper introduces two novel clustering methods for point pattern data, addressing a gap in existing algorithms, and demonstrates their effectiveness on simulated and real datasets.

Contribution

It proposes a non-parametric distance-based method and a model-based approach using finite set theory for clustering point patterns.

Findings

Methods perform well on simulated data

Methods are effective on real-world data

New approaches outperform existing techniques

Abstract

Clustering is one of the most common unsupervised learning tasks in machine learning and data mining. Clustering algorithms have been used in a plethora of applications across several scientific fields. However, there has been limited research in the clustering of point patterns - sets or multi-sets of unordered elements - that are found in numerous applications and data sources. In this paper, we propose two approaches for clustering point patterns. The first is a non-parametric method based on novel distances for sets. The second is a model-based approach, formulated via random finite set theory, and solved by the Expectation-Maximization algorithm. Numerical experiments show that the proposed methods perform well on both simulated and real data.

Figures

Equations

$$d_{\mathtt{H}}(X,Y)\triangleq \max \left\{\underset{x\in X}{\max }\underset{y\in Y}{\min }d(x,y),\underset{y\in Y}{\max }\underset{x\in X}{\min }d(x,y)\right\}.$$

$$d_{\mathtt{W}}^{(p)}(X,Y)\triangleq \underset{C}{\min }{\left(\sum _{i=1}^m\sum _{j=1}^n{c}_{i,j}d{\left(x_i,y_j\right)}^p\right)}^{1/p},$$

$$\begin{array}[]{ccc}c_{i,j}\geq 0&\mbox{ for }1\leq i\leq m,&1\leq j\leq n\\ \sum_{j=1}^{n}c_{i,j}=\frac{1}{m}&\mbox{ for }1\leq i\leq m,\\ \sum_{i=1}^{m}c_{i,j}=\frac{1}{n}&\mbox{for }1\leq j\leq n.\end{array}$$

$$d_{\mathtt{O}}^{(p,c)}(X,Y)\triangleq$$

$${\left(\frac{1}{n}\left(\underset{\pi \in {\Pi }_n}{\min }\sum _{i=1}^m{\left(\min \left(c,d\left(x_i,y_{\pi (i)}\right)\right)\right)}^p+c^p\left(n-m\right)\right)\right)}^{1/p},$$

$$\mu (\mathcal{T})=\sum _{m=0}^{\infty }\frac{1}{m!U^m}\int 1_{\mathcal{T}}\left(\{x_1,...,x_m\}\right)d\left(x_1,...,x_m\right)$$

$$p(X)=p_c(m)m!U^m{f}_m\left(x_1,...,x_m\right),$$

$$p(X)=p_c(|X|)|X|![U{p}_f{]}^X$$

$$p(X)={\lambda }^{|X|}{e}^{-\lambda }[U{p}_f{]}^X$$

$$p(X\mid \Theta )=\sum _{k=1}^{N_{\mathrm{comp}}}{w}_{k}p(X\mid {\mathcal{C}}_k,{\mathcal{F}}_k)$$

$$p(X\mid {\mathcal{C}}_k,{\mathcal{F}}_k)=p_c(|X|\mid {\mathcal{C}}_k)|X|!U^{|X|}\prod _{x\in X}{p}_f(x\mid {\mathcal{F}}_k)$$

$$p(k\mid X_n,{\Theta }^{(i\text{mbox}-1)})=\frac{w_k^{(i\text{mbox}-1)}p(X_n\mid {\mathcal{C}}_k^{(i\text{mbox}-1)},{\mathcal{F}}_k^{(i\text{mbox}-1)})}{{\sum }_{\ell =1}^{N_{\mathrm{comp}}}{w}_{\ell }^{(i\text{mbox}-1)}p(X_n\mid {\mathcal{C}}_{\ell }^{(i\text{mbox}-1)},{\mathcal{F}}_{\ell }^{(i\text{mbox}-1)})}$$

$$w_k^{(i)}=\frac{1}{N_{\mathrm{data}}}\sum _{n=1}^{N_{\mathrm{data}}}p(k\mid X_n,{\Theta }^{(i\text{mbox}-1)})$$

$$q_{k,m}^{(i)}=\frac{{\sum }_{n=1}^{N_{\mathrm{data}}}{\delta }_{(m\text{mbox}-|X_n|)}p(k\mid X_n,{\Theta }^{(i\text{mbox}-1)})}{{\sum }_{\ell =0}^{N_{\mathrm{card}}}{\sum }_{n=1}^{N_{\mathrm{data}}}{\delta }_{(\ell \text{mbox}-|X_n|)}p(k\mid X_n,{\Theta }^{(i\text{mbox}-1)})}$$

$${\mathbf{\mu }}_k^{(i)}=$$

$${\mathbf{\Sigma }}_k^{(i)}=$$

$${\mathcal{C}}_k$$

$${\mathcal{F}}_k$$

$$Q\left(\Theta ,{\Theta }^{(i-1)}\right)$$

$$=\sum _{k=1}^{N_{\mathrm{comp}}}\sum _{n=1}^{N_{\mathrm{data}}}\log \left(w_{k}p(X_n\mid {\mathcal{C}}_k,{\mathcal{F}}_k)\right)p(k\mid X_n,{\Theta }^{(i\text{mbox}-1)})$$

$${\hat{k}}_n=\underset{k\in \mathbb{K}}{\text{mbox}argmax}p(k\mid X_n,\Theta ),$$

$$p(k\mid X_n,\Theta )=\frac{w_{k}p(X_n\mid {\mathcal{C}}_k,{\mathcal{F}}_k)}{{\sum }_{\ell =1}^{N_{\mathrm{comp}}}{w}_{\ell }p(X_n\mid {\mathcal{C}}_{\ell },{\mathcal{F}}_{\ell })}.$$

Full text

††footnotetext: Dinh Phung gratefully acknowledges support from the Air Force Office of Scientific Research under award number FA2386-16-1-4138.

I Introduction

Clustering is a data analysis task that groups similar data items together [1] and can be viewed as an unsupervised classification problem since the class (or cluster) labels are not given [2, 3]. Clustering is a fundamental problem in machine learning with a long history dating back to the 1930s in psychology [4]. Today, clustering is widely used in a host of application areas including genetics [5], medical imaging [6], market research [7], social network analysis [8], and mobile robotics [9]. Excellent surveys can be found in [2, 10].

Clustering algorithms can broadly be categorized as hard or soft. In hard clustering, each datum could only belong to one cluster [2, 11], i.e. a hard clustering algorithm outputs a partition of the dataset111More specifically, the hard clustering could be dichotomized as either hierarchical or partitional. The hierarchical clustering has the output as a nested tree of partitions, whereas the output of partitional clustering is only one partition [1, 2].. K-means is a typical example of hard clustering. Soft clustering, on the other hand, allows each datum to belong to more than one cluster with certain degrees of membership. The Gaussian mixture model [3] is an example of soft clustering wherein the degree of membership of a data point to a cluster is given by its mixing probability. Hard clustering can be obtained from soft clustering results by simply assigning the cluster with the highest membership degree to each data point [2].

In many applications, each datum is a point pattern, i.e. set or multi-set of unordered points (or elements). For example, in natural language processing and information retrieval, the ‘bag-of-words’ representation treats each document as a collection or set of words [12, 13]. In image and scene categorization, the ‘bag-of-visual-words’ representation – the analogue of the ‘bag-of-words’ – treats each image as a set of its key patches [14, 15]. In data analysis for the retail industry as well as web management systems, transaction records such as market-basket data [16, 17, 18] and web log data [19] are sets of transaction items. Other examples of point pattern data could be found in drug discovery [20], and protein binding site prediction [21]. In multiple instance learning [22, 23], the ‘bags’ are indeed point patterns. Point patterns are also abundant in nature, such as the coordinates of trees in a forest, stars in a galaxy, etc. [24, 25, 26].

While point pattern data are abundant, the clustering problem for point patterns has received very limited attention. To the best of our knowledge, there are two clustering algorithms for point patterns: the Bag-level Multi-instance Clustering (BAMIC) algorithm [27]; and the Maximum Margin Multiple Instance Clustering (M3IC) algorithm [28]. BAMIC adapts the k-medoids algorithm for the clustering of point pattern data (or multiple instance data) by using the Hausdorff metric as a measure of dissimilarity between two point patterns [27]. M3IC, on the other hand, poses the point pattern clustering problem as a non-convex optimization problem which is then relaxed and solved via a combination of the Constrained Concave-Convex Procedure and Cutting Plane methods [28].

In this paper, we propose a non-parametric approach and a model-based approach to the clustering problem for point pattern data:

• •

Our non-parametric approach uses Affinity Propagation (AP) [29] with a novel measure of dissimilarity, known as the Optimal Sub-Pattern Assignment (OSPA) metric. This metric alleviates the insensitivity of the Hausdorff metric (used by BAMIC) to cardinality differences. Moreover, AP is known to find clusters faster and with much lower error compared to other methods such as k-medoids (used by BAMIC) [29].

• •

In our model-based approach, point patterns are modeled as random finite sets. Moreover, for a class of models known as independently and identically distributed (iid) cluster random finite sets, we develop an Expectation-Maximization (EM) technique [30] to learn the model parameters. To the best of our knowledge, this is the first model-based framework for bag-level clustering of multiple instance data.

II Non-parametric clustering for

point pattern data

In AP, the notion of similarity/dissimilarity between data points plays an important role. There are various measures of similarity/dissimilarity, for example distances between observations, joint likelihoods of observations, or even manually setting the similarity/dissimilarity for each observation pair [29]. In this paper, we use distances for sets as a measures of dissimilarity. In the following section, we describe several distances for point patterns.

II-A Set distances

Let $X=\left\{x_{1},...,x_{m}\right\}$ and $Y=\left\{y_{1},...,y_{n}\right\}$ denote two finite subsets of a metric space $(\mathcal{S},d)$. Note that when $\mathcal{S}$ is a subset of $\mathbb{R}^{n}$, $d$ is usually the Euclidean distance, i.e., $d(x,y)=||x-y||$.

The Hausdorff distance is defined by

[TABLE]

and $d_{\mathtt{H}}(X,Y)=\infty$, when $X=\emptyset$ and $Y\neq\emptyset$.

While the Hausdorff distance is a rigorous measure of dissimilarities between point patterns, it is relatively insensitive to differences in cardinalities [31, 32]. Consequently, clustering based on this distance has the undesirable tendency to group together point patterns with large differences in cardinality.

The Wasserstein distance of order $p\geq 1$ is defined by

[TABLE]

where each $C=\left(c_{i,j}\right)$ is an $m\times n$ transportation matrix, i.e. $C$ satisfies [31, 32]:

[TABLE]

Note that when $X=\emptyset$ and $Y\neq\emptyset$, $d_{\mathtt{W}}^{(p)}(X,Y)=\infty$ by convention.

The Wasserstein distance is more sensitive to differences in cardinalities than the Hausdorff distance and also has a physically intuitive interpretation when the point patterns have the same cardinality. However, it does not have a physically consistent interpretation when the point patterns have different cardinalities, see [31] for further details.

The OSPA distance of order $p$ with cut-off $c$ is defined by

[TABLE]

if $m\leq n,$ and $d_{\mathtt{O}}^{(p,c)}(X,Y)\triangleq d_{\mathtt{O}}^{(p,c)}(Y,X)$ if $m>n$, where $1\leq p<\infty$, $c>0$, and $\Pi_{n}$ is the set of permutations on $\left\{1,2,...,n\right\}$. The two adjustable parameters $p$, and $c$, are interpreted as outlier the sensitivity and cardinality penalty, respectively. The OSPA metric allows for a physically intuitive interpretation even if the cardinalities of two sets are not the same, see [31] for further details.

II-B AP clustering with set distances

The AP algorithm first considers all data points as potential exemplars, i.e., centroids of clusters. Progressively better sets of exemplars and corresponding clusters are then determined by passing “responsibility” and “availability” messages based on similarities/dissimilarities between data points. Further details on the AP algorithm can be found in [29]. The dissimilarities measures considered in this work are the Hausdorff, Wasserstein, OSPA distances. In the following section, we will evaluate AP clustering performance with various set distances.

II-C Numerical experiments

We evaluate the proposed clustering algorithm with both simulated and real data222For AP clustering algorithm, we use the Frey Lab’s code, http://www.psi.toronto.edu/index.php?q=affinity%20propagation. Retrieved August 4, 2015. For computing Hausdorff distance, we use Zachary Danziger’s code, http://www.mathworks.com/matlabcentral/fileexchange/26738. Retrieved February 1, 2016. For computing OSPA distance, we use Ba-Ngu Vo’s code, http://ba-ngu.vo-au.com/vo/ospa_dist.zip. Retrieved February 1, 2016.. The performance are measured by Rand index [33]. From a performance perspective, we consider AP clustering with the Hausdorff metric as version of BAMIC since the only difference is that latter uses k-medoids instead of AP.

In the following experiments, we report the performance of AP clustering with different set distances, namely Hausdorff, Wasserstein and OSPA. We use $p=2$ for the Wasserstein and OSPA metrics, and additionally $c=20$ or $c=60$ for OSPA.

II-C1 AP clustering with simulated data

The simulated point pattern data are generated from Poisson RFS models with 2-D Gaussian feature distributions (see section III-A), with each Poisson RFS representing a cluster. There are 3 clusters, each with 100 point patterns, in the simulated data (Fig. 2).

The test is run 10 times with 10 different (simulated) datasets. The averaged results are shown in Fig. 1, while individual results for certain datasets are show in Fig. 2. Observe that the OSPA’s performance could be improved if one choose a suitable cut-off $c$ (see section II-B).

II-C2 AP clustering with real data

This experiment involves clustering images from the classes “T14_brick1” and “T15_brick2” of the Texture images dataset [34]. Fig. 3 visualizes some example images from these classes.

Features are extracted from each image by applying the SIFT algorithm333Using the VLFeat library [35]., followed by Principal Component Analysis (PCA) to convert the 128-D SIFT features into 2-D features. Thus each image is compressed into a point pattern of 2-D features. Fig. 4 plots the 2-D features of the images in the Texture images dataset. The clustering algorithm is used to group together point patterns extracted from images of the same class.

The AP clustering results for various set distances are shown in Fig. 5. Observe that for this dataset, the OSPA distance achieves a better performance than the Hausdorff and Wasserstein distances.

III Model-based clustering for

point pattern data

In this section, we present generative models for clusters of point pattern data. In particular, the point patterns are assumed to be distributed according a finite mixture of iid cluster random finite set distributions. We then derive an Expectation-Maximization (EM) algorithm to learn the parameters of these models, which are then used to cluster the data.

III-A Random Finite Set

Point patterns can be modeled as random finite sets (RFSs), or simple finite point processes. The likelihood of a point pattern of discrete features is straightforward since this is simply the product of the cardinality distribution and the joint probability of the features given the cardinality. The difficulties arise in continuous feature spaces. In this work, we only consider continuous feature spaces.

Let $\mathcal{F}(\mathcal{X})$ denote the space of finite subsets of a space $\mathcal{X}$. A random finite set (RFS) $X$ of $\mathcal{X}$ is a random variable taking values in $\mathcal{F}(\mathcal{X})$. An RFS $X$ can be completely specified by a discrete (or categorical) distribution that characterizes the cardinality $|X|$, and a family of symmetric joint distributions that characterizes the distribution of the points (or features) of $X$, conditional on the cardinality.

Analogous to random vectors, the probability density of an RFS (if it exists) is essential in the modeling of point pattern data. The probability density $p:\mathcal{F}(\mathcal{X})\rightarrow[0,\infty)$ of an RFS is the Radon-Nikodym derivative of its probability distribution relative to the dominating measure $\mu$, defined for each (measurable) $\mathcal{T}$ $\subseteq$ $\mathcal{F}(\mathcal{X})$, by [24, 26, 36, 37]:

[TABLE]

where $U$ is the unit of hyper-volume in $\mathcal{X}$, and $\mathbf{1}_{\mathcal{T}}(\cdot)$ is the indicator function for $\mathcal{T}$. The measure $\mu$ is the unnormalized distribution of a Poisson point process with unit intensity $u=1/U$ when $\mathcal{X}$ is bounded. Note that $\mu$ is unitless and consequently the probability density $p$ is also unitless.

In general the probability density of an RFS, with respect to $\mu$, evaluated at $X=\{x_{1},...,x_{m}\}$ can be written as [38, p. 27] ((Eqs. (1.5), (1.6), and (1.7)), [26]:

[TABLE]

where $p_{c}(m)=\mbox{Pr}(|X|=m)$ is the cardinality distribution, and $f_{m}\left(x_{1},...,x_{m}\right)$ is a symmetric joint probability density of the points $x_{1},...,x_{m}$ given the cardinality.

Imposing the independence assumption among the features on the model in (5) reduces to the iid-cluster RFS model [25]

[TABLE]

where $p_{f}$ is a probability density on $\mathcal{X}$, referred to as the* feature density, *and $h^{X}\triangleq\prod_{x\in X}h(x)$, with $h^{\emptyset}=1$ by convention, is the finite-set exponential notation.

When $p_{c}$ is a Poisson distribution we have the celebrated Poisson point process (aka, Poisson RFS)

[TABLE]

where $\lambda$ is the mean cardinality. The Poisson model is completely determined by the intensity function $u=\lambda p_{f}$ [36, 37]. Note that the Poisson cardinality distribution is described by a single non-negative number $\lambda$, hence there is only one degree of freedom in the choice of cardinality distribution for the Poisson model.

III-B Mixture of iid-cluster RFSs

A mixture of iid-cluster RFSs is a probability density of the form

[TABLE]

where $k\in\left\{1,...N_{\mathrm{comp}}\right\}$ is the component label, $w_{k}=\mbox{Pr}(C=k)$ is the component weight (the probability of an observation belonging to the $k^{th}$ component), and

[TABLE]

is iid-cluster density of the $k^{th}$ component with cardinality distribution parameter $\mathscr{C}_{k}$, and feature distribution parameter $\mathscr{F}_{k}$.

Note that $\Theta=\left\{\left(w_{k},\mathscr{C}_{k},\mathscr{F}_{k}\right):k=1,...,N_{\mathrm{comp}}\right\}$ is the complete collection of parameters of the iid-cluster mixture model. The probability density (8) is a likelihood function for iid point patterns arising from $N_{\mathrm{comp}}$ clusters.

III-C EM clustering using mixture of iid-cluster RFSs

Given a dataset $\mathcal{D}=\left(X_{1},...,X_{N_{\mathrm{data}}}\right)$ where each datum is a point pattern $X_{n}=\left\{x_{n,1},...,x_{n,|X_{n}|}\right\}$. Assume that $\mathcal{D}$ is generated from an underlying iid-cluster RFS mixture model with $N_{\mathrm{comp}}$ components (8) where each component represents a data cluster. The EM clustering for the given data includes two main steps:

Model learning: Estimate the parameters of the underlying model using EM method. 2. 2.

Cluster assignment: Assign observations to clusters using maximum a posteriori (MAP) estimation.

**Model learning step: **In this paper, we present EM learning for the iid-cluster mixture model (8) with categorical cardinality distribution and Gaussian feature distribution, i.e., the parameters of the $k^{th}$ component are

[TABLE]

where $N_{\mathrm{card}}$ is the maximum cardinality of point patterns, $\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k}$ are the means and covariances of the Gaussian distributions. The EM algorithm for learning the model parameters $\Theta=\left\{\left(w_{k},\mathscr{C}_{k},\mathscr{F}_{k}\right):k=1,...,N_{\mathrm{comp}}\right\}$ proceeds as shown in Algorithm 1.

In each iteration of Algorithm 1, the parameters are estimated by $\Theta^{(i)}=\mbox{argmax}_{\Theta}Q\left(\Theta,\Theta^{(i-1)}\right)$ [39], where

[TABLE]

and $p(X_{n}\mid\mathscr{C}_{k},\mathscr{F}_{k})$, $p(k\mid X_{n},\Theta^{(i\mbox{-}1)})$ are given by (9), (10) respectively.

Cluster assignment step:* *The cluster label $k_{n}$ of an observation $X_{n}$ can be estimated using MAP estimation:

[TABLE]

where $\mathbb{K}=\left\{1,...,N_{\mathrm{comp}}\right\}$, parameters $\Theta$ are learned by Algorithm 1, and $p(k\mid X_{n},\Theta)$ is the posterior probability:

[TABLE]

III-D Numerical experiments

In this section, we evaluate the proposed EM clustering with both simulated and real data. For comparison with AP clustering, we use the same datasets from section II-C. Note that our experiments assume a mixture of Poisson RFSs model.

III-D1 EM clustering with simulated data

In this experiment we apply EM clustering on the same (simulated) datasets from section II-C1. The performance is shown in Fig. 6. EM clustering performs very well for these datasets with the average Rand index of 0.95.

To further investigate performance of the proposed algorithm, we illustrate in Fig. 7 the distributions learned by the EM algorithm for the 3 datasets shown in Fig. 2. Observe that the model-based method has better performance than the non-parametric method when there is little overlap between the feature distributions or cardinality distributions of the data model (e.g., the $1^{st}$ and $2^{nd}$ simulated dataset). However, as expected, if there is significant overlap in both the feature distributions and the cardinality distributions, then model-based clustering performs poorly since there is not enough information to separate the data (e.g., the $3^{rd}$ simulated dataset).

III-D2 EM clustering with real data

In this experiment, we cluster the Texture dataset (described in section II-C2) using EM clustering with various iterations $N_{\mathrm{iter}}=\{1,...,10\}$. The performance is shown in Fig. 8. The best performance is Rand index of 0.86, which is equal the best performance of AP clustering using OSPA with $c=60$ (Fig. 5). Fig. 9 shows the learned distribution after 10 EM iterations.

IV Conclusion

This paper has detailed a non-parametric approach (based on set distance) and a model-based approach (based on random finite set) to the clustering problem for point pattern data (aka ‘bags’ or multiple instance data). Experiments with both simulated and real data indicate that, in the non-parametric method, the choice of distance has a big influence on the clustering performance. The experiments also indicate that the model-based method has better performance than the non-parametric method when there is little overlap between the feature distributions or cardinality distributions of the data model. However, as expected, if there is significant overlap in both the feature distributions and the cardinality distributions, then model-based clustering performs poorly since there is not enough discriminative information to separate the data.

Future research directions may include developing more complex RFS models such as RFSs with Gaussian mixture feature distribution which can capture better multiple modes feature data (e.g., section II-C2). Another promising development is adapting the proposed clustering algorithms into data stream mining – an emerging research topic dealing with rapidly and continuously generated data such as search or surveillance data [40].