Python,成对"距离",需要一种快速的方法来实现它



在我博士学位的一个附带项目中,我参与了用Python建模一些系统的任务。就效率而言,我的程序在下面的问题中遇到了瓶颈,我将在一个最小工作示例中揭示这个问题。

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

我需要计算一个成对的最小分段间距离。两个线段之间的最小距离的解析表达式在该源中找到。至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个分量。我认为我无法将距离计算转化为二进制函数。

欢迎提供任何意见。

您可以使用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()

返回CCD_ 2。我也很确定算法中隐藏着矩阵乘法,所以应该有一些进一步加速(以及清理)的潜力。

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

这更像是一个元答案,至少对初学者来说是这样。你的问题可能已经在"我的程序遇到瓶颈"one_answers"我意识到这是非常低效的"。

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

那么,如果您对以上所有问题都有答案,并且需要进一步优化您的代码——您当前代码的瓶颈在哪里??你对它进行了简介吗?环路的主体可能比环路本身的评估更重吗?如果是这样,那么"循环"就不是你的瓶颈,你一开始就不需要担心嵌套的循环。一开始可以优化身体,可能是通过对数据提出非正统的矩阵表示,这样你就可以一步完成所有这些计算——例如,通过矩阵乘法。如果你的问题不能通过有效的线性代数运算来解决,你可以开始写一个C扩展,或者使用Cython或PyPy(它最近才得到一些基本的numpy支持!)。优化有无限的可能性——真正的问题是:你离实际解决方案有多近,你需要优化多少,以及你愿意投入多少精力。

免责声明:在我的博士学位上,我也用scipy/numpy做过非规范的成对距离的事情;-)。对于一个特定的距离度量,我最终用简单的Python编码了"成对"部分(即,我也使用了双嵌套循环),但我花了一些精力使主体尽可能高效(使用I)我的问题的神秘矩阵乘法表示和ii)使用bottleneck的组合)。

您可以这样使用它:

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))

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

最新更新