使用 cupy 从 GPU 上的另一个矩阵创建距离矩阵答案

【问题标题】：Using cupy to create a distance matrix from another matrix on GPU使用 cupy 从 GPU 上的另一个矩阵创建距离矩阵
【发布时间】：2020-10-22 06:42:28
【问题描述】：

我使用 numpy 编写了代码，它采用大小为 (m x n) 的数组...行 (m) 是由 (n) 个特征组成的单个观察值...并创建一个大小为 (m x m) 的平方距离矩阵。该距离矩阵是给定观测值与所有其他观测值的距离。例如。第 0 行第 9 列是观测值 0 和观测值 9 之间的距离。

import numpy as np
#import cupy as np

def l1_distance(arr):
    return np.linalg.norm(arr, 1)

X = np.random.randint(low=0, high=255, size=(700,4096))

distance = np.empty((700,700))

for i in range(700):
    for j in range(700):
        distance[i,j] = l1_distance(X[i,:] - X[j,:])

我通过 umcommenting 第二个 import 语句在 GPU 上使用 cupy 尝试了此操作，但显然双 for 循环效率极低。 numpy 大约需要 6 秒，但 cupy 需要 26 秒。我明白为什么，但我并不清楚如何并行化这个过程。

我知道我需要编写某种归约内核，但我想不出如何通过对另一个数组元素的迭代操作来构造一个 Cupy 数组。

【问题讨论】：

你为什么不使用 scipy.spatial.distance.cdist？这必须比 CPU 循环快 5 到 10 倍
我明确地编码了我试图为读者理解做的事情。在现实生活中，我不会为 CPU 实现做双循环。然而，据了解，GPU 实现中没有像 scipy 这样的 elementwiseKernel 或 reductionKernel 的 Cuda 内核。这只是一个玩具数据集。我的真实数据集比这大得多，因此我为什么要尝试并行化它并将其移植到 GPU

标签： python gpu cupy

【解决方案1】：

在 A100 GPU 中使用广播 CuPy 需要 0.10 秒，而 NumPy 需要 6.6 秒

    for i in range(700):
        distance[i,:] = np.abs(np.broadcast_to(X[i,:], X.shape) - X).sum(axis=1)

这将向量化并使一个向量与所有其他向量的距离平行。

【讨论】：

我要睡觉了，我早上试试这个，但你能解释一下 broadcast_to() 在做什么吗？我想了解数据如何在 CPU 和 GPU 之间移动
broadcast_to 创建一个形状为 X 但其所有条目都是 X[i,:] 向量的数组，广播数组是一个视图，它使用步幅始终从原始数组，因此不需要额外的内存。数据根本没有从 CPU 移动到 GPU。它一直存在于 GPU 中。
我正在头脑中编写代码，这个实现看起来很棒。这么简单的事情，却让我难以捉摸哈哈。你有什么好的 GPU 编程参考资料吗？
在使用cupy时，你应该尽量向量化所有你可以做的操作，避免单项赋值，广播和摆弄步幅是非常重要的，但这与NumPy非常相似。你可以阅读ipython-books.github.io/46-using-stride-tricks-with-numpy 如果你想做你的自定义内核，CuPy 文档有一些例子docs.cupy.dev/en/stable/tutorial/kernel.html
我想知道这与 scikit-learn pairwise_distances、scipy cdist 等相比如何。