k-means|| a.k.a. Scalable K-Means++

Bahman Bahmani et al. introduced k-means||, which is a faster version of k-means++.

This algorithm is taken from page 4 of their paper, Bahmani, B., Moseley, B., Vattani, A., Kumar, R., & Vassilvitskii, S. (2012). Scalable k-means++. Proceedings of the VLDB Endowment, 5(7), 622-633.

Unfortunately I don't understand those fancy Greek letters, so I need some help to understand how this works. As far as I understand this algorithm is an improved version of k-means++, and it uses oversampling, to reduce the number of iterations: k-means++ has to iterate $k$ times, where $k$ is the number of desired clusters.

I got a very good explanation through a concrete example of how k-means++ works, so I will use the same example again.

Example

I have the following dataset:

(7,1), (3,4), (1,5), (5,8), (1,3), (7,8), (8,2), (5,9), (8,0)

$k = 3$ (number of desired clusters)

$\ell = 2$ (oversampling factor)

I started to calculate it, but I'm not sure if I got it right, and have no idea about Steps 2, 4, or 5.

• Step 1: $\mathcal{C} \leftarrow$ sample a point uniformly at random from $X$

Let's say the first centroid is $(8,0)$ (same as in k-means++)

• Step 2: $\psi \leftarrow \phi_X(\mathcal{C})$

no idea

• Step 3:

• $d^2(x, \mathcal{C}) = [2, 41, 74, 73, 58, 65, 4, 90]$

We calculate the squared distances to the closest center to each point. In this case we have only one center so far, $(8,0)$.

• $\ell \cdot d^2(x, \mathcal{C}) = [4, 81, 148, 146, 116, 130, 8, 180]$

(Because $\ell = 2$ in this case.)

• $\text{cumulative } \ell \cdot d^2(x, \mathcal{C}) = [4, 85, 233, 379, 495, 625, 633, 813]$

Pick $\ell = 2$ random numbers in the interval $[0, 813)$. Say you pick $246.90$ and $659.42$. They fall in the ranges $[379, 495)$ and $[633, 813)$ which correspond to the 4th and 8th items respectively.

• Repeat it $\mathcal{O}(\log \psi)$ times, but what is $\psi$ (calculated in Step 2) in this case?

• Step 4: For $x \in \mathcal{C}$, set $w_x$ to be the number of points in $X$ closer to $x$ than any other point in $\mathcal{C}$.
• Step 5: Recluster the weighted points in $\mathcal{C}$ into $k$ clusters.

Any help in general or in this particular example would be great.

data points: (7,1),(3,4),(1,5),(5,8),(1,3),(7,8),(8,2),(5,9),(8,0)

l = 2 // oversampling factor

k = 3 // no. of desired clusters

Step 1:

Suppose the first centroid is $\mathcal{C}$ is $\{ c_1\} = \{ (8,0) \}$. $X = \{x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8\}=\{(7,1),(3,4),(1,5),(5,8),(1,3),(7,8),(8,2),(5,9)\}$

Step 2:

$\phi_X(\mathcal{C})$ is the sum of all smallest 2-norm distances (euclidean distance) from all points from the set $X$ to all points from $\mathcal{C}$. In other words, for each point in $X$ find the distance to the closest point in $\mathcal{C}$, in the end compute the sum of all those minimal distances, one for each point in $X$.

Denote with $d^2_{\mathcal{C}}(x_i)$ as the distance from $x_i$ to the closest point in $\mathcal{C}$. We then have $\psi = \sum_{i=1}^{n}d^2_{\mathcal{C}}(x_i)$.

At step 2, $\mathcal{C}$ contains a single element (see step 1), and $X$ is the set of all elements. Thus in this step the $d^2_{\mathcal{C}}(x_i)$ is simply the distance between the point in $\mathcal{C}$ and $x_i$. Thus $\phi = \sum_{i=1}^{n}{||x_i-c||^2}$.

$\psi = \sum_{i=1}^nd^2(x_i,c_1) = 1.41+6.4+8.6+8.54+7.61+8.06+2+9.4 = 52.128$ $log(\psi) = log(52.128) = 3.95 = 4 (rounded)$

Note however that in step 3, the general formula is applied since $\mathcal{C}$ will contain more than one point.

Step 3:

The for loop is executed for $log(\psi)$ previously computed.

The drawings are not like you understood. The drawings are independent, which means you will execute a draw for each point in $X$. So, for each point in $X$, denoted as $x_i$, compute a probability from $p_x = l d^2(x,\mathcal{C})/\phi_X(\mathcal{C})$. Here you have $l$ a factor given as parameter, $d^2(x,\mathcal{C})$ is the distance to the closest center, and $\phi_X(\mathcal{C})$ is explained at step 2.

The algorithm is simply:

• iterate in $X$ to find all $x_i$
• for each $x_i$ compute $p_{x_i}$
• generate an uniform number in $[0, 1]$, if is smaller than $p_{x_i}$ select it to form $\mathcal{C'}$
• after you done all draws include selected points from $\mathcal{C'}$ into $\mathcal{C}$

Note that at each step 3 executed in iteration (line 3 of the original algorithm) you expect to select $l$ points from $X$ (this is easily shown writing directly the formula for expectation).

for(int i=0; i<4; i++) {

// compute d2 for each x_i
int[] psi = new int[X.size()];
for(int i=0; i<X.size(); i++) {
double min = Double.POSITIVE_INFINITY;
for(int j=0; j<C.size(); j++) {
if(min>d2(x[i],c[j])) min = norm2(x[i],c[j]);
}
psi[i]=min;
}

// compute psi
double phi_c = 0;
for(int i=0; i<X.size(); i++) phi_c += psi[i];

// do the drawings
for(int i=0; i<X.size(); i++) {
double p_x = l*psi[i]/phi;
if(p_x >= Random.nextDouble()) {
X.remove(x[i]);
}
}
}
// in the end we have C with all centroid candidates
return C;


Step 4:

A simple algorithm for that is to create a vector $w$ of size equals to the number of elements in $\mathcal{C}$, and initialize all its values with $0$. Now iterate in $X$ (elements not selected in as centroids), and for each $x_i \in X$, find the index $j$ of the closest centroid (element from $\mathcal{C}$) and increment $w[j]$ with $1$. In the end you will have the vector $w$ computed properly.

double[] w = new double[C.size()]; // by default all are zero
for(int i=0; i<X.size(); i++) {
double min = norm2(X[i], C[0]);
double index = 0;
for(int j=1; j<C.size(); j++) {
if(min>norm2(X[i],C[j])) {
min = norm2(X[i],C[j]);
index = j;
}
}
// we found the minimum index, so we increment corresp. weight
w[index]++;
}


Step 5:

Considering the weights $w$ computed at the previous step, you follow kmeans++ algorithm to select only $k$ points as starting centroids. Thus, you will execute $k$ for loops, at each loop selecting a single element, drawn randomly with probability for each element being $p(i) = w(i)/\sum_{j=1}^m{w_j}$. At each step you select one element, and remove it from candidates, also removing its corresponding weight.

for(int k=0; k<K; k++) {
// select one centroid from candidates, randomly,
// weighted by w
// see kmeans++ and you first idea (which is wrong for step 3)
...
}


All the previous steps continues, as in the case of kmeans++, with the normal flow of the clustering algorithm

I hope is clearer now.

[Later, later edit]

I found also a presentation made by authors, where you can not clearly that at each iteration multiple points might be selected. The presentation is here.

[Later edit @pera's issue]

It is obvious that $log(\psi)$ depends on data and the issue you raised would be a real problem if the algorithm would be executed on a single host/machine/computer. However you have to note that this variant of kmeans clustering is dedicated to large problems, and for running on distributed systems. Even more, the authors, in the following paragraphs above the algorithm description state the following:

Notice that the size of $C$ is significantly smaller than the input size; the reclustering can therefore be done quickly. For instance, in MapReduce, since the number of centers is small they can all be assigned to a single machine and any provable approximation algorithm (such as k-means++) can be used to cluster the points to obtain k centers. A MapReduce implementation of Algorithm 2 is discussed in Section 3.5. While our algorithm is very simple and lends itself to a natural parallel implementation (in $log(\psi)$ rounds ), the challenging part is to show that it has provable guarantees.

Another thing to note is the following note on the same page which states:

In practice, our experimental results in Section 5 show that only a few rounds are enough to reach a good solution.

Which means you could run the algorithm not for $log(\psi)$ times, but for a given constant time.

• could you please extend your answer with the calculation for my example? – user1930254 Feb 4 '15 at 14:29
• I am a programmer, I think I can write that in code faster than typing here :). Hope it explains the algo. – rapaio Feb 4 '15 at 15:16
• Can you please explain what is the idea with log(Ksi) number of iterations? I don't understand the idea beneath it, it seems that the number of iterations will depend on the range of values of objects, which doesn't seem reasonable. For example, if objects have attribute values of about 1000, that could, for example, result that the error is about 1000, which means that there will be 3 iterations. On the other hand, if values are in range of 10 that could result that the error is of about 10 which results in 1 iteration. Shouldn't the number of iterations depend on number of objects? – Marko Apr 15 '15 at 10:16
• @pera I update the answer to clarify the issue you raised – rapaio Apr 21 '15 at 6:10
• @rapaio Thank you for your answer, I'm already going for the solution that will determine the number of iterations based on the number of medoids. Where x can be increased to get better initialization at the cost of a couple of iterations more. Do you fine that ok, based on the second part that you have given? Thanks, again. – Marko Apr 21 '15 at 7:41