Numpy，通过迭代替换广播答案

【问题标题】：Numpy, replace a broadcast by iterationNumpy，通过迭代替换广播
【发布时间】：2015-11-17 16:30:28
【问题描述】：

我有以下代码sn-p

def norm(x1, x2):
    return np.sqrt(((x1 - x2)**2).sum(axis=0))

def call_norm(x1, x2):
    x1 = x1[..., :, np.newaxis]
    x2 = x2[..., np.newaxis, :]
    return norm(x1, x2)

据我了解，每个x 代表 N 维空间中的一个点数组，其中 N 是数组的最终维度的大小（因此对于 3 空间中的点，最终维度的大小为 3）。它插入额外的维度并使用广播来生成这些点集的笛卡尔积，从而计算所有点对之间的距离。

x = np.array([[1, 2, 3],[1, 2, 3]])
call_norm(x, x)
array([[ 0.        ,  1.41421356,  2.82842712],
       [ 1.41421356,  0.        ,  1.41421356],
       [ 2.82842712,  1.41421356,  0.        ]])

（所以[1,1]和[2,2]之间的距离是1.41421356，正如预期的那样）

我发现对于中等大小的问题，这种方法可能会使用大量内存。我可以轻松地“去矢量化”问题并通过迭代替换它，但我希望这会很慢。我有一个（合理的）简单的折衷解决方案，我可以拥有矢量化的大部分速度优势，但没有内存损失？一些花哨的生成器技巧？

【问题讨论】：

你可以使用pdist & squareform: squareform(pdist(x.T))。
您是否总是计算同一数组中的对的距离，即call_norm 的两个输入是否相同？
@Divakar ：感谢您的回复！不，我有时需要有不同的论点
我想然后看看cdist。

标签： python numpy generator broadcast

【解决方案1】：

如果没有 numpy 向量化的内存损失，就无法进行这种计算。对于有效计算成对距离矩阵的特定情况，包倾向于通过在 C（例如scipy.spatial.distance）或 Cython（例如sklearn.metrics.pairwise）中实现事物来解决这个问题。

如果您想“手动”执行此操作，可以说，使用 numpy 样式的语法但不会导致内存损失，当前最佳选择可能是 dask.array，它可以自动构建和执行灵活任务使用 numpy 样式语法的批处理执行图。

以下是使用 dask 进行此计算的示例：

import dask.array as da

# Create the chunked data. This can be created
# from numpy arrays as well, e.g. x_dask = da.array(x_numpy)
x = da.random.random((100, 3), chunks=5)
y = da.random.random((200, 3), chunks=5)

# Compute the task graph (syntax just like numpy!)
diffs = x[:, None, :] - y[None, :, :]
dist = da.sqrt((diffs ** 2).sum(-1))

# Execute the task graph
result = dist.compute()
print(result.shape)
# (100, 200)

您会发现 dask 的内存效率比 NumPy 高得多，通常在计算上比 NumPy 更高效，并且还可以相对简单地并行/核外计算。

【讨论】：