V
(例如一个64维的向量,如[1,2,...64]
)。从向量集合 V
中,我想找到一个由 n
个不同向量组成的子集 V1
,其中:
- 子集
V1
中的任意两个向量之间的最小距离为d
;并且 - 对于向量集合
V
中的任何元素w
,都存在一个元素w'
在子集V1
中,使得它与w
的距离小于d
。
什么算法能够有效地确定一个包含 n
个不同向量的子集,其中 n
的值是最小的?
距离可以是欧几里得距离或余弦相似度(当然不能混用)。
def greedy(V, d):
cover = set()
uncovered = set(V)
# invariant: all(distance(x,y) > d for x in cover for y in uncovered)
while uncovered:
x = uncovered.pop()
cover.add(x)
uncovered = {y for y in uncovered if distance(x,y) > d}
uncovered.pop()
调用来使这个贪心算法更好,例如选择一个能够“覆盖”最多剩余点的向量。v
中的每个点都可以被至少一个距离 <r
的 v1
中的点到达,而且我们还希望v1
中的所有点相互之间的距离都大于 r
。 - Pierre D首先,我们可以推导出一个无向图,其中顶点是向量v
的元素,边缘是不超过r
的点对。这不需要距离成为一个适当的度量:它不必满足三角不等式,只需对称(dist(a,b) == dist(b,a)
)并且为0,如果两个点相等。
在完成这一步之后,我们可以完全忽略距离,并专注于对图进行分区。正如其他人所指出的那样,该分区是一个Vertex Cover问题。但是,它有一个扭曲的地方:我们要求覆盖中的所有顶点都是不相交的(即:v1
中的向量不能彼此在距离r
以内)。
为了获得图形,我们可以使用高效的KDTree
仅计算一次最近邻结构。特别地,我们将使用kd.sparse_distance_matrix(kd, r)
来获取所有距离小于r
的点对:
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
(注意:如果我们想要使用其他距离,.sparse_distance_matrix()
有一个参数p
)。
图的分区可以通过多种方式完成。@DanR展示了一种贪婪的方法,非常快速,但通常在找到v1
的大小时不够优秀。
下面展示了一种暴力方法,如果存在最小解,则保证找到最小解。对于相对较少的向量来说是足够的,并且可以为其他启发式算法的最优解提供基准。
为了加速组合搜索,我们首先观察到通常一些点与任何其他点的距离不超过r
(即上面找到的图并未表示v
的所有点)。我们可以将这些点放在一边,因为它们必然是v1
的一部分,但不会干扰搜索的其余部分(下面称之为“单例”)。其次,如果它们不完全不相交,我们会削减无望的“线索”(组合扩展中的前缀)。这通过切断整个搜索空间的大片区域来显著加速完整搜索。如果没有削减前缀,则完整搜索的时间呈指数增长。
这个的速度是多少?它取决于v
向量的分布,以及与距离r
(强烈)相关。实际上,它取决于找到的聚类数。为了给出一个想法,在2D中均匀选择v
20个向量,我观察到大约30毫秒。对于40个向量,通常在100毫秒左右,但对于某些r
值,它可能会跳到超过2秒。
我们实现了combinations
的变体,它具有一个check
函数来削减不太有希望的前缀:
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
例子:
# 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)]
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
示例:
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();
idx
是 v
向量列表的索引,用于生成 v1
。在向量量化(VQ)和聚类术语中,这样的解决方案 v1
被称为码本或原型。每张图片显示了给定 r
值的解决方案,它绘制了每个原型的绿色圆圈(以其为中心,半径为 r
),以及 v
向量作为蓝点,以便通过视觉验证解决方案。 - Pierre D我的看法:
A B C D | min,max
A 0 3 2 1 | 1,3
B 3 0 5 4 | 3,5
C 2 5 0 6 | 2,6
D 1 4 6 0 | 1,6 # note: D-C violates the triangle inequality
距离矩阵,重新标记:
A D B C | min,max
A 0 1 3 2 | 1,3
D 1 0 4 6 | 1,6
B 3 4 0 5 | 3,5
C 2 6 5 0 | 2,6
现在,对于每一行(或列),取最小值的最大值(或最大值的最小值...) 距离为3时,西北部分完全连接,其余部分仍然部分连接 (对于d=4,B也将进入子集)
A D | B C | min,max
A 0 1 | 3 2 | 1,3
D 1 0 | 4 6 | 1,6
----+ ----------
B 3 4 | 0 5 | 3,5
C 2 6 | 5 0 | 2,6
v1
的向量是什么? - Pierre Dkd = KDTree(v)
,然后kd.sparse_distance_matrix(kd,r)
是一种非常有效的方法,用于查找彼此之间距离更近的所有点对。然后这只是一个带有扭曲的图形覆盖问题:覆盖中相互排斥。 - Pierre D
v1
成为一组新向量(不一定是v
的成员)。我所知道的一个例外是k-medoids。但这并不完全符合您的问题,只是探索的一个领域。其他感兴趣的领域包括分层聚类。 - Pierre D