T. Kamishima and S. Akaho, "Efficient Clustering for Orders", In D. A. Zighed, S. Tsumoto, Z.
W. Ras, and H. Hacid, editors, Mining Complex Data, Vol.165 of Studies in Computational
Intelligence, chapter 15, pp. 261–280. Springer (2009)
DRAFT VERSION

Efﬁcient Clustering for Orders
Toshihiro Kamishima and Shotaro Akaho
National Institute of Advanced Industrial Science and Technology (AIST),AIST Tsukuba Central 2, Umezono 1–1–1, Tsukuba, Ibaraki, 305–8568 Japan,
mail@kamishima.net hhttp://www.kamishima.net/i and s.akaho@aist.go.jp

Abstract. Lists of ordered objects are widely used as representational forms.
Such ordered objects include Web search results or best-seller lists. Clustering isa useful data analysis technique for grouping mutually similar objects. To cluster
orders, hierarchical clustering methods have been used together with dissimilarities deﬁned between pairs of orders. However, hierarchical clustering methods
cannot be applied to large-scale data due to their computational cost in terms ofthe number of orders. To avoid this problem, we developed an k-o’means algorithm. This algorithm successfully extracted grouping structures in orders, andwas computationally efﬁcient with respect to the number of orders. However,it was not efﬁcient in cases where there are too many possible objects yet. We
therefore propose a new method (k-o’means-EBC), grounded on a theory of order statistics. We further propose several techniques to analyze acquired clustersof orders.

1Introduction

The term order indicates a sequence of objects sorted according to some property. Suchorders are widely used as representational forms. For example, the responses from Web
search engines are lists of pages sorted according to their relevance to queries. Bestseller lists, which are item-sequence sorted according to sales volume, are used on many
E-commerce sites.
Orders have also been exploited for sensory test of human respondents’ sensations,impressions, or preference. For such a kind of surveys, it is typical to adopt a scoring
method. In this method, a respondents’ sensation is measured using a scale on which
extremes are represented by antonymous words. One example is a ﬁve-point-scale onwhich 1 and 5 indicate don’t prefer and prefer, respectively. If one very much prefers
an apple, he/she rates the apple as 5. Though this scoring method is widely used, it
is not the best way for all types of sensory test. For example, as pointed out in [1],a trained expert, e.g., a wine taster, can maintain a consistent mapping from his/her
sensation level to rating score throughout a given session. However, users’ mappings
generally change for each response, especially if the intervals between responses are
long. Hence, even if two respondents rate the same item at the same score, their true
degrees of sensation may not be the same. When effects of such demerits cannot beignored, a ranking method is used. In this method, respondents show their degree of
sensation by orders, i.e., object sequences according to the degree of a target sensation.
In this case, respondents’ sensation patterns are represented by orders, and analysis
techniques for orders are required.

2

Orders are also useful when the absolute level of observations cannot be calibrated.
For example, when analyzing DNA microarray data, in order that the same ﬂuorescelevel represents the same level of gene expression, experimental conditions must be
calibrated. However, DNA databases may consist of data sampled under various conditions. Even in such cases, the higher level of ﬂuoresce surely corresponds to the higher
level of gene expression. Therefore, by treating the values in the microarray data as ordinal values, non-calibrated data would be processed. Fujibuchi et al. adopted such useof orders in searching a gene expression database for similar cell types [2].
And clustering is the task of partitioning a sample set into clusters having the properties of internal cohesion and external isolation [3]. This method is a basic tool for
exploratory data analysis. Clustering methods for orders are useful for revealing thegroup structure of data represented by orders such as those described above.
To cluster a set of orders, classical clustering has been mainly used [4, chapter 2].
In these studies, clustering methods were applied to ordinal data of a social survey,
sensory test, etc. These data sets have been small in size; the number of objects to besorted and the length of orders are at most ten, and the number of orders to be clusteredare at most thousands. This is because an scoring method has been used to acquire
responses for a large-scale survey. Responses can easily be collected by requesting for
respondents to mark on rating scales that are printed on paper questionnaire forms.
On the other hand, using printed questionnaire forms is not appropriate for ranking
method, because respondents must rewrite entire response orders when they want tocorrect them. Therefore, in a ranking method, respondents generally reply by sorting
real objects. For example, respondents are requested to sort glasses of wine according
to their preference. However, it would be costly to prepare so many glasses. Due to thisreason, ranking method has been used for a small-scale survey, even if its advantage
to an scoring method is known as described above. But now, adoption of computerinterface clear this obstacle in using a ranking method. Respondents can sort virtual
objects instead of real objects. Further, methods to implicitly collect preference orders
have proposed [5, 6]. These technical progress has made it easier to collect the large
number of ordinal data.
We can now collect a large-scale data that consist of orders. However, current techniques for clustering orders are not fully scalable. For example, to cluster a set of orders,dissimilarities are ﬁrst calculated for all pairs of orders, and agglomerative hierarchicalclustering techniques are applied. This approach is computationally inefﬁcient, because
computational cost of agglomerative hierarchical clustering is O(N 2 log(N ))) under
non-Euclidean metric [7], where N is the number of orders to be clustered. To alleviate
this inefﬁciency in terms of N , we proposed a k-means-type algorithm k-o’means in
our previous work [8]. The computational complexity was reduced to O(N ) in terms ofthe number of orders. Though this method successfully extracted a grouping structure ina set of orders, it was not efﬁcient yet, if the number of possible objects to be sorted waslarge. In this paper, to alleviate this inefﬁciency, we propose a new method, k-o’meansEBC. Note that EBC means Expected Borda Count, which is a classic method to ﬁnd
an order so as to be as concordant as possible with a given set of orders. And incompleteness in orders are processed based on a theory of order statistics. Additionally, we
propose several methods for interpreting the clusters of orders.

3

We formalize this clustering task in Section 2. Our previous and new clustering
methods are presented in Section 3. The experimental results are shown in Sections 4
and 5. Section 6 summarizes our conclusions.

2

Clustering Orders

In this section, we formalize the task of clustering orders. We start by deﬁning our basic notations regarding orders. An object, entity, or substance to be sorted is denotedby xj . The universal object set, X ∗ , consists of all possible objects, and L∗ is deﬁned
as |X ∗ |. The order is denoted by O = xa Â · · · Âxj Â · · · Âxb . Note that subscript jof x doesn’t mean “The j-th object in this order,” but that “The object is uniquely indexed by j in X ∗ .” The order x1 Âx2 represents “x1 precedes x2 .” An object set X(Oi )or simply Xi is composed of all objects in the order Oi . The length of Oi , i.e., |Xi |,is shortly denoted by Li . An order of all objects, i.e., Oi s.t. X(Oi )=X ∗ , is called
a complete order; otherwise, the order is incomplete. Rank, r(Oi , xj ) or simply rij ,is the cardinal number that indicates the position of the object xj in the order Oi . For
example, for Oi =x1 Âx3 Âx2 , r(Oi , x2 ) or ri2 is 3. Two orders, O1 and O2 , are concordant if ordinal relations are consistent between any object pairs commonly contained inthese two orders; otherwise, they are discordant. Formally, for two orders, O1 and O2 ,consider an object pair xa and xb such that xa , xb ∈X1 ∩X2 , xa 6=xb . We say that theorders O1 and O2 are concordant w.r.t. xa and xb if the two objects are placed in thesame order, i.e., (r1a − r1b )(r2a − r2b ) ≥ 0; otherwise, they are discordant. Further,
O1 and O2 are concordant if O1 and O2 are concordant w.r.t. all object pairs such thatxa , xb ∈X1 ∩X2 , xa 6=xb .
A pair set Pair(Oi ) is composed of all the object pairs xa Âxb , such that xa precedes xb in the order Oi . For example, from the order O1 =x3 Âx2 Âx1 , three object
pairs, x3 Âx2 , x3 Âx1 , and x2 Âx1 , are extracted. For a set of orders S, the Pair(S) is
composed of all pairs in Pair(Oi ) of Oi ∈ S. Note that if the same object pairs arecontained in numbers of Pair(Oi ), these pairs are multiply added into the Pair(S). For
example, if the same object pairs x1 Âx2 are extracted from O5 and O7 in S, both twoordered pairs x1 Âx2 are multiply included in Pair(S).
The task of clustering orders is as follows. A set of sample orders, S = {O1 , O2 , . . . ,ON }, N ≡ |S|, is given. Note that sample orders may be incomplete, i.e., Xi 6=Xj , i 6=
j. In addition, Oi and Oj can be discordant. The aim of clustering is to divide the S
into a partition. The partition, π = {C1 , C2 , . . . , CK }, K = |π|, is a set of all clusters. Clusters are mutually disjoint and exhaustive, i.e., Ck ∩ Cl = ∅, ∀k, l, k 6= l andS = C1 ∪ C2 ∪ · · · ∪ CK . Partitions are generated such that orders in the same cluster
are similar (internal cohesion), and those in different clusters are dissimilar (external
isolation).
2.1 Similarity Between Two Orders
Clusters are deﬁned as a collection of similar orders; thus, the similarity measures between two orders are required. Spearman’s ρ [9, 4] is one such measure, signifying the4correlation between ranks of objects. The ρ between two orders, O1 and O2 , consistingof the same objects (i.e., X ≡ X(O1 ) = X(O2 )) is deﬁned as:()()∑
xj ∈X r1j − r̄1 r2j − r̄2ρ = √∑
(
)2 √∑
(
)2 ,
xj ∈X r1j −r̄1
xj ∈X r2j −r̄2∑where r̄i = (1/L) xj ∈X rij , L=|X|. If no tie in rank is allowed, this can be calculatedby the simple formula:
ρ=1−

6

(

∑

xj ∈X r1j − r2j
L3 − L)2.

(1)

The ρ becomes 1 if the two orders are concordant, and −1 if one order is the reverse of
the other order. Observing Equation (1), this similarity depends only on the term∑2dS (O1 , O2 ) =(r1j − r2j ) .
(2)
xj ∈X

This is called Spearman’s distance. If two or more objects are tied, we give the same
midrank to these objects [4]. For example, consider an order x5 Âx2 ∼x3 (“∼” denotes
a tie in rank), in which x2 and x3 are ranked at the 2nd or 3rd positions. In this case, themidrank 2.5 is assigned to both objects.
Another widely used measure of the similarity of orders is Kendall’s τ . Intuitively,
this is deﬁned as the number of concordant object pairs subtracted by that of discordantpairs, and then it is normalized. Formally, Kendall’s τ is deﬁned as∑()
1
τ = L(L−1)/2
sgn (r1a −r1b )(r2a −r2b ) ,(3)xa Âxb ∈Pair(O1 )where sgn(x) is a sign function that takes 1 if x>0, 0 if x=0, and −1 otherwise. Manyother types of similarities between orders have been proposed (see [4, chapter 2]), but
the above two are widely used and have been well studied.
In this paper, we adopt Spearman’s ρ rather than Kendall’s τ because of the following reasons: First, these two measures have similar properties. Both measures ofsimilarities between two random orders asymptotically follow normal distribution as
the length of the orders grows. Additionally, these are highly correlated, because thedifference between the two measures is bounded by Daniels’ inequality [9]:−1 ≤2(L + 1)3(L + 2)τ−
ρ ≤ 1.
L−2
L−2

Second, Spearman’s ρ can be calculated more quickly. All of the object pairs have to
be checked to derive Kendall’s τ , so O(L2 ) time is required. In the case of Spearman’s ρ, the most time consuming task is sorting objects to decide their ranks; thus, thetime complexity is O(L log L). Further, the central orders under Spearman distance istractable, but the derivation under Kendall’s distance is NP-hard [10].

5

For the clustering task, distance or dissimilarity is more useful than similarity. We
deﬁned a dissimilarity between two orders based on ρ:dρ (O1 , O2 ) = 1 − ρ(O1 , O2 ).

(4)

Since the range of ρ is [−1, 1], this dissimilarity ranges [0, 2]. This dissimilarity becomes 0 if the two orders are concordant.

3

Methods

Here, we describe exiting clustering methods and our new clustering method.
3.1

Hierarchical Clustering Methods

In the literature of psychometrics, questionnaire data obtained by a ranking method havebeen processed by traditional clustering techniques [4]. First, for all pairs of orders in
S, the dissimilarities in Section 2.1 are calculated, and a dissimilarity matrix for S is
obtained. Next, this matrix can be clustered by standard hierarchical clustering methods,
such as the group average method. In these survey researches, the size of the processeddata set is rather small (N < 1000, L∗ < 10, Li < 10). Therefore, hierarchicalclustering methods could cluster order sets, even though the time complexity of thesemethods is O(N 2 log(N )) under non-Euclidean metric [7] and is costly. However, thesemethod cannot be applied to a large-scale data, due to their computational cost.
Additionally, when the number of objects, L∗ , is large, it is hard for respondents tosort all objects in X ∗ . Therefore, sample orders are generally incomplete, i.e., X(Oi ) ⊂X ∗ , the dissimilarities cannot be calculated because the dissimilarity measures are deﬁned between two orders consisting of the same objects. One way to deal with incomplete orders is to introduce the notion of an Incomplete Order Set (IOS)1 [4], which
is deﬁned as a set of all possible complete orders that are concordant with the given
incomplete order. Given the incomplete order O that consists of the object set X, an
IOS is deﬁned as
ios(O) = {Oi∗ |Oi∗ is concordant with O, X(Oi∗ ) = X ∗ }.
This idea is not ﬁt for large-scale data sets because the size of the set is (L∗ !/L!), which
grows exponentially in accordance with L∗ . Additionally, there are some difﬁculties indeﬁning the distances between the two sets of orders. One possible deﬁnition is to adopt
the arithmetic mean of the distances between orders in each of the two sets. However,this is not distance because d(iosa , iosa ) may not be 0. Therefore, more complicateddistance, i.e., Hausdorff distance, has to be adopted.
Since the above IOS cannot be derived for a large-scale data set, we adopted thefollowing heuristics in this paper. In such cases, the dissimilarity between the orders is1

In [4], this notion is referred by the term incomplete ranking, but we have adopted IOS to insist
that this is a set of orders.

6determined based on the the objects included in both. Take, for example, the following
two orders:O1 =x1 Âx3 Âx4 Âx6 ,O2 =x5 Âx4 Âx3 Âx2 Âx6 .

From these orders, all objects that are not included in both orders are eliminated. The
generated orders become:O10 = x3 Âx4 Âx6 ,O20 = x4 Âx3 Âx6 .

The ranks of objects in these orders are:r(O10, x3 )=1, r(O10, x4 )=2, r(O10, x6 )=3;r(O20, x3 )=2, r(O20, x4 )=1, r(O20, x6 )=3.
Consequently, the Spearman’s ρ becomes()
6 (1−2)2 + (2−1)2 + (3−3)2
ρ=1−= 0.5.
33 − 3
If no common objects exists between the two orders, ρ = 0 (i.e., no correlation).
3.2

k-o’means-TMSE (Thurstone Minimum Square Error)

In [8], we proposed a k-o’means algorithm as a clustering method designed to processorders. To differentiate our new algorithm described in detail later, we call it by a ko’means-TMSE algorithm.

Algorithm k-o’means(S, K, maxIter)
S = {O1 , . . . , ON }: a set of orders
K: the number of clusters
maxIter: the limit of iteration times
1) S is randomly partitioned into a set of clusters: π = {C1 , . . . , CK },π 0 := π, t := 0.
2) t := t + 1, if t > maxIter goto step 6.
3) for each cluster Ck ∈ π,derive the corresponding central order Ōk .
4) for each order Oi in S,
assign it to the cluster: arg minCk d(Ōk , Oi ).
5) if π = π 0 then goto step 6; else π 0 := π, goto step 2.
6) output π.
Fig. 1. k-o’means algorithm

A k-o’means-TMSE in Figure 1 is similar to the well-known k-means algorithm[11]. Speciﬁcally, an initial cluster is reﬁned by the iterative process of estimating

7

new cluster centers and the re-assigning of samples. This process is repeated until nochanges in the cluster assignment is detected or the pre-deﬁned iteration time is reached.
However, different notions of dissimilarity and cluster centers have been used to handle
orders. For the dissimilarity d(Ōk , Oi ), equation (4) was used in step 4. As a cluster
center in step 3, we used the following notion of a central order [4]. Given a set oforders Ck and a dissimilarity measure between orders d(Oa , Ob ), a central order Ōk is
deﬁned as the order that minimizes the sum of dissimilarities:∑d(O, Oi ).
(5)Ōk = arg minO

Oi ∈CkNote that the order Ōk consists of all the objects in Ck , i.e., XCk = ∪Oi ∈Ck X(Oi ). The
dissimilarity d(Ōk , Oi ) is calculated over common objects as in Section 3.1. However,
because Xi ⊆ X(Ōk ), the dissimilarity can always be calculated over Xi . Unfortunately, the optimal central order is not tractable except for a special cases. For example,
if using a Kendall distance, the derivation of central orders is NP-hard even if all sample
orders are complete [10].
Therefore, many approximation methods have been developed. However, to use asa sub-routine in a k-o’means algorithm, the following two constraints must be satisﬁed.
First, the method must deal with incomplete orders that consist of objects randomly
sampled from X ∗ . In [12], they proposed a method to derive a central order of top klists, which are special kinds of incomplete orders. Top k list is an order that consistsof the most preferred k objects, and the objects that are not among the top k list are
implicitly ranked lower than these k objects. That is to say, the top k objects of a hidden
complete order are observed. In our case, objects are randomly sampled, and such a
restriction is not allowed. Second, the method should be executed without using iterative
optimization techniques. Since central orders are derived K times in each loop of the
k-o’means algorithm, the derivation method of central orders would seriously affectefﬁciency if it adopts the iterative optimization.
To our knowledge, the method satisfying these two constraints is the following one
to derive the minimum square error solution under a generative model of Thurstone’s
law of comparative judgment [13]. Because we used this method to derive central orders, we call this clustering algorithm by the k-o’means-TMSE (Thurstone Minimum
Square Error) algorithm. We describe this method for deriving central orders. First, theprobability Pr[xa Â xb ] is estimated. The pair set of Pair(Ck ) in Section 2 is generated
from Ck in step 3 of k-o’means-TMSE. Next, we calculate the probabilities for every
pair of objects in Ck :Pr[xa Â xb ] =|xa Â xb | + 0.5,|xa Â xb | + |xb Â xa | + 1

where |xa Â xb | is the number of the object pairs, xa Â xb , in the Pair(Ck ). These
probabilities are applied to a model of Thurstone’s law of comparative judgment. Thismodel assumes that scores are assigned to each object xl , and an order is derived bysorting according to these scores. Scores follow a normal distribution; i.e., N (µl , σ),where µl is the mean score of the object xl , and σ is a common constant standard

8

deviation. Based on this model, the probability that object xa precedes the xb is∫ ∞∫ tt − µau − µbPr[xa Âxb ] =φ()φ()du dt
σ
σ
−∞
−∞
(
)
µa − µb
√
=Φ
,
2σ

(6)

where φ(·) is a normal distribution density function, and Φ(·) is a normal cumulative
distribution function. Under the minimum square error criterion of this model [14], µ0l ,which is a linearly transformed image of µl , is analytically derived as
µ0l =

∑
(
)
1
Φ−1 Pr[xl Â x] ,|XCk |

(7)

x∈XCk

∪
where XCk = Oi ∈Ck Xi . The value of µ0l is derived for each object in XCk . Finally,
the central order Ōk can be derived by sorting according to the corresponding µ0l . Because the resultant partition by k-o’means-TMSE is dependent on the initial cluster, thisalgorithm is run multiple times, randomly changing the initial cluster; then, the partition
minimizing the following total error is selected:∑ ∑d(Oi , Ōk ).
(8)
Ck ∈π Oi ∈Ck

This k-o’means-TMSE could successfully ﬁnd the cluster structure in a set of incomplete orders due to the following reason: Because the dissimilarity in Section 3.1 wasmeasured between two orders, the precision of the dissimilarities was unstable. On theother hand, in the case of k-o’means-TMSE, central orders are calculated based on
the |Ck | orders. |Ck | is generally much larger than two, and much more information
is available; thus, the central order can be stably calculated. The dissimilarity between
the central orders and each sample order can be stably measured, too, because all ofobjects in a sample order always exist in the corresponding central order and so the full
information in the sample orders can be considered.
However, the k-o’means-TMSE is not so efﬁcient in terms of time and memory
complexity. Time or memory complexity in N and K is linear, and these are efﬁcient.
However, complexity in terms of L∗ is quadratic, and further, the constant factor israther large due to the calculation of the inverse function of a normal distribution. Dueto this inefﬁciency, this algorithm cannot be used if L∗ is large. To overcome this inefﬁciency, we propose a new method in the next section.
3.3

k-o’means-EBC (Expected Borda Count)

To improve efﬁciency in computation time and memory requirement, though we used
the k-o’means framework in Figure 1 and the dissimilarity measure dρ of equation (4)in step 4 of Figure 1, we employed other types of derivation procedures for the central
orders.
Below, we describe this derivation method for a central order Ōk of a cluster Ck instep 3 of Figure 1. We call this the Expected Borda Count(EBC) method, and our new

9clustering method is called a k-o’means-EBC algorithm. The Borda Count method isused to derive central orders from complete orders; we modiﬁed this so as to make itapplicable to incomplete orders. The Borda Count method [15] was originally developed for determining the order of candidates in an election from a set of ranking votes.
A set of complete orders, Ck , is given. First, for each object xj in X ∗ , the vote count iscalculated:∑ ()
vote(xj ) =L∗ − rij + 1 .
Oi ∈CkThen, a central order is derived by sorting objects xj ∈ X ∗ in descending order ofvote(xj ). Clearly, this method is equivalent to sorting the objects in ascending order ofthe following mean ranks:
1 ∑
rij .
(9)
r̄j =
|Ck |
Oi ∈CkIf all sample orders are complete and Spearman’s distance is used, it is known that thecentral order derived by the above Borda Count optimally minimizes Equation (5) [4,
theorem 2.2].
Because all sample orders are complete, Spearman’s distance is proportional to the
distance dρ . Therefore, even in the case that dρ is used as dissimilarity, the optimal
central order can be derived by this Borda Count method. This optimal central order
can also be considered as a maximum likelihood estimator of the Mallows-θ model[16]. The Mallows-θ model is a distribution model of the complete order O, and isdeﬁned as
Pr[O; O0 , θ] ∝ exp(θdS (O0 , O)),
(10)
where the parameters θ and O0 are called a dispersion parameter and a modal order,
respectively.
Unfortunately, this original Borda Count method cannot be applied to incomplete
orders. To cope with incomplete orders, we must show the facts known in the order
statistics literature. First, we assume that there is hidden complete order Oh∗ which is
randomly generated. A sample order Oi ∈ Ck is generated by selecting objects fromthis Oh∗ uniformly at random. That is to say, from a universal object set X ∗ , Li objects
are sampled without replacement; then, Oi is generated by sorting these objects so as
to be concordant with Oh∗ . Now we are given Oi generated through this process. In thiscase, the complete order Oh∗ follows the distribution:{
Pr[Oh∗ |Oi ] =

Li !
L∗ !

0

if Oh∗ and Oi are concordant,
otherwise.

(11)

Based on the theory of order statics from a without-replacement sample [17, section 3.7],
if an object xj is contained in Xi , the conditional expectation of ranks of the object xj
in the order Oh∗ given Oi is
E[rj∗ |Oi ] = rij

L∗ + 1,
Li + 1

if xj ∈ Xi ,

(12)

10

where the expectation is calculated over all possible complete orders, Oh∗ , and rj∗ ≡r(Oh∗ , xj ). If an object xj is not contained in Xi , the object is at any rank in the hidden
complete order uniformly at random; thus, an expectation of ranks is
E[rj∗ |Oi ] =1 ∗(L + 1),
2

if xj ∈/ Xi .

(13)

Next, we turn to the case where a set of orders, Ck , consists of orders independently
generated through the above process. Each Oi ∈ Ck is ﬁrst converted to a set of all complete orders; thus, the total number of complete orders is L∗ !|Ck |. For each complete
order, we assign weights that follow equation (11). By the Borda Count method, an
optimal central order for these weighted complete orders can be calculated. The meanrank of xj (equation (9)) for these weighted complete orders is
E[r̄j ] =

1 ∑
|Ck |

∑

Pr[Oh∗ |Oi ]r(Oh∗ , xj )∗ ∈S(L∗ )Oi ∈Ck Oh=

1 ∑E[rj∗ |Oi ],
|Ck |

(14)

Oi ∈Ckwhere S(L∗ )2 is a set of all complete orders. A central order is derived by sorting
objects xj ∈ XCk in ascending order of the corresponding E[r̄j ]. Since objects are
sorted according to the means of expectation of ranks, we call this method an Expected
Borda Count (EBC).
A central order derived by an EBC method is optimal if the distance d(Oi , Ōk ) ismeasured by∑(15)
Pr[Oh∗ |Oi ]dS (Oh∗ , Ōk ).
∗ ∈S(L∗ )Oh

Hence, in step 4 of Figure 1, not dρ , but this equation (15) should be used. However, it is intractable to compute equation (15), because its computational complexity is O(L∗ (L∗ !/Li !)). Therefore, we adopt dρ , and it empirically performed well, as
is shown later. Furthermore, if all sample orders are complete, dρ is compatible withequation (15). Note that we also tried∑2d(Ō, Oi )(r(Ōk , xj ) − E[rj∗ |Oi ]) ,xj ∈X ∗

but empirically, it performed poorly.
The time complexity of a k-o’means-EBC is()
O K max(N L̄ log(L̄), L∗ log L∗ ) ,(16)

where L̄ is the mean of Li over S. First, in step 3 of Figure 1, the K central orders
are derived. For each cluster, O((N/K)L̄) time is required for the means of expected
ranks and O(L∗ log L∗ ) time for sorting objects. Hence, the total time required for2

S(L∗ ) is equivalent to a permutation group of order L∗

11deriving K central orders is O(max(N L̄, KL∗ log L∗ )). Second, in step 4, N orders
are classiﬁed into K clusters. Because O(L̄ log L̄) time is required for calculating one
dissimilarity, O(N L̄ log(L̄)K) time is required in total. The number of iterations isconstant. Consequently, the total complexity becomes equation (16).
Note that the uniformity assumption of missing objects might look too strong. However, in the case of a questionnaire survey by ranking methods, the objects to be ranked
by respondents can be controlled by surveyors.
Further, if all the sample orders are ﬁrst converted into the expected rank vectors,∗hE[r1∗ |Oi ], . . . , E[rL
∗ |Oi ]i, then an original k-means algorithm is applied to these vectors. One might suppose that this k-means is equivalent to our k-o’means-EBC, but thisis not the case. A k-means is different from this k-o’means-EBC in terms of the derivation of centers; In the k-means case, the mean vectors of the expected ranks are directly
used as cluster centers; in a k-o’means case, these means are sorted and converted torank values. Therefore, in the k-means case, the centers that correspond to the same
central orders are simultaneously kept during clustering. For example, two mean rank
vectors h1.2, 1.5, 4.0i and h1, 5, 10i, correspond to the same central order x1 Â x2 Â
x3 , but these two vectors are not differentiated. On the other hand, in a k-o’means-EBC
algorithm, they are considered as equivalent, and thus we suppose that the k-o’meansEBC algorithm can ﬁnd the cluster structure reﬂecting the ordinal similarities among
data.

4

Experiments on Artiﬁcial Data

We applied the algorithms in Section 3 to two types of data: artiﬁcially generated data
and real questionnaire survey data. In the the former experiment, we examined the characteristics of each algorithm. In the latter experiment of the next section, we analyzed aquestionnaire survey data on preferences in sushi.
4.1

Evaluation Criteria

The evaluation criteria for partitions was as follows. The same object set was divided
into two different partitions: a true partition π ∗ and an estimated one π̂. To measure thedifference of π̂ from π ∗ , we adopted the ratio of information loss (RIL) [18], whichis also called the uncertainty coefﬁcient in numerical taxonomy literature. The RIL is
the ratio of the information that is not acquired to the total information required for
estimating a correct partition. This criterion is deﬁned based on the contingency tablefor indicator functions [11]. The indicator function I((xa , xb ), π) is 1 if an object pair(xa , xb ) are in the same cluster; otherwise, it is 0. The contingency table is a 2 × 2
matrix consisting of elements, ast , that are the number of object pairs satisfying thecondition I((xa , xb ), π ∗ )=s and I((xa , xb ), π̂)=t, among all the possible object pairs.
RIL is deﬁned as∑1 ∑1 asta·t
s=0
t=0 a log2 a
RIL =
(17)
∑1 as· ·· a·· st ,s=0 a·· log2 as·∑∑∑where a·t = s ast , as· = t ast , and a·· = s,t ast . The range of the RIL is [0, 1];
it becomes 0 if two partitions are identical.

12

4.2 Data Generation Process
Test data were generated in the following two steps: In the ﬁrst step, we generated the
K orders to be used as central orders. One permutation (we called it a pivot) consisting of all objects in X ∗ was generated. The other K − 1 centers were generated bytransforming this pivot. Two adjacent objects in the pivot were randomly selected andexchanged. This exchange was repeated at speciﬁed times. By changing the number of
exchanges, the inter-cluster isolation could be controlled.
In the second step, for each cluster, constituent orders were generated. From thecentral order, Li objects were randomly selected. These objects were sorted so as tobe concordant with the central order. Again, two adjacent object pairs were randomly
exchanged. By changing the number of times that objects were exchanged, the intracluster cohesion could be controlled. Note that the sizes of clusters are equal.

Table 1. Parameters of experimental data
1) the number of sample orders: N = 10002) the length of the orders: Li = 10
3) the total number of objects: L∗ = 10, 1004) the number of clusters: K = {2, 5, 10}5) the inter-cluster isolation: {0.5, 0.2, 0.1, 0.001}
6) the intra-cluster cohesion: {1.0, 0.999, 0.99, 0.9}

The parameters of the data generator are summarized in Table 1. The differences
between orders cannot be statistically tested if Li is too short; on the other respondents
cannot sort too many objects. Therefore, we set the order length to Li = 10. Param 1–
2 are common for all the data. The total number of objects (Param 3) is set to 10
or 100. All the sample orders are complete if L∗ = 10, and these are examined in
Section 4.3. We examine the incomplete case (L∗ = 100) in Section 4.4. Param 4
was the number of clusters. It is difﬁcult to partition if this number is large, since thesizes of the clusters then decrease. Param 5 was the inter-cluster isolation that could
be tuned by the number of times that objects are exchanged in the ﬁrst step of the data
generation process. This isolation is measured by the probability that the ρ between apivot and another central order is smaller than that between a pivot and a random order.
The larger the isolation, the more easily clusters are separated. Param 6 was the theintra-cluster cohesion indicating the number of times that objects are exchanged in thesecond step of the data generation process. This cohesion is measured by the probability
that the ρ between the central order and a sample one is larger than that between thecentral order and a random one. The larger the cohesion, the more easily a cluster could
be detected.
For each setting, we generated 100 sample sets. For each sample set, we ran thealgorithms ﬁve times using different initial partitions; then the best partition in terms ofEquation (8) was selected. Below, we show the means of RIL over these sets.

13
RIL

RIL

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

AVE
TMSE
EBC

0
0.5

0.2

0.1

AVE
TMSE
EBC

0
0.001

1.0

0.999

0.99

0.9

(a) K = 2
RIL

RIL

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

AVE
TMSE
EBC

0
0.5

0.2

0.1

AVE
TMSE
EBC

0
0.001

1.0

0.999

0.99

0.9

(b) K = 5
RIL

RIL

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

AVE
TMSE
EBC

0
0.5

0.2

0.1

AVE
TMSE
EBC

0
0.001

1.0

0.999

0.99

0.9

(c) K = 10
Fig. 2. Experimental results on artiﬁcial data of complete orders
NOTE: The left charts show the variation of RIL in the inter-cluster isolation when theintra-cluster cohesion is ﬁxed to 0.999. The right charts show the variation of RIL in theintra-cluster cohesion when the inter-cluster isolation is ﬁxed to 0.2.

4.3

Complete Order Case

We analyzed the characteristics of the methods in Section 3 by applying these to artiﬁcial data of complete orders. The two k-o’means methods were abbreviated to TMSE
and EBC, respectively. Additionally, a group average hierarchical clustering method
using dissimilarity as described in Section 2.1 was tested, and we denoted this result by
AVE. The experimental results on artiﬁcial data of complete orders (i.e, L∗ = 10) are
shown in Figure 2. In Figures 2(a), (b), and (c), the means of RIL are shown in cases of
K = 2, 5, and 10, respectively. The left three charts show the variation of RIL in the14inter-cluster isolation when the intra-cluster cohesion is ﬁxed to 0.999. The right three
charts show the variation of RIL in the intra-cluster cohesion when the inter-cluster
isolation is ﬁxed to 0.2.
As expected, the more inappropriate clusters were obtained when the inter-cluster
isolation or the intra-cluster cohesion decreased and the number of clusters increased.
We begin with the variation of estimation performance according to the decrease of
intra-cluster cohesion. If the cohesion is 1, sample orders are exactly concordant with
their corresponding true central orders. In this trivial case, the AVE method succeeds
almost perfectly in recovering the embedded cluster structure. Because the dissimilarities between sample orders are 0 if and only if they are in the same cluster, this method
could lead to perfect clusters. Though both the EBC and TMSE methods found almost
perfect clusters in the K=2 case, the performance gradually worsened when K increased. In a k-o’means clustering, a central order is chosen from the ﬁnite set, S(L∗ ).
This is contrasted to the fact that a domain of centers is an inﬁnite set in the clusteringof real value vectors. Hence, the central orders of two clusters happen to agree, and oneof these clusters is diminished during execution of the k-o’means. As the increase ofK, clusters are merged with higher probability. For example, in the EBC case, whenK = 2 and K = 5, clusters are merged in 7% and 35% of the trials, respectively.
Such occurrence of merging degrades the ability of recovering clusters. As the cohesion increases, the performance of AVE became more drastically worse than the othertwo methods. Furthermore, in terms of the inter-cluster isolation, the performance of
AVE became drastically worse as K increased, except for the trivial case in which thecohesion was 1. In the AVE method, the determination to merge clusters is based on
local information, that is, a pair of clusters. Hence, the chance that orders belonging todifferent clusters would happen to be merged increases when orders are broadly distributed. When comparing EBC and TMSE, these two methods are almost completely
the same.

4.4

Incomplete Order Case

We move to the experiments on artiﬁcial data of incomplete orders (i.e, L∗ = 100). The
results are shown in Figure 3. The meanings of the charts are the same as in Figure 2.
TMSE was slightly better than EBC when K = 2 and K = 5 cases; but EBC
overcame TMSE when K = 10. AVE was clearly the worst. We suppose that thisadvantage of the k-o’means is due to the fact that the dissimilarities between order
pairs could not be measured precisely if the number of objects commonly included inthese two orders is few. Furthermore, the time complexity of AVE is O(N 2 log N ),while the k-o’means algorithms are computationally more inexpensive as in Equation (16). When comparing TMSE and EBC, TMSE would be slightly better. However, in terms of time complexity, TMSE’s O(N L∗ max(L∗ , K)) is much worse than
EBC’s O(K max(N L̄ log(L̄), L∗ log L∗ ) if L∗ is large. In addition, while the required
memory for TMSE is O(L∗ 2 ), EBC demands far less O(KL∗ ). Therefore, it is reasonable to conclude that k-o’means-EBC is an efﬁcient and effective method for clusteringorders.

15
RIL

RIL

1
0.8

1
0.8

0.6

0.6

AVE
TMSE
EBC

0.4

0.4

0.2

0.2

0

AVE
TMSE
EBC

0
0.5

0.2

0.1

0.001

1.0

0.999

0.99

0.9

(a) K = 2
RIL

RIL

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

AVE
TMSE
EBC

0
0.5

0.2

0.1

AVE
TMSE
EBC

0
0.001

1.0

0.999

0.99

0.9

(b) K = 5
RIL

RIL

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

AVE
TMSE
EBC

0
0.5

0.2

0.1

AVE
TMSE
EBC

0
0.001

1.0

0.999

0.99

0.9

(c) K = 10
Fig. 3. Experimental results on artiﬁcial data of incomplete orders
NOTE: See note in Figure 2.

5 Experiments on Real Data
We applied our two k-o’means to questionnaire survey data, and proposed a method tointerpret the acquired clusters of orders.
5.1 Data Sets
Since the notion of true clusters is meaningless for real data sets, we used the k-o’means
as tools for exploratory analysis of a questionnaire survey of preference in sushi (aJapanese food). This data set was collected by the procedure in our previous works [19,168]. In this data set, N = 5000, Li = 10, and L∗ = 100; in the survey, the probability
distribution of sampling objects was not uniform as in equation (11). We designed it sothat the more frequently supplied sushi in restaurants were more frequently shown to
respondents. Objects were selected independently with probabilities ranging from 3.2%
to 0.13%. Therefore, the assumption of the uniformity of the sampling distribution, introduced by the EBC method, was violated. The best result in terms of Equation (8)ware selected from 10 trials. The number of clusters, K, was set to 2. Note that responses of both authors were clustered into Cluster 1.
5.2

Qualitative Analysis of Order Clusters

Table 2. Relations between clusters and attributes of objects
Attribute
A1
A2
A3
A4

Cluster 1
Cluster 2
EBC
TMSE
EBC
TMSE
0.0999 0.0349 0.3656 0.2634−0.5662 −0.7852 −0.4228 −0.6840−0.0012 −0.0724 −0.4965 −0.6403−0.1241 −0.4555 −0.1435 −0.5838

In [8], we proposed a technique to interpret the acquired clusters based on the relation between attributes of objects and central orders. We applied this method to clusters
derived by the EBC and TMSE methods. Table 2 shows Spearman’s ρ between central
orders of each cluster and an order of objects sorted according to the speciﬁc object
attributes. For example, the third row presents the ρ between the central order and thesorted object sequence according to their price. Based on these correlations, we wereable to learn what kind of object attributes affected the preferences of the respondents
in each cluster. We will comment next on each of the object attributes.
Almost the same observations were obtained by both EBC and TMSE. The attribute A1 shows whether the object tasted heavy (i.e., high in fat) or light (i.e., low
in fat). The positive correlation indicate a preference for heavy testing. The cluster 2
respondents preferred heavy-tasting sushi. The attribute A2 shows how frequently the
respondent eats the sushi. The positive correlation indicates a preference for the sushi
that the respondent infrequently eats. Respondents in both clusters preferred the sushi
they usually eat. No clear difference was observed between clusters. The attribute A3
is the prices of the objects. The positive correlation indicates a preference for economical sushi. The cluster 2 respondents preferred more expensive sushi. The attribute A4shows how frequently the objects are supplied at sushi shops. The positive correlation
indicates a preference for the objects that fewer shops supply. Though the correlationof cluster 1 was rather larger, the difference was not very clear. Roughly speaking, themembers of cluster 2 preferred more heavy-tasting and expensive sushi than those ofcluster 1.
In this paper, we propose a new technique based on the changes in object ranks.
First, a central order of all the sample orders was calculated, and was denoted by Ō∗ .

17
Table 3. The top 10 ranked up and the worst 10 ranked down sushi

#
1
2
3
4
5
6
7
8
9
10
91
92
93
94
95
96
97
98
99
100

Cluster 1
2313
egg ♣
cucumber roll ♣fermented bean roll ♣octopusdeep-fried tofu ♣salad ♣
pickled plum & perilla leaf roll ♣fermented bean ♣perilla leaf roll ♣raw beef..
.
ﬂying ﬁsh ♦young yellowtail ♦battera ♦sea bassamberjack ♦hardtail ♦ﬂuke ﬁnabalone ♥sea urchin ♠salmon roe

+74
+62
+38
+36
+33
+29
+28
+26
+24
+21
-10
-12
-13
-14
-37
-41
-46
-63
-84
-85

Cluster 22687
ark shell ♥crab liver ♠turban shell ♥sea bass
abalone ♥tsubu shellangler liver ♠sea urchin ♠clam ♥
hardtail ♦..
.
chili cod roe roll ♣pickled plum roll ♣shrimptuna roll ♣egg ♣salad roll ♣deep-fried tofu ♣salad ♣
octopus
squid

+63
+39
+26
+23
+22
+16
+16
+15
+13
+13
-15
-15
-17
-19
-19
-27
-30
-32
-57
-82

NOTE: Sushi in each cluster derived by k-o’means-EBC were sorted in descending order ofrankup(xj ) (Equation (18)). In top row labeled “#”, the sizes of clusters were listed. The upper
half of the tables show the ranked up sushi, and the bottom half show the ranked down sushi.
Just to the right of each sushi name, the rankup(xj ) values are shown.

Next, for each cluster, the central orders were also calculated, and were denoted by Ōk .
Then, for each object xj in X ∗ , the difference of ranks,
rankup(xj ) = r(Ō∗ , xj ) − r(Ōk , xj ),(18)

was derived. We say that xj is ranked up if rankup(xj ) is positive, and that it is ranked
down if rankup(xj ) is negative. If the object xj was ranked up, it was ranked higher
in cluster center Ōk than in the entire center Ō∗ . By observing the sushi whose theabsolute values of rankup(xj ) were large, we investigated the characteristics of eachcluster. Table 3 list the most 10 ranked up and the most 10 ranked down sushi in clusters
derived by k-o’means-EBC. That is to say, we show the objects whose rankup(xj ) werethe 1st to 10th largest, and were the 1st to 10th smallest. The upper half of the tables
shows the ranked up sushi, and the bottom half shows the ranked down sushi. In the top
row labeled “#”, the sizes of the clusters are listed. Sushi names that we were not able
to translate into English were written using their original Japanese names in italics. Just
to the right of each sushi name, the rankup(xj ) values are shown.

18

We interpreted this table qualitatively. In this table, the mark ♠ indicates objects
whose internal organs, such as liver or sweetbread, are eaten. The sushi marked by ♦are so-called blue ﬁsh, and those marked by ♥ are clams or shells. These sushi wererather substantial and oily, as revealed in the A1 row of Table 2. However, we could
not conclude that the respondents in cluster 2 preferred simply oily sushi. For example,
sushi categorized as a red ﬁsh meat, e.g., fatty tuna, were not listed in the table, because
the preference of sushi in this category were similar in both clusters. We can say thatthe respondents in cluster 2 preferred rather oily sushi, especially blue ﬁsh, clam/shell,or liver. The sushi marked by ♣ are very economical. Though these sushi were fairly
ranked up in cluster 1, this would not indicate a preference for economical sushi. These
would be ranked up because these respondents had sushi that they disliked more thanthese inexpensive types of sushi. Therefore, to interpret the acquired cluster of orders,
not only should the values of equation (18) be observed, but also the kind of objects
that were ranked up or ranked down.

6

Conclusions

We developed a new algorithm for clustering orders called the k-o’means-EBC method.
This algorithm is far more efﬁcient in computation and memory usage than k-o’meansTMSE. Therefore, this new algorithm can be applied even if the number of objects L∗ islarge. In the experiments on artiﬁcial data, our k-o’means outperformed the traditional
hierarchical clustering. For artiﬁcial data, the prediction ability of k-o’means-TMSE
is almost equal to that of k-o’means-EBC. Therefore, by taking computational cost
into account, it could be concluded that the k-o’means-EBC method was superior to
the k-o’means-TMSE for clustering orders. Additionally, we advocated the method tointerpret the acquired ordinal clusters.
We plan to improve this method in the following ways. During clustering orders,
undesired merges of clusters more frequently occur than in clustering of real value
vectors. To overcome this defect, it is necessary to improve the initial clusters. For
applying ordinal clustering to DNA microarray data, the curse of dimensionality must
be solved. We want to develop a dimension reduction technique for orders like PCA. In
the case of an Euclidean space, there are many points far from one point. However, in acase of a space of orders (a permutation group), the order most distant from one orderis unique, i.e., the reverse order. Therefore, there are biases for central orders to become
exact reversals of themselves. We also would like to lessen this bias.
Acknowledgments: This work is supported by the grants-in-aid 14658106 and 16700157
of the Japan society for the promotion of science.

References
1. Luaces, O., Bayón, G.F., Quevedo, J.R., Dı́ez, J., del Coz, J.J., Bahamonde, A.: Analyzing
sensory data using non-linear preference learning with feature subset selection. In: Proc. of
the 15th European Conf. on Machine Learning. (2004) 286–297 [LNAI 3201].
2. Fujibuchi, W., Kiseleva, L., Horton, P.: Searching for similar gene expression proﬁles acrossplatforms. In: Proc. of the 16th Int’l Conf. on Genome Informatics. (2005) P143

19
3. Everitt, B.S.: Cluster Analysis. third edn. Edward Arnold (1993)
4. Marden, J.I.: Analyzing and Modeling Rank Data. Volume 64 of Monographs on Statistics
and Applied Probability. Chapman & Hall (1995)
5. Branting, L.K., Broos, P.S.: Automated acquisition of user preference. Int’l Journal of
Human-Computer Studies 46 (1997) 55–77
6. Joachims, T.: Optimizing search engines using clickthrough data. In: Proc. of The 8th Int’l
Conf. on Knowledge Discovery and Data Mining. (2002) 133–142
7. Olson, C.F.: Parallel algorithms for hierarchical clustering. Parallel Computing 21 (1995)
1313–1325
8. Kamishima, T., Fujiki, J.: Clustering orders. In: Proc. of The 6th Int’l Conf. on Discovery
Science. (2003) 194–207 [LNAI 2843].
9. Kendall, M., Gibbons, J.D.: Rank Correlation Methods. ﬁfth edn. Oxford University Press(1990)
10. Dwork, C., Kumar, R., Naor, M., Sivakumar, D.: Rank aggregation methods for the Web. In:
Proc. of The 10th Int’l Conf. on World Wide Web. (2001) 613–622
11. Jain, A.K., Dubes, R.C.: Algorithms for Clustering Data. Prentice Hall (1988)
12. Fligner, M.A., Verducci, J.S.: Distance based ranking models. Journal of The Royal Statistical Society (B) 48(3) (1986) 359–369
13. Thurstone, L.L.: A law of comparative judgment. Psychological Review 34 (1927) 273–286
14. Mosteller, F.: Remarks on the method of paired comparisons: I — the least squares solution
assuming equal standard deviations and equal correlations. Psychometrika 16(1) (1951) 3–9
15. de Borda, J.C.: On elections by ballot (1784). In McLean, I., Urken, A.B., eds.: Classics ofSocial Choice. The University of Michigan Press (1995) 81–89
16. Mallows, C.L.: Non-null ranking models. I. Biometrika 44 (1957) 114–130
17. Arnold, B.C., Balakrishnan, N., Nagaraja, H.N.: A First Course in Order Statistics. John
Wiley & Sons, Inc. (1992)
18. Kamishima, T., Motoyoshi, F.: Learning from cluster examples. Machine Learning 53 (2003)
199–233
19. Kamishima, T.: Nantonac collaborative ﬁltering: Recommendation based on order responses.
In: Proc. of The 9th Int’l Conf. on Knowledge Discovery and Data Mining. (2003) 583–588