We present an algorithm for the 2-clustering problem with cluster size constraints in the plane assuming `1-norm, that works in O(n3 log n) time and O(n) space. Such a procedure also solves a full version of the problem, computing the optimal solutions for all possible constraints on cluster sizes. The algorithm is based on a separation result concerning the clusters of any optimal solution of the problem and on an extended version of red-black trees to maintain a bipartition of a set of points in the plane.
Clustering is one of the most used techniques in statistical data analysis and
machine learning [
A central problem in distance clustering is the well-known Minimum
Sum-ofSquares Clustering (MSSC), also called Variance-based Clustering [
2 ||x − CAi ||2 CAi = argmin X ||μ − x||22 = μ∈Rd x∈Ai
1
|Ai| x∈Ai
MSSC is difficult in general; indeed for an arbitrary dimension d (given in the
instance) the problem is NP-hard even if the number of clusters is fixed to k = 2
[
Moreover, a known variant of MSSC, called k-Medians Problem, is given by
substituting the squared-euclidean norm ||·||22 with the `1-norm ||·||1 [
In many applications, often people have some information on the clusters [
In this work we study the 2-clustering problem in the plane, with cluster size constraints, assuming `1-norm. The relevance of the 2-clustering problems is also due to the wide spread of hierarchical clustering techniques [11, Sec.14.3.12], that repeatedly apply the 2-clustering as the key step. This leads to natural size-constrained versions of so-called divisive hierarchical clustering.
We recall that the 2-clustering problem with cluster size constraints has been
studied in [
Here we prove that, assuming dimension 2 and `1-norm, the 2-clustering problem with cluster size constraints can be solved in O(n3 log n) time and O(n) space. Our procedure actually solves a full version of the problem, computing the optimal solutions for all possible sizes of cluster. Clearly this yields the solution to the general unconstrained 2-clustering problem in the plane with `1-norm.
The algorithm is based on a separation result stating that, in our hypotheses, the clusters of any optimal solution are separated by curves of some special forms, which can be analysed in a suitable linear order. The algorithm also make use of a particular extension of red-black trees to maintain a bipartition of a set of points in the plane.
We remark that in this work we propose an efficient method for obtaining a
solution that is globally optimal, instead of a locally optimal solution as yielded
by various heuristics [
In this section we introduce some basic notions of distance-based clustering problems on the plane assuming the Manhattan norm (also called `1-norm), which is defined by kxk1 = |x1| + |x2| for any x = (x1, x2) ∈ R2.
Given a set X ⊂ R2 of n points in the plane, a k-partition of X is a family k {A1, A2, ..., Ak} of k nonempty subsets of X such that ∪i=1Ai = X and Ai∩Aj = ∅, for i 6= j. Each Ai is called cluster. The centroid CA of a cluster A ⊂ X is CA = argmin X kx − μk1 = argmin X (|x1 − μ1| + |x2 − μ2|)
μ∈R2 x∈A μ∈R2 x∈A Since the objective function in the minimization is not strictly convex, the centroid CA is not necessarily unique. Indeed it is well-known that the centroid CA is the component-wise median, i.e. the abscissa (ordinate) of the centroid is the median of the abscissae (ordinates) of the points in A. The weight W (A) of a cluster A is
W (A) = X
kx − CAk1 x∈A while the weight of a k-partition {A1, A2, ..., Ak} (also called k-clustering ) is k W (A1, A2, · · · , Ak) = X W (Ai).
1 The classical Clustering Problem in the plane assuming the Manhattan norm is stated as follows.
Definition 1 (Clustering Problem in R2 under `1-norm). Given a point set X ⊂ R2 of cardinality n and an integer k, 1 < k < n, find a k-clustering {A1, A2, ..., Ak} that minimizes the weight W (A1, A2, · · · , Ak) = P1k W (Ai). Note that here k is included in the instance of the problem. If k is fixed the problem is called k-Clustering in R2. Our main contribution in this paper concerns the 2-Clustering in R2 under `1-norm.
We are interested in a version of clustering problem where the cluster sizes are constrained. Formally, the problem can be stated as follows: Definition 2 (Size Constrained Clustering Problem in R2 under `1norm). Given a point set X ⊂ R2 of cardinality n, an integer k > 1 and k positive integers m1, m2, ..., mk such that P1k mi = n, find a k-clustering {A1, A2, ..., Ak} with |Ai| = mi for i = 1, ..., k, that minimizes the weight W (A1, A2, · · · , Ak) = P1k W (Ai).
We denote this problem by SCC-2 (under `1-norm). We stress that in the SCC-2 problem the integers n, k, m1, . . . , mk are part of the instance. On the contrary, if k is fixed and does not belong to the instance, the problem is denoted by k-SCC-2.
In the following sections we present an algorithm for the 2-Clustering problem in R2 under `1-norm, which solves at the same time the full version of the 2SCC-2 problem, i.e. that yields the solutions to all 2-SCC-2 problems for every value of m1 ∈ {1, 2, . . . , bn/2c} and the same point set X = {x1, x2, ..., xn} ⊂ R2 in input. The procedures works in O(n3 log n) time and is based on a separation result concerning the optimal solutions of 2-Clustering in the plane, presented in Section 3, and on an augmented version of red-black trees used to maintain bipartitions of a set of real points, described in Section 4. 3
In this section we present a separation result concerning the optimal solution of
2-SCC-2 problem under `1-norm. It turns out that the clusters of the optimal
solutions of this problem are separated by plane curves of 8 different types, we
will call Ci and Si, i = 1, 2, 3, 4, respectively. Similar results are obtained in [
First, we show that in an optimal 2-clustering the swapping of two elements in different clusters does not decrease the solution value.
Lemma 1 (Swapping Lemma). Let {A, B} be an optimal solution of the 2SCC-2 problem. Then, for any a ∈ A and b ∈ B, it holds
||a − CA||1 + ||b − CB||1 ≤ ||a − CB||1 + ||b − CA||1 .
Proof. For the sake of simplicity, we omit the subscript 1 in denoting norm. By contradiction, let us assume that for some a ∈ A and b ∈ B
||a − CA|| + ||b − CB|| > ||a − CB|| + ||b − CA|| .
W (A, B) = X ||x − CA|| + X ||x − CB|| = > = x∈A
X x∈Ar{a}
X x∈Ar{a}∪{b}
x∈B ||x − CA|| +
X x∈Br{b}
X x∈Br{b}∪{a}
||x − CB|| . ||x − CA|| + ||b − CA|| + ||x − CB|| + ||a − CB|| = Now, it is clear that the 2-clustering {A0, B0}, with A0 = A r {a} ∪ {b} and B0 = B r {b} ∪ {a}, is a feasible solution since |A0| = |A| and |B0| = |B|. Furthermore, by definition of centroid it holds: X ||x − CA||+ X ||x − CB|| ≥ ||x − CA0 ||+
X ||x − CB0 || = W (A0, B0) x∈A0 x∈B0 x∈B0 and hence W (A, B) > W (A0, B0). However, this is in contradiction since {A, B} is an optimal solution. tu X x∈A0 We are now able to obtain the general form of the equation for the curves separating the clusters in an optimal 2-clustering.
Theorem 1 (Separation Result). In an optimal solution {A, B} of the 2SCC-2 problem, the clusters A and B are separated by a curve of equation ||z − α||1 − ||z − β||1 = g (1) where the variable z ranges over R2, while α, β ∈ R2 and g ∈ R are suitable constants.
Proof. Let CA and CB be the centroids of A and B respectively. Then, by Lemma 1 the inequality
||a − CA||1 − ||a − CB||1 ≤ ||b − CA||1 − ||b − CB||1
is satisfied for all a ∈ A and b ∈ B, and hence it also holds by taking the
maximum over a’s on the left and the minimum over b’s on the right:
cA := ma∈aAx{||a − CA||1 − ||a − CB||1} ≤ mb∈iBn{||b − CA||1 − ||b − CB||1} =: cB
Therefore, choosing g ∈ [cA, cB], all the points of A and B are contained in the
region defined respectively by ||z − CA||1 − ||z − CB||1 ≤ g and ||z − CA||1 −
||z − CB||1 ≥ g. The common boundary of the two regions has equation
||z − CA||1 − ||z − CB||1 = g, and thus the proof is concluded by setting α = CA
and β = CB. tu
Remark 1. It is noteworthy to observe that both the Swapping Lemma and the
Separation Result can be extended to the case of `p-norm with any p ≥ 1 and
dimension d ≥ 2, provided the norm is raised to the p-th power [
In this section we describe a natural data structure to maintain a bipartition of a set of real numbers, with possible multiple values, able to test membership, to execute insertion and deletion, and to perform two further operations: selecting the i-th smallest element in either of the two clusters of the bipartition, for any integer i, and moving an element from one cluster to the other.
Let X be a finite multiset of real numbers, i.e. a finite set X of reals where any value may occur several times, represented as a finite ordered sequence
X = (x1, x2, . . . , xn), xi ∈ R for every i, and x1 ≤ x2 ≤ · · · ≤ xn. Moreover, consider a bipartition {A, B} of X, i.e. two non-empty subsets A ⊆ X, B ⊆ X such that A ∩ B = ∅ and A ∪ B = X. In the following A and B are called clusters. Operations Member, Insert and Delete can be defined in order to specify the cluster the elements belong to. More precisely, for every z ∈ R, Member(z) = (j1, j2), where j1 (resp. j2) is the number of xi’s in A (resp. in B) such that xi = z, Insert(z, A) adds a new element z to A, Insert(z, B) does the same to B, while Delete(z, A) and Delete(z, B) execute the delete operations. S1 C2 C3 S3 f d d f α α
C1
C2 C3 f d d f
S2
C1 α Moreover, for every i ∈ {1, 2, . . . , n}, we define
SelectA(i) = SelectB(i) = a if i ≤ |A| and a is the i-th smallest element of A ⊥ if i > |A| b if i ≤ |B| and b is the i-th smallest element of B .
⊥ if i > |B| Finally, for any z ∈ R, if A contains an element xi = z then MoveAB(z) deletes a value z from A and adds it to B, otherwise MoveAB(z) = ⊥. The operation MoveBA(z) is defined symmetrically.
A natural data structure, able to executes the above operations in time
logarithmic with respect to n = |X|, is an augmented version of the well-known
red-black trees [
For each internal node v, we denote by key (v), left (v) and right (v) the key contained in v, the left son and the right son of v, respectively. Clearly, v can be implemented as a record including the values mA(v), mB(v)), cA(v), cB(v).
Procedures Member, Insert and Delete can be defined as in the standard
redblacks trees [
Operations Select and Move can be executed as follows. To compute SelectA(i) one checks whether i ≤ cA(r), where r is the root of T : in the affirmative case a standard searching procedure (we avoid to define for lack of space) is called, otherwise the value ⊥ is returned. SelectB(i) is defined in a similar way.
To compute MoveAB(z) one first determines the node v containing z by calling a rather usual procedure, then the commands mA(v) = mA(v) − 1 and mB(v) = mB(v) + 1 are executed. Procedure MoveBA(z) is defined symmetrically.
It is clear that all previous procedures can be executed in O(log n) time. 5
Our main result, presented in the next section, is based on a procedure that updates a current bipartition by moving one point from a cluster to the other. Here we want to show how to update the centroids and the weights of the clusters of a current bipartition after executing such a moving operation.
To this end we first deal with the problem of updating centroid and weight of a set of real numbers subject to insert and delete operations. The results will be easily extended to a set of points in R2.
Let A be a cluster in R, i.e. a finite sequence of ordered real numbers, A = {a1, a2, . . . , ak}, where a1 ≤ a2 ≤ · · · ≤ ak. Assuming `1-norm it is well-known that the centroid of A is the median of the set, that is the value
M (A) = ( a k+1 if k is odd
2 1 a k + a k2 +1 if k is even 2 2 Thus, if A is a cluster of a bipartition (of a finite multiset X ⊂ R of n elements), maintained by a Bipartition Red-Black Tree, then M (A) can be computed in O(log n) time by calling procedure SelectA(j) for suitable j’s.
Now, let us recall the definition of the left and right part of the weight of A; given m = k +21 , we have
X L(A) := ai ,
R(A) :=
X m<i≤k ai It is not difficult to see that the weight of A, given by W (A) = Pk i=1 |ai −M (A)|, can be computed as the difference between R(A) and L(A).
Proposition 1. For every ordered sequence of real numbers A = {a1, a2, . . . , ak}, where a1 ≤ a2 ≤ · · · ≤ ak, we have
W (A) = R(A) − L(A) Proof. If k is even then m = k+1 is not integer and hence
2 k W (A) = X |ai − M (A)| = i=1
X (M (A) − ai) +
X (ai − M (A)) 1≤i<m = −L(A) + R(A) m<i≤k
As a consequence, in order to maintain W (A) it is sufficient to update L(A) and R(A). This can be done by a straightforward procedure that implements, as a chain of if-then-else instructions, the rules for computing the left and right part of the weight of a current cluster upon insertion and deletion. For lack of space here we omit their cumbersome formal description. These rules together with tu Proposition 1, define procedures UpdateW-ins(A, x) and UpdateW-del(A, aj ), which update the weight of a current cluster A ⊂ R, respectively upon insertion and deletion of an element of value x ∈ R.
They can be executed by using a Bipartition Red-Black tree with the corresponding Select operation, introduced in the previous section. Note that, for updating both L(A) and R(A), the Select procedure is necessary to find in the tree the i-th value of A when required. As a consequence also UpdateW-ins(A, x) and UpdateW-del(A, aj ) can be executed in logarithmic time. 6
In this section we apply the previous results, in particular the separation property of Section 3, to define an algorithm for the 2-SCC-2 problem in the full version, i.e. for all possible sizes of cluster. Thus the same algorithm also yields a solution of the general 2-Clustering problem in the plane.
By our Theorem 1, the optimal solution of a 2-SCC-2 problem for a set X ⊂ R2 of n points and a cluster size m, consists of a bipartition (A, B), with |A| = m and |B| = n − m, where A and B are separated by a curve of the form Ci or Si, i = 1, 2, 3, 4.
Note that each curve of type Ci can be obtained from a curve of any other form Cj , j 6= i, by means of a simple rotation of the plane. Similarly each curve of type Si can be obtained from another curve of any different type Sj by means of a rotation and/or a reflection through a line.
For these reasons we treat in detail only the case of bipartitions separated by a curve of type C1. The cases of curves of type C2, C3 and C4 are solved by the same algorithm through a rotation of the input points. The same algorithm can be easily adapted to the cases of separating curves of type Si, i = 1, 2, 3, 4.
Let X = {p1, p2, . . . , pn} ⊂ R2 be a set of n points and let us set pi = (xi, yi), for any i = 1, 2 . . . , n. Without loss of generality we may assume that all pi’s lie in the first quadrant, i.e. xi > 0 and yi > 0 for all i.
We define a procedure that searches an optimal bipartition of X whose clusters are separated by a curve of type C1. The procedure first sorts X according with 3 parameters: the two coordinates and their sum, i.e. the values xi’s, yi’s and xi + yi’s, respectively. Then the elements of X are processed by three for-loops, each nested into the other, corresponding to the three parameters considered above. The algorithm evaluates all possible bipartitions of X whose clusters are separated by curves of type C1 and, for every admissible cluster size, only the bipartition with minimum weight is maintained. These bipartitions are considered in a linear order and, at each step, a new bipartition of X is obtained from the current one {A, B} by swapping a point from cluster B to cluster A.
Such a current bipartition {A, B} is maintained by two bipartition red-black trees: one for the 2-clustering {Ax, Bx} of the set {xi | pi = (xi, yi) ∈ X}, where Ax = {xi | pi ∈ A} and Bx = {xi | pi ∈ B}, and the other for the 2-clustering (Ay, By) of {yi | pi = (xi, yi) ∈ X}, where Ay = {yi | pi ∈ A} and By = {yi | pi ∈ B}.
An Update subroutine can be defined which, for an input (A, B, pi), where {A, B} is the current bipartition of X and pi ∈ B, computes the new bipartition moving pi from B to A and returns the weights of its clusters. The subroutine applies the procedures described in Sections 4 and 5 and it is defined by the following scheme: Procedure Update(A, B, pi) begin
MoveBxAx (xi) MoveByAy (yi) w(Ax) := UpdateW-ins(Ax, xi) w(Bx) := UpdateW-del(Bx, xi) w(Ay) := UpdateW-ins(Ay, yi) w(By) := UpdateW-del(By, yi) return (w(Ax) + w(Ay), w(Bx) + w(By)) end
By the discussion presented in Sections 4 and 5 also this subroutine only requires O(log n) time.
Thus, the overall algorithm is described in detail by the procedure given below, where three bipartitions are actually maintained, denoted by (A, B), (A0, B0) and (A00, B00), respectively. Each of them need two Bipartition Red-Black trees, one for the abscissae and the other for the ordinates.
Moreover, instruction Sortx (resp., Sorty and Sortx+y) yields an ordered list of the indices of the points in X, sorted in non-decreasing order with respect to the values xi (resp., yi and xi + yi).
Taking into account the previous discussion, it is easy to see that the algorithm works in O(n3 log n) time; moreover, it has a space complexity O(n) since each Bipartition Red-Black tree only requires linear space.
begin (i1, i2, . . . , in) := Sortx(1, 2, . . . , n) (j1, j2, . . . , jn) := Sorty(1, 2, . . . , n) (k1, k2, . . . , kn) := Sortx+y(1, 2, . . . , n) A := ∅ ; w(A) := 0 B := {p1, p2, . . . , pn} compute the weight w(B) of cluster B i0 := 0 ; j0 := 0 ; x0 := 0 ; y0 := 0 for r = 0, 1, . . . , n do
A0 := A ; B0 := B ; w(A0) := w(A) ; w(B0) := w(B) for s = 0, 1, . . . , n do
A00 := A0 ; B00 := B0 ; w(A00) := w(A0) ; w(B00) := w(B0) for k = k1, k2, . . . , kn do if xk ≥ xir ∧ yk ≥ yjs then (w(A00), w(B00)) = Update(A00, B00, pk) h := min(|A00|, |B00|) z := w(A00) + w(B00) if W [h] > z then
W [h] := z Π[h] := (ir, js, k) if s 6= n ∧ xjs+1 ≥ xir then (w(A0), w(B0)) = Update(A0, B0, pjs+1 ) h := min(|A0|, |B0|) z := w(A0) + w(B0) if W [h] > z then if r 6= n then (w(A), w(B)) = Update(A, B, pir+1 ) h := min(|A|, |B|) z := w(A) + w(B)
W [h] := z Π[h] := (ir, js+1, 0) if W [h] > z then
W [h] := z Π[h] := (ir+1, 0, 0) return (W, Π) end 7
In this paper we have shown a O(n3 log n) time algorithm for the size-constrained 2-Clustering problem in the plane with `1-norm. It is clear that our results strictly depend on those conditions and it is not evident (at least at the present time) whether they hold under different hypotheses or for more general problems. Moreover, being our algorithm exact and working in polynomial time, we do not think it is necessary to produce here, in our hypotheses (2-clustering in R2 with `1-norm), an experimental comparison with traditional heuristics like k-Means, which instead yields an approximate solution, tackles a more general problem with larger number of clusters, and works in exponential time in the worst case.
However, we think that the present contribution could be the starting point
for a further analysis, still based on the separation result introduced in [
We also remark we have already obtained recently some results under other hypotheses. In particular, in a companion paper we present more efficient algorithms for solving exactly the same problem with Euclidean norm, based on rather different techniques of computational geometry. More precisely, it turns out that the 2-Clustering problem in R2 under Euclidean norm, with size con√ straint k, can be solved by an exact algorithm in O(n 3 k log2 n) time, while the full version of the same problem (i.e. for all k = 1, 2, . . . , bn/2c) is solvable in O(n2 log n) time.