使用pairwise_distances_chunked计算最近邻搜索

Question

我有一个细长的数据矩阵(大小: 250,000 x 10)，我将其命名为X。我还有一个矢量p来测量我的数据点的质量。我的目标是为数据矩阵X中的每一行x计算以下函数

r(x) = min{ ||x-y| py>px，y in X}

在较小的数据集上，我将使用sklearn.metrics.pairwise_distances预先计算距离，如下所示：

from sklearn import metrics
n = len(X);

D_full = metrics.pairwise_distances(X);
r = np.zeros((n,1));
for i in range(n):
    r[i] = (D_full[i,p>p[i]]).min();

但是，上面的方法占用大量内存，因为我需要存储D_full：一个完整的n x n矩阵。看起来sklearn.metrics.pairwise_distances_chunked是解决这类问题的好工具，因为距离矩阵一次只存储一个块。我希望在如何使用它方面得到一些帮助，因为我目前还不熟悉生成器对象。假设我调用以下代码：

from sklearn import metrics
D = metrics.pairwise_distances_chunked(X);
D_chunk = next(D)

上面的代码产生D (一个生成器对象)和D_chunk (一个536xn数组)。根据我之前的方法，D_chunk是否对应于矩阵D_full的前536行？如果是，next(D_chunk)是否对应于接下来的536行？

非常感谢你的帮助。

Answer 1

这是一个可能的解决方案的概要，但缺少细节。简而言之，我会做以下几件事：

创建一个要查询的BallTree，并将大小为250000的min_quality_distance初始化为0。

对于k=2

1. 对于每个矢量，查找最近的k邻域(包括其自身)。
2. 如果找到k内距离最大的矢量，则具有足够的质量，请更新该点的min_quality_distance。
3. 对于剩余，请使用k=k+1

重复

在每次迭代中，我们必须查询较少的向量。这个想法是，在每一次迭代中，你都会在合适的条件下一点点地蚕食最近的邻居，每一步都会变得更容易。(容易50%？)我将展示如何进行第一次迭代，这样就可以构建循环了。

你可以做到；

import numpy as np
size = 250000

X = np.random.random( size=(size,10))
p = np.random.random( size=size)

并使用以下命令创建BallTree

from sklearn.neighbors import BallTree

tree = BallTree(X, leaf_size=10, metric='minkowski')

并查询它进行第一次迭代(这将需要大约5分钟)。

k_nearest = 2

distances, indici = tree.query(X, k=k_nearest, return_distance=True, dualtree=False, sort_results=True)

最近的k中最远的点的indici是

most_far_away_indici = indici[:,-1:]

它的质量

p[most_far_away_indici]

所以我们可以

quality_closeby = p[most_far_away_indici]

并使用以下命令检查是否有效

indici_sufficient_quality = quality_closeby > np.expand_dims(p, axis=1)

我们有

found_closeby = np.all( indici_sufficient_quality, axis=1 )

这是真的，我们在附近找到了足够的质量。

我们可以使用以下命令更新向量

distances_nearby = distances[:,-1:]

rx = np.zeros(size)
rx[found_closeby] = distances_nearby[found_closeby][:,0]

我们现在需要注意剩下的不幸之处，这些是

~found_closeby

所以

indici_not_found = ~found_closeby

和

distances, indici = tree.query(X[indici_not_found], k=3, return_distance=True, dualtree=False, sort_results=True)

等等。

我确信前几次循环将花费几分钟，但经过几次迭代后，速度将很快上升到几秒钟。

这是一个使用np.argwhere()等的小练习，以确保正确的indicis得到更新。

它可能不是最快的，但它是一种可行的方法。

Answer 2

由于无法知道某些块的大小，因此我建议使用np.ones_like而不是np.zeros。