Introduction

Clustering is one of the basic problems of data mining aimed to organizing a set of data objects into subsets (clusters) such that objects in a cluster are similar to one another, yet dissimilar to objects in other clusters. Similarity is commonly defined in terms of how close the objects are and is based on a specified distance metric.

The most fundamental method of clustering is partitioning, which organizes the objects of a set into several exclusive groups. More formally, given a set of n objects, a partitioning algorithm constructs k partitions of the data, where each partition represents a cluster and k ≤ n. The algorithm divides the data objects into k clusters. An object is assigned to a closest cluster based on the distance measure between the object and the cluster center. Then algorithm iteratively improves the within-cluster variation by computing the new cluster center using the objects assigned to the cluster in the previous iteration. After

This work was financially supported by the Ministry of education and science of the Russian Federation ("Research and development on priority directions of scientific-technological complex of Russia for 2014-2020" Federal Program, contract No. 14.574.21.0035).

this cluster centers are updated and all the objects are then reassigned using the new cluster centers. The iterations continue until the assignment is stable, that is, the clusters formed in the current round are the same as those formed in the previous round.

Partitioning clustering algorithms differ in a way of calculation cluster centers, e.g. k-Means [11] and k-Modes [5] algorithms uses mean and mode values of clustered objects respectively, whereas k-Medoids algorithm uses an object of clustered data set (called medoid ).

The Partition Around Medoids (PAM) [18] is a variation of k-Means, which is used in a wide spectrum of applications, e.g. text analysis, bioinformatics, intelligent transport systems, etc. The complexity of each iteration in the PAM algorithm is O(k(n − k) 2 ). For large values of n and k computations are very costly. That is why there are approaches to speed up k-Means and PAM algorithms by means of GPU [3,10]. At the same time there none for modern accelerators based on the Intel Many Integrated Core (MIC) [8] architecture. In spite of many recent developments for manycore platforms in data mining [13,15,23] and databases [1,7] there are few ones for the Intel MIC architecture.

In this paper we present a parallel version of PAM for MIC accelerators. The remaining part of the paper is organized as follows. Section 2 gives an overview of serial PAM algorithm and discusses related work. In section 3 we describe parallelization of PAM adapted for the Intel MIC architecture. The results of the experiments evaluating the algorithm are presented in section 4. Section 5 contains summary and directions for future research.

Background of the Research

Serial PAM Algorithm

To provide formal description of the PAM [9] algorithm we will use the following notation. Let O = {o 1 , o 2 , . . . , o n } is a set of objects to be clustered where each object is a tuple consisting of p real-valued attributes. Let k is the number of clusters, k n, and C = {c 1 , c 2 , . . . , c k } is a set of medoids, C⊂O, and

ρ : O × C → R is a distance metric.

The algorithm takes the form of a steepest ascent hill climber, using a simple swap neighbourhood operation. In each iteration medoid object c i and nonmedoid object o j are selected that produce the best clustering when their roles are switched. The objective function used is the sum of the distances from each object to the closest medoid:

E = n j=1 min 1≤i≤k ρ(c i , o j ). (1)

Algorithm 1 depicts PAM pseudocode. PAM consists of two phases, namely BUILD and SWAP. In the first phase an initial clustering is obtained by the successive selection of representative objects until k objects have been found. The first object c 1 is the one for which the sum of the distances to all other objects is as small as possible:

c 1 = arg min 1≤h≤n n j=1 ρ(o h , o j ).

(

)2

Object c 1 is the most centrally located in O set. Subsequently, at each step another object is selected, which decreases the objective function as much as possible. This object is the one for which the minimal distance to all selected medoids and distance to this object is as small as possible:

c 2 = arg min 1≤h≤n n j=1 min(ρ(c 1 , o j ), ρ(o h , o j )),(3)c 3 = arg min 1≤h≤n n j=1 min( min 1≤l≤2 (ρ(c l , o j )), ρ(o h , o j )),(4)

...

c k = arg min 1≤h≤n n j=1 min( min 1≤l≤k−1 (ρ(c l , o j )), ρ(o h , o j )). (5)

This process is continued until k objects have been found.

In the second phase of the algorithm, it is attempted to improve C (i.e. set of medoids) and therefore also to improve the clustering yielded by this set. Algorithm searches for a pair of objects (c min , o min ), which minimizes the objective function. This is done by considering all pairs of objects (c i , o h ) where c i is a medoid and o h is not a medoid. It is determined what effect is obtained on the objective function when a swap is carried out, i.e., when object c i is no longer selected as a medoid but object o h is. Let denote this effect as T ih , then minimum value of T min is achieved with (c min , o min ) pair. If T min > 0 then C set can not be improved so the algorithm stops.

Let us consider calculation of the T ih effect using the following notation. Let D = {d 1 , d 2 , . . . , d n } is a set of distances from each object to the closest medoid. Let S = {s 1 , s 2 , . . . , s n } is a set of distances from each object to second closest medoid. Let C jih is a contribution of non selected object o j to the effect T ih of a swap between c i and o h on the objective function. In this case T ih is the sum of the contributions C jih :

T ih = n j=1 C jih .(6)

Algorithm 2 [9] depicts pseudocode of calculating C jih .

Input : oj, ci, o h , dj, sj Output: C jih 1 if ρ(oj, ci) > dj and ρ(oj, o h ) > dj then 2 C jih ← 0 3 else if ρ(oj, ci) = dj then 4 if ρ(oj, o h ) < sj then 5 C jih ← ρ(oj, o h ) − dj 6 else 7 C jih ← sj − dj 8 end 9 else if ρ(oj, o h ) < dj then 10 C jih ← ρ(oj, o h ) − dj 11 end Fig. 2. Calculating C jih

Related Work

A significant amount of work has been done in the area of cluster analysis. The classical k-Means and k-Medoids algorithms was suggested in [5,11]. The original PAM algorithm was proposed in [9].

The research devoted to accelerating clustering algorithms using parallel hardware includes the following. In [6] FPGA and GPU implementations of k-Means are compared. Authors of [20] describe improvements of k-Means reducing data transfers between CPU and GPU. In [21] a technique improving data distribution among GPU threads in k-Means is suggested. k-Means implementation for Hadoop framework with GPUs is described in [22]. In [3] several clustering methods on GPU including k-Medoids are implemented. A GPU-based framework for clustering genetic data using k-Medoids described in [10].

In our opinion currently the potential of the Intel MIC accelerators for cluster analysis is underestimated. Paper [16] proposes modification of the DBSCAN density-based clustering algorithm for the Intel MIC architecture. In [19] a version of k-Means for CPU and Intel MIC heterogeneous architecture is presented, where authors used vectorization and sophisticated layout scheme to improve data locality. The contribution of this paper is technique of acceleration of the Partitioning Around Medoids clustering algorithm with the Intel Xeon Phi many-core coprocessor.

Parallel PAM Algorithm for MIC Accelerators

In this section we describe an approach to implementation of PAM algorithm for the Intel Xeon Phi coprocessor [17]. The Intel Xeon Phi coprocessor is an x86based SMP-on-a-chip with over fifty cores. It supports 4× hardware threads per core and contains 512-bit wide vector processor unit (VPU). Each core has two levels of cache memory: a 32 Kb L1 data cache, a 32 Kb L1 instruction cache, and a core-private 512 Kb unified L2 cache. The Intel Xeon Phi coprocessor is connected to other devices via the PCIe bus. Intel Xeon Phi coprocessor is based on Intel x86 architecture and it supports the same programming tools and models as a regular Intel Xeon processor. Our approach is based on the following principles.

Data parallelism and vectorization. Using OpenMP technology we perform simultaneous execution on multiple cores of the same function across the elements of a dataset. Most loops of the original PAM algorithm with arithmetic operations were implemented to provide conversion of such operations from scalar form to vector form to be effectively computed by the coprocessor's VPUs.

Our implementation strives to provide data locality as much as possible, i.e. the program uses data close to recently accessed locations. Since the coprocessor loads a chunk of memory around an accessed location into the cache, locations close to recently accessed locations are also likely to be in the cache so finally it increases algorithm's performance.

Algorithm 3 depicts PAM pseudocode adapted for use on the Intel Xeon Phi many-core coprocessor. The summary of parallel PAM subalgorithms is presented in Tab. 1. The PAM algorithm deals with a lot of data arrays which are not fit into Intel Xeon Phi L2 memory cache. We process data by chunks of L bytes to satisfy data locality requirement. It is recommended [8] to set L to 16 and try multiplying or dividing by 2 and use n divisible by L. In our work we use L = 32.

To improve performance we use precomputing technique by means of calculating distances between all objects of O set in advance. There is no need for

The PrepareDistanceMatrix subalgorithm initializes distance matrix (see Algorithm 4). Unlike in [9] we store matrix in full form (not in upper triangular form) to provide better data locality for the rest of subalgorithms. To achieve better performance of this subalgorithm we use tiling technique [8]. Tiling is a technique for improving data reuse in cache architectures. Cache architectures generally employ least recently used (LRU) methods to determine which data is evicted from the cache as new data is requested. Therefore, the longer data remains unused, the more likely it will be evicted from the cache and no longer available immediately when needed. Tiling the access pattern can exploit data that remains in the cache from recent, previous iterations.

Input : Set of objects O Output: Distance matrix M 1 parallel for oi such that 1 ≤ i ≤ n do 2 for j = 1 to n step L do 3 for k = 1 to p do 4 for l such that j ≤ l ≤ j + L do /* vectorized */ 5 m il ← m il + (oi[k] − o l [k]) 2 ; /*

The BuildMedoids subalgorithm implements BUILD phase (see Alg. 5) according to formulas (2)-( 5). The FindBestSwap subalgorithm implements SWAP phase (see Alg. 6). It checks all pairs of (c i , o h ) objects where c i is a medoid and o h is not a medoid, calculates the effect for each T ih swapping and returns the minimal one. T ih ; Fig. 6. SWAP phase SWAP phase executes many logical operations in (see Alg. 2). By this reason two versions of the PAM algorithm were implemented. PAM-1 executes more logical operations with lesser temporary data. PAM-2 executes lesser logical op-erations but stores more temporary data. Preparing matrix function and BUILD phase are the same in both versions.

Input : Distance matrix M Output: Set of medoids C 1 parallel for i = 1 to n do 2 if n j=1 mij is minimal then /* sum is vectorized */ 3 c1 ← oi;

Experimental Evaluation

To evaluate the developed algorithm we performed experiments on the hardware specified in Tab. 2. Experiments were performed on single precision data, the coprocessor was used in offload mode. We measured PAM runtime while varying number of clustered objects and investigated the influence of dataset properties on runtime of PAM subalgorithms. Experimental results for FCS Human dataset are introduced in Fig. 7(a). FCS Human dataset has large dimension so the most time is taken by calculation of distance matrix. Calculation of distance matrix on the Intel Xeon Phi is two times faster then on the Intel Xeon. There is no significant difference between PAM-1 and PAM-2 for this dataset.

Experimental results for Corel Image Histogram dataset are introduced in Fig. 7(b). Data dimension is small so preparing distance matrix does not require much time. PAM-1 shows similar performance on both CPU and Intel Xeon Phi. The PAM-2 algorithm is two times slower on the Intel Xeon than on the Intel Xeon Phi.

Experimental results for MixSim dataset are introduced in Fig. 7(c). Again PAM-1 shows similar performance on both CPU and the Intel Xeon Phi. PAM-2 on the Intel Xeon Phi shows best result on this dataset. Experimental results for Letter Recognition dataset are introduced in Fig. 7(d). PAM-1 shows the best result on the Intel Xeon. Both PAM-1 and PAM-2 shows similar results on the Intel Xeon Phi.

Intuitively PAM-2 is a better implementation for the Intel Xeon Phi. This suggestion is confirmed by experiments. In all tests PAM-2 is twice better on the Intel Xeon Phi than the Intel Xeon. In the same time PAM-1 is the best with the Intel Xeon only once. In other tests there is no significant difference.

To investigate this fact deeper we made more experiments to see contribution of every PAM subalgorithm in Fig. 8. Figures 8(a Experiments show that PAM performance depends on clustered data nature. The most complex thing for large dimension data is calculation of distance matrix. In case of small dimension data the rest of the PAM subalgorithms take significantly larger part of runtime than distance matrix calculation. BUILD phase is more effective on the Intel Xeon Phi. SWAP phase perform better on the Intel Xeon. PAM execution on Letter Recognition dataset requires more iterations than MixSim experiment. By this reason PAM-1 shows best result with Letter Recognition and PAM-2 shows best result with MixSim.

Conclusion

The paper has described a parallel version of Partitioning Around Medoids clustering algorithm for the Intel Xeon Phi many-core coprocessor. An optimized version of PAM for the Intel Xeon Phi coprocessor is introduced where OpenMP parallelizing technology, loop vectorization, tiling technique and efficient distance matrix computation for Euclidean metric are used. Algorithm stores data in continuous arrays and process data by chunks to achieve data locality for better performance.

Experimental results show effectiveness of suggested approach. Experiments show that PAM performance depends on clustered data nature. The most complex thing for large dimension data is calculation of distance matrix. In case of small dimension data the rest of the PAM subalgorithms take significantly larger part of runtime than distance matrix calculation. BUILD phase is more effective on the Intel Xeon Phi. SWAP phase perform better on the Intel Xeon. PAM-1 shows best result with Letter Recognition dataset and PAM-2 shows best result with MixSim.

As future work we plan to extend our research in the following directions: implement our algorithm for the cases of several coprocessors and cluster system based on nodes equipped with the Intel Xeon Phi coprocessor(s).