大型NumPy数组的成对距离(分块?)

6
问题: 我有一个大约为 [350000,1] 的向量,我希望计算成对距离。这会导致一个整数数据类型的[350000, 350000]矩阵,无法适应RAM。最终我想要得到一个布尔值(适合RAM),所以我目前正在一次处理一个元素,但效率不高。

编辑:标准的sklearn和scipy函数由于数据的大小而无法使用 - 但是如果我可以将其分块以使用硬盘,则应该能够使用这些函数。

可视化问题: [a_1,a_2,a_3] ^ t -> [[a_1-a_1,a_1-a_2,a_1-a_3],[a_2-a_1,a_2-a_2,a_2-a_3],[a_3-a_1,a_3-a_2,a_3-a_3]]
请注意,只有上三角需要计算,因为在取绝对值时对称。 需要分块或其他解决方案的矢量化代码: 我已经找到了一种计算所有点之间距离(差)的方法,并在小矩阵上使用广播进行了计算,但需要一种方法能够在更大的矩阵上进行计算而不会超过 RAM 限制。

或者可能有更快的改进 MWE(最小工作示例)的方法?

distMatrix = np.absolute((points[np.newaxis, :, :] - points[:, np.newaxis, :])[:, :, 0])

其他尝试: 我尝试使用dask和memmap,但仍然遇到内存错误,因此一定是做了些错误的操作。 我还尝试了memmap并手动分块数据,但无法获得完整的结果,因此需要任何帮助。

当前方法的最小可复现示例:


## Data ##
#Note that the datatype and code may not match up exactly as just creating to demonstrate. Essentially want to take first column and create distance matrix with itself through subtracting, and then take 2nd and 3rd column and create euclidean distance matrix.

data = np.random.randint(1, 5, size=[350001,3])
minTime = 3
maxTime = 4
minDist = 1
maxDist = 2

### CODE ###
n = len(data)

for i in trange(n):
    for j in range(i+1, n):
        #Within time threshold?
        if minTime <= (data[j][idxT] - data[i][idxT]) <= maxTime:
            #Within distance threshold?
            xD = math.pow(data[j][idxX] - data[i][idxX], 2)
            yD = math.pow(data[j][idxY] - data[i][idxY], 2)
            d = math.sqrt(xD + yD)
            #If within  threshold then
            if minDist <= d <= maxDist:
                #DO SOMETHING

原因: 我有大约350000个点的时间,x轴和y轴坐标向量。我想计算所有时间点之间的距离(简单减法)以及每个(x,y)点之间的欧几里得距离。然后,我想能够识别所有在时间和距离阈值内的点对,并生成一个布尔值。


这个回答解决了你的问题吗?Python中最快的成对距离度量 - taha
1
你对欧几里得距离还是基于绝对值的距离感兴趣? - Divakar
@Divakar 对两者都感兴趣,但如果我能让基于绝对值的算法正常工作,那么我就可以从那里计算欧几里得距离。 - Daniel J
你能添加一个带有示例数据的最小可运行代码吗? - Divakar
@Divakar 已修改。 - Daniel J
idxTidxXidxY 的示例代表值是什么? - Divakar
1个回答

2
你可以将数组分成较小的部分,并单独计算每对之间的距离。
splits = np.array_split(data, 10)
for i in range(len(splits)):
    for j in range(i, len(splits)):
        m = scipy.spatial.distance.cdist(splits[i], splits[j])
        # do something with m

作为大部分计算发生在scipy中,Python循环的开销将最小化。
如果您的布尔数组适合内存,并且您尝试查找某个范围内的值,则可以这样做。
import numpy as np
import scipy.spatial.distance


boolean = np.zeros((350, 350), dtype=np.bool_)
a = np.random.randn(350, 2)
splits = np.array_split(a, 10)
shift = splits[0].shape[0]
minDist = -0.5
maxDist = +0.5
for i in range(len(splits)):
    for j in range(i, len(splits)):
        m = scipy.spatial.distance.cdist(splits[i], splits[j])
        masked = (minDist <= m) & (m <= maxDist)
        boolean[i * shift: (i + 1) * shift, j * shift : (j + 1) * shift] = masked
        boolean[j * shift : (j + 1) * shift, i * shift: (i + 1) * shift] = masked.T

@DanielJ 我更新了答案。这是你需要的吗? - V. Ayrat
@v-ayrat 当我使用一个长度为350,000的向量运行时,速度非常慢 - 我计算出按这个速度需要大约90小时,而且这仅仅是针对一个布尔值,不包括我的其他操作,这些操作总共目前需要20小时。我也尝试过调整块的大小。有没有什么方法可以加快速度?逐个迭代会更快一些 :( - Daniel J
我怀疑这不是计算时间的问题。根据您的系统,最多应该需要半个小时。您可以尝试注释掉带有 boolean = ...boolean[j * ... 的行,这意味着您只进行计算而不存储值在内存中。如果速度快得多,则问题出在内存上。即使是那么大的布尔矩阵也应该约为100 Gb,因此它不会存储在RAM中,在计算过程中会不断交换。在这种情况下,您有两个选择:关闭所有侧面程序,并希望此数组适合您的RAM - 我怀疑这样做有帮助。第二:您是否继续对每个进行进一步计算 - V. Ayrat
分块处理,无需存储所有数据。 - V. Ayrat
让我们在聊天中继续这个讨论 - Daniel J
显示剩余3条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接