TD 1 - Clustering via Mode Seeking

The goal of this lab session is to implement the graph-based hill-climbing algorithm of Koontz et al., together with its variant ToMATo based on topological persistence. Both were presented during the lecture.

Graph-based hill-climbing

Class Point

• an instance attribute coords (a tuple), to store the coordinates of the point,
• a method sqrt_distance(point_a, point_b), which computes the squared Euclidean distance between two points,
• a method __str__ for pretty printing on the console... :-)

Class HillClimbing

• cloud, a list to store the data points,
• neighbors, a list to store the indices of the neighbors of each data point in the neighborhood graph (one list per point),
• density, a list to store the density values at the data points (one value per point),
• parent, a list to store the parent of each point in its cluster tree,
• label, a list to store the cluster label of each point.

• read_data(filename), to read the data points from a file and store it in cloud,
• compute_neighbors(k), to compute and store the k-NN graph of the data points. For this you can apply a naive strategy and compute all the pairwise distances (watch out that a point does not belong to its k-NN set).
• compute_density(k), to compute and store a density estimate at every data point. For this you can use the so-called distance to measure, defined as the average of the squared distances of the data point to its k nearest neighbors. The value of the density estimate is then the inverse of the square root of that quantity:
  
  \( f(p) = \frac{1}{\sqrt{\tfrac{1}{k}\sum\limits_{q \in k\text{-NN}(p)} \lvert p-q \rvert^2}} \)
  
• compute_forest(k), to connect each data point to its parent, which is its neighbor in the k-NN graph with highest density value among the ones with strictly higher density (if they exist) or itself (otherwise).
• compute_labels(), to extract the clusters' labels from the array of parents. Note that the labels do not have to be contiguous nor to start at 0. Hence, you can simply take the index of the root of the tree for each cluster.

Tests

cluster(cloud, labels, neighbors, k)

We recommend you test your code on the data set test first, to have fast runs and easy debugging. Then you can proceed with crater (takes a few dozens of seconds to handle on a decent machine), and then finally spirals (takes about half an hour). Try the following sets of parameters:

// parameters for test.xy
kDensity = 10;
kGraph = 5;

// parameters for crater.xy
kDensity = 50;
kGraph = 15;

// parameters for spirals.xy
kDensity = 100;
kGraph = 30;

ToMATo

• compute_persistence(k), which sorts the data points by decreasing estimated density values, then computes the 0-dimensional persistence of the superlevel sets of the density estimator via a union-find on the k-NN graph. At this stage you do not have to implement a full-fledged union-find, just use the trees stored in the field parent to do the merges. Whenever a cluster C is to be merged into another cluster C', simply hang the root of the tree of C to the root of the tree of C'.

• compute_persistence(k, tau), which merges only the clusters whose corresponding density peaks have prominence less than tau.

// parameters for test.xy
String filename = "test.xy";
kDensity = 10;
kGraph = 5;
tau = 0.35;

// parameters for crater.xy
String filename = "crater.xy";
kDensity = 50;
kGraph = 15;
tau = 2;

// parameters for spirals.xy
kDensity = 100;
kGraph = 30;
tau = 0.03;

Optimizations

1. Implement your own data structure for efficient proximity queries, kd-tree for instance. For this you can refer to Section 1 of the following TD from INF562 for building and queryig a kd-tree, and if you have more time you can even do the whole TD which is about Mean-Shift, the other mode-seeking algorithm mentioned in today's lecture. :-)
2. If your data points lie in 2D or 3D, then you can build their Delaunay triangulation. Then, the k-NN of each data point are found among its neighbors, their own neighbors, etc. within the triangulation, hence by a simple graph traversal. For this you can refer to the session on Delaunay triangulations in the course INF562.
3. You can also convert your code to C++ and use one of the existing libraries for NN search, such as ANN (which is based on a variant of the kd-trees and is extremely efficient in low and medium dimensions). This is probably the best option for optimizing this part of you code.