Python，成对“距离”，需要一种快速的方法来做到这一点答案

【问题标题】：Python, Pairwise 'distance', need a fast way to do itPython，成对“距离”，需要一种快速的方法来做到这一点
【发布时间】：2015-02-24 10:36:42
【问题描述】：

对于我博士期间的一个副项目，我从事了用 Python 对某些系统进行建模的任务。在效率方面，我的程序在以下问题中遇到了瓶颈，我将在一个最小的工作示例中展示。

我处理大量由它们的 3D 起点和终点编码的片段，因此每个片段由 6 个标量表示。

我需要计算成对的最小段间距离。两个线段之间最小距离的解析表达式见此source。致 MWE：

import numpy as np
N_segments = 1000
List_of_segments = np.random.rand(N_segments, 6)

Pairwise_minimal_distance_matrix = np.zeros( (N_segments,N_segments) )
for i in range(N_segments):
    for j in range(i+1,N_segments): 

        p0 = List_of_segments[i,0:3] #beginning point of segment i
        p1 = List_of_segments[i,3:6] #end point of segment i
        q0 = List_of_segments[j,0:3] #beginning point of segment j
        q1 = List_of_segments[j,3:6] #end point of segment j
        #for readability, some definitions
        a = np.dot( p1-p0, p1-p0)
        b = np.dot( p1-p0, q1-q0)
        c = np.dot( q1-q0, q1-q0)
        d = np.dot( p1-p0, p0-q0)
        e = np.dot( q1-q0, p0-q0)
        s = (b*e-c*d)/(a*c-b*b)
        t = (a*e-b*d)/(a*c-b*b)
        #the minimal distance between segment i and j
        Pairwise_minimal_distance_matrix[i,j] = sqrt(sum( (p0+(p1-p0)*s-(q0+(q1-q0)*t))**2)) #minimal distance

现在，我意识到这非常低效，这就是我在这里的原因。我已经广泛研究了如何避免循环，但我遇到了一些问题。显然，这种计算最好用 python 的cdist 来完成。但是，它可以处理的自定义距离函数必须是二进制函数。在我的情况下这是一个问题，因为我的向量的长度特别是 6，并且必须按位拆分为它们的第一个和最后 3 个分量。我不认为我可以将距离计算转换为二进制函数。

感谢任何输入。

【问题讨论】：

可以在这里提供octree（common uses 提到维基百科上的最近邻搜索）的帮助吗？
你读过 Lumelsky 的On fast computation of distance between line segments [PDF 警告] 吗？你的实现与它相比如何？（可以找到更通用的方法here）
谢谢你的链接，我去看看。
@Michael Foukarakis，我很快阅读了这篇论文，这是我自己分析得出的。它没有具体说明如何加速计算。它只是概述了一种处理特殊情况的聪明方法。虽然很好读

标签： python performance binary distance

【解决方案1】：

您可以使用 numpy 的矢量化功能来加快计算速度。我的版本一次计算距离矩阵的所有元素，然后将对角线和下三角形设置为零。

def pairwise_distance2(s):
    # we need this because we're gonna divide by zero
    old_settings = np.seterr(all="ignore")

    N = N_segments # just shorter, could also use len(s)

    # we repeat p0 and p1 along all columns
    p0 = np.repeat(s[:,0:3].reshape((N, 1, 3)), N, axis=1)
    p1 = np.repeat(s[:,3:6].reshape((N, 1, 3)), N, axis=1)
    # and q0, q1 along all rows
    q0 = np.repeat(s[:,0:3].reshape((1, N, 3)), N, axis=0)
    q1 = np.repeat(s[:,3:6].reshape((1, N, 3)), N, axis=0)

    # element-wise dot product over the last dimension,
    # while keeping the number of dimensions at 3
    # (so we can use them together with the p* and q*)
    a = np.sum((p1 - p0) * (p1 - p0), axis=-1).reshape((N, N, 1))
    b = np.sum((p1 - p0) * (q1 - q0), axis=-1).reshape((N, N, 1))
    c = np.sum((q1 - q0) * (q1 - q0), axis=-1).reshape((N, N, 1))
    d = np.sum((p1 - p0) * (p0 - q0), axis=-1).reshape((N, N, 1))
    e = np.sum((q1 - q0) * (p0 - q0), axis=-1).reshape((N, N, 1))

    # same as above
    s = (b*e-c*d)/(a*c-b*b)
    t = (a*e-b*d)/(a*c-b*b)

    # almost same as above
    pairwise = np.sqrt(np.sum( (p0 + (p1 - p0) * s - ( q0 + (q1 - q0) * t))**2, axis=-1))

    # turn the error reporting back on
    np.seterr(**old_settings)

    # set everything at or below the diagonal to 0
    pairwise[np.tril_indices(N)] = 0.0

    return pairwise

现在让我们试一试。以你的例子，N = 1000，我得到了一个时间

%timeit pairwise_distance(List_of_segments)
1 loops, best of 3: 10.5 s per loop

%timeit pairwise_distance2(List_of_segments)
1 loops, best of 3: 398 ms per loop

当然，结果是一样的：

(pairwise_distance2(List_of_segments) == pairwise_distance(List_of_segments)).all()

返回True。我也很确定在算法的某个地方隐藏了一个矩阵乘法，所以应该有进一步加速（以及清理）的潜力。

顺便说一句：我尝试过简单地先使用 numba，但没有成功。不过不知道为什么。

【讨论】：

非常感谢。当我读到你的第一句话时，我想“当然”，但是看到实施是如何完成的，我想你节省了我几个小时的工作：为此我感谢你！据我了解，Numba 是基于 GPU 计算的。根据您的经验，在矢量大小相似（N_segments ~= 1E3-E4）的情况下，将数据推送到显卡的开销是否超过了加速，或者在实施之前无法判断？
我很快就遇到了内存问题（N_segments ~= 10k）。我会调查什么是内存密集型并回复你（我可以让它在 matlab 上运行高达 30k）
@Mathusalem 一个带有N=10000 的双精度矩阵N*N 需要大约0.75GB 的内存。这个函数创建了一堆。我的方法是（因为我对这个问题的了解还不够多，无法设计出更好的算法）将问题拆分成可管理大小的块并一个接一个地计算它们（或者甚至并行计算，如果你愿意的话）。
我也尝试从程序中提取循环部分并进行 JIT 编译，但我没有获得任何速度。你有同样的问题吗？我怀疑
对于您的 numba 评论。我一直在研究它。我认为缺乏加速的原因在于 numba 中仍然没有太多对 numpy 的支持。特别是，numba 难以处理许多 numpy 函数，如 numpy.dot、numpy.sqrt，因为它不知道这些函数的返回类型，因此无法优化。链接：github.com/numba/numba/issues/251

【解决方案2】：

这更像是一个元答案，至少对于初学者来说。您的问题可能已经在“我的程序遇到瓶颈”和“我意识到这非常低效”。

效率极低？用什么衡量标准？你有比较吗？您的代码是否太慢而无法在合理的时间内完成？什么对你来说是合理的？你能在这个问题上投入更多的计算能力吗？同样重要——您是否使用适当的基础架构来运行您的代码（使用供应商编译器编译的 numpy/scipy，可能支持 OpenMP）？

那么，如果您对上述所有问题都有答案，并且需要进一步优化您的代码——您当前代码的瓶颈究竟在哪里？你有介绍吗？它的循环体可能比循环本身的评估重得多？如果是这样，那么“循环”不是您的瓶颈，您首先不必担心嵌套循环。首先优化主体，可能通过提出数据的非正统矩阵表示，以便您可以一步执行所有这些单一计算 - 例如通过矩阵乘法。如果您的问题无法通过有效的线性代数运算解决，您可以开始编写 C 扩展或使用 Cython 或使用 PyPy（最近才获得一些基本的 numpy 支持！）。优化的可能性无穷无尽——问题实际上是：您离实际解决方案有多近，您需要优化多少，以及您愿意投入多少努力。

免责声明：我也为我的博士学位使用 scipy/numpy 完成了非规范的成对距离的工作；-)。对于一个特定的距离度量，我最终用简单的 Python 编写了“成对”部分（即，我也使用了双重嵌套循环），但花了一些努力使主体尽可能高效（结合 i）a我的问题的神秘矩阵乘法表示，ii）使用bottleneck）。

【讨论】：

【解决方案3】：

你可以像这样使用它：

def distance3d (p, q):
    if (p == q).all ():
        return 0

    p0 = p[0:3]
    p1 = p[3:6]
    q0 = q[0:3]
    q1 = q[3:6]

    ...  # Distance computation using the formula above.

print (distance.cdist (List_of_segments, List_of_segments, distance3d))

不过，它似乎并没有更快，因为它在内部执行相同的循环。

【讨论】：