First, we can derive an undirected graph where vertices are the vectors of v and edges are the pairs of points that are no farther apart than r. This doesn't require the distances to be a proper metric: it doesn't have to satisfy the triangle inequality, only that it is symmetrical (dist(a,b) == dist(b,a)) and that it is 0 iff two points are equal.
After that step, we can forget about the distances altogether and focus entirely on partitioning the graph. As others have noted, that partition is a Vertex Cover problem. However, it has a twist: we require all the vertices in the cover to be disjoint (i.e.: no vectors in v1 can be within distance r of each other).
To obtain the graph, we can use the efficient KDTree to compute just once the nearest-neighbor structure. In particular, we'll use kd.sparse_distance_matrix(kd, r) to obtain all pairs of points separated by r or less:
from scipy.spatial import KDTree
def get_graph(v, r):
kd = KDTree(v)
sm = kd.sparse_distance_matrix(kd, r)
graph = {i: {i} for i, _ in sm.keys()}
for i, j in sm.keys():
graph[i] |= {j}
return graph
(note: .sparse_distance_matrix() has a parameter p if we want other distances than Euclidean).
The partitioning of the graph can be done in various ways. @DanR shows a greedy approach that is very fast, but often suboptimal in the size of v1 that it finds.
The following shows instead a brute-force approach that is guaranteed to find a minimal solution if there is one. It is adequate for relatively small number of vectors and can provide a ground-truth for the optimal solution, when researching other heuristics.
To speed up the combinatorial search, we first observe that often some points are not within distance r of any other (i.e. the graph found above doesn't represent all the points of v). We can leave those other points aside, since they will necessarily be part of v1 but don't interfere with the rest of the search (below, we call them "singletons"). Second, we cut unpromising "leads" (prefixes in the combinatorial expansion) if they are not fully disjoint. This speeds up considerably the full search, by cutting down entire swathes of search space. If no prefix is cut, the full search is time-exponential.
What is the speed of this? It depends on the distribution of the v vectors and (strongly) on the distance r. In fact, it depends on the number of clusters found. To give an idea, with v 20 vectors picked in uniform in 2D, I observe roughly 30ms. For 40 vectors, typically around 100ms, but it can jump to over 2 seconds for certain values of r.
We implement a variation of combinations that has a check function to cut down unpromising prefixes:
from itertools import combinations
def _okall(tup, *args, **kwargs):
return True
def combinations_all(iterable, n0=1, n1=None, check=None, *args, **kwargs):
pool = tuple(iterable)
n = len(pool)
if n0 > n:
return
n1 = n if n1 is None else n1
check = _okall if check is None else check
if n0 < 1:
yield ()
n0 = 1
seed = list(combinations(range(n), n0-1))
for r in range(n0, n1+1):
prev_seed = seed
seed = []
for prefix in prev_seed:
for j in range(max(prefix, default=-1)+1, n):
indices = prefix + (j,)
tup = tuple(pool[i] for i in indices)
if check(tup, *args, **kwargs):
seed.append(indices)
yield tup
Example:
# list all combinations of [0,1,2,3] (of all sizes), excluding those starting
# by (0,1)
>>> list(combinations_all(range(4), check=lambda tup: tup[:2] != (0,1)))
[(0,),
(1,),
(2,),
(3,),
(0, 2),
(0, 3),
(1, 2),
(1, 3),
(2, 3),
(0, 2, 3),
(1, 2, 3)]
Now, we use that to find the minimal disjoint cover and return the indices of v1:
def check(tup, graph):
# refuses tup if any one is within reach of another
# optimization: tup[:-1] has already passed this test, so just check the last one
if len(tup) < 2:
return True
reach_of_last = graph[tup[-1]]
prefix = set(tup[:-1])
return prefix.isdisjoint(reach_of_last)
def brute_vq(v, r):
n = v.shape[0]
graph = get_graph(v, r)
to_part = set(graph)
singletons = set(range(n)).difference(to_part)
if not to_part:
return sorted(singletons)
for tup in combinations_all(to_part, check=check, graph=graph):
# here tup are indices that are far apart enough (guaranteed disjoint in reach)
# check if they fully cover to_part
cover = {j for i in tup for j in graph[i]}
if cover == to_part:
v1_idx = sorted(singletons.union(tup))
return v1_idx
Examples:
v = np.random.uniform(size=(20, 2))
r_s = [[0.2, 0.3], [0.4, 0.5]]
fig, axes = plt.subplots(nrows=len(r_s), ncols=len(r_s[0]),
figsize=(12,12), sharex=True, sharey=True)
for r, ax in zip(np.ravel(r_s), np.ravel(axes)):
idx = brute_vq(v, r)
ax.set_aspect('equal')
ax.scatter(v[:, 0], v[:, 1])
for i, p in enumerate(v):
ax.text(*p, str(i))
for i in idx:
circ = plt.Circle(v[i], radius=r, color='g', fill=False)
ax.add_patch(circ)
plt.show();
