My Solutions
I will try to elaborate on what solutions I found, their implementations and those implementation's running time.
As I don't have any real computer science background, please point out any mistakes/unclarities you find.
First: The problem, in general terms, is to partition a set of points from
$A \subset \mathbb{R}^d$ into groups such that every point in them is $\epsilon$-reachable
from every other point in the group.
Reachable,for $x_1, x_n \in A$, means there is a sequence of points $x_1,\dots,x_n \in A$ s.t. $d(x_i,x_{i+1})<\epsilon$.
In other terms: If two points are closer than $\epsilon$ to each other, the algorithm must assign them to the same group/cluster.
DBSCAN:
DBSCAN is a super general algorithm. I will use $B_{\epsilon}(x)$ for the $d$-dimensional sphere around $x$ with
respect to some arbitrary metric $d(\cdot,\cdot)$. DBSCAN finds clusters by defining:
- $x \in A$ is a core point, if at least $N$ other points of $A$ lie in $B_{\epsilon}(x)$
- $x \in A$ is directly reachable from $y \in A$ if both are core points and $d(x,y)<\epsilon$
- $x \in A$ is reachable from $y \in A$ if there exists a sequence $x_1 =x, \dots, x_n =y \in A$ s.t. $x_{i+1}$ is reachable from $x_i$
- non reachable points are outliers
Now, as you might see, for $N=2$ (if you include every point itself when counting) and $d(\cdot, \cdot)$ the usual
euclidean metric
this defines exactly the problem I stated above. DBSCAN clusters together core-points deterministically,
but also adds non-core points to clusters. So if $x \in A$ is reachable from core point $y$, it might
end up in the same cluster without being a core point itself. I say might, because non-core points
may be reachable from different clusters. The assignation is then determined by the order in which
the algorithm runs through the clusters. This ambiguity however resolves with $N=2$, as then every point
that is reachable from another is automatically a core point, that is reachability as stated
above is then a symmetric relation. Clusters are then unambiguous.
In case of an indexing structure that does neighborhood querys in $O(logn)$, the runtime complexity is
apparently $O(n \; logn)$, but don't take my word for it.
I used scikit-learn's implementation of DBSCAN (http://scikit-learn.org/stable/modules/generated/sklearn.cluster.DBSCAN.html).
I failed to find out what runtime complexity this implementation has, but the docs state memory complexity
of $O(n.d)$.
Intersection graph:
As it has been proposed in a different answer: One can also think of the problem as a graph-theoretical one.
The graphs vertices are represented by the points in $A$, with two vertices being adjacent iff their
points are in $\epsilon$-distance of each other. On this graph the problem translates to finding the
connected components.
The paper of
Bentley, Jon L.; Stanat, Donald F.; Williams, E. Hollins, Jr. (1977), "The complexity of finding fixed-radius near neighbors", Information Processing Letters, 6 (6): 209–212
provides an algorithm for constructing this graph that runs in linear time and is essentially identical to the one proposed
by @D.W.. Finding the connected components can also be done in linear time, e.g. by breadth- or depth-first algortihms.
Python Code:
For the Code, I mainly used the DBSCAN implementation of scikit-learn, networkx for graphs and finding
connected components (it is open source but I couldnt find any reference on the runtime complexity)
and scipy.spatial's cKDTree implementation (which I do not understand in detail).
from sklearn.cluster import DBSCAN
from sklearn.datasets.samples_generator import make_blobs
import networkx as nx
import scipy.spatial as sp
def cluster(data, epsilon,N): #DBSCAN, euclidean distance
db = DBSCAN(eps=epsilon, min_samples=N).fit(data)
labels = db.labels_ #labels of the found clusters
n_clusters = len(set(labels)) - (1 if -1 in labels else 0) #number of clusters
clusters = [data[labels == i] for i in range(n_clusters)] #list of clusters
return clusters, n_clusters
Is my minimal implementation of the previously described goal. The method returns clusters and their size.
For the graph method, I wrote a little class. The class describes a graph whose vertices are the indices
of the corresponding points in $A$ as list. Using the cKDTree's data structure, finding pairs of
reachable points is lightning fast. You can then slice out the elements via the indices.
class IGraph:
def __init__(self, nodelst=[], radius = 1):
self.igraph = nx.Graph()
self.radii = radius
self.nodelst = nodelst #nodelst is array of coordinate tuples, graph contains indices as nodes
self.__make_edges__()
def __make_edges__(self):
self.igraph.add_edges_from( sp.cKDTree(self.nodelst).query_pairs(r=self.radii) )
def get_conn_comp(self):
ind = [list(x) for x in nx.connected_components(self.igraph) if len(x)>1]
return [self.nodelst[indlist] for indlist in ind]
def graph_cluster(data, epsilon):
graph = IGraph(nodelst = data, radius = epsilon)
clusters = graph.get_conn_comp()
return clusters, len(clusters)
I also did a little benchmark, but it is in no way representative. I create a cluster with:
centers = [[1, 1,1], [-1, -1,1], [1, -1,1]]
X,_ = make_blobs(n_samples=N, centers=centers, cluster_std=0.4,
random_state=0)
For different N, the methods compare differently, on my machine for example:
- $N = 20000$: "cluster": 293.711 ms, "graph_cluster": 422.089 ms
- $N = 7000$: "cluster": 87.455 ms, "graph_cluster": 95.322 ms
- $N = 700$: "cluster": 3.764 ms, "graph_cluster": 1.689 ms
For small datasets, and this specific clustering, the graph solution wins. But apparently
the DBSCAN shows better scaling behavior. I haven't looked at different dimensions of the data,
but that may also be interesting. Also: the results are exactly the same for both algortihms.
Additional Note: The graph solution scales heavily with $\epsilon$, which is 0.1 in my code. Larger values make it far worse! In other words: If you have small clusters, the graph solution outperforms DBSCAN by ages, for large clusters it is the other way around (apparently).
