更有效的方法来计算距离在纽比？

Question

我有一个问题，关于如何尽可能快地计算距离

def getR1(VVm,VVs,HHm,HHs):
    t0=time.time()
    R=VVs.flatten()[numpy.newaxis,:]-VVm.flatten()[:,numpy.newaxis]
    R*=R
    R1=HHs.flatten()[numpy.newaxis,:]-HHm.flatten()[:,numpy.newaxis]
    R1*=R1
    R+=R1
    del R1
    print "R1\t",time.time()-t0, R.shape, #11.7576191425 (108225, 10500) 
    print numpy.max(R) #4176.26290975
    # uses 17.5Gb ram
    return R


def getR2(VVm,VVs,HHm,HHs):
    t0=time.time()
    precomputed_flat = numpy.column_stack((VVs.flatten(), HHs.flatten()))
    measured_flat = numpy.column_stack((VVm.flatten(), HHm.flatten()))
    deltas = precomputed_flat[None,:,:] - measured_flat[:, None, :]
    #print time.time()-t0, deltas.shape # 5.861109972 (108225, 10500, 2)
    R = numpy.einsum('ijk,ijk->ij', deltas, deltas)
    print "R2\t",time.time()-t0,R.shape, #14.5291359425 (108225, 10500)
    print numpy.max(R) #4176.26290975
    # uses 26Gb ram
    return R


def getR3(VVm,VVs,HHm,HHs):
    from numpy.core.umath_tests import inner1d
    t0=time.time()
    precomputed_flat = numpy.column_stack((VVs.flatten(), HHs.flatten()))
    measured_flat = numpy.column_stack((VVm.flatten(), HHm.flatten()))
    deltas = precomputed_flat[None,:,:] - measured_flat[:, None, :]
    #print time.time()-t0, deltas.shape # 5.861109972 (108225, 10500, 2)
    R = inner1d(deltas, deltas)
    print "R3\t",time.time()-t0, R.shape, #12.6972110271 (108225, 10500)
    print numpy.max(R) #4176.26290975
    #Uses 26Gb
    return R


def getR4(VVm,VVs,HHm,HHs):
    from scipy.spatial.distance import cdist
    t0=time.time()
    precomputed_flat = numpy.column_stack((VVs.flatten(), HHs.flatten()))
    measured_flat = numpy.column_stack((VVm.flatten(), HHm.flatten()))
    R=spdist.cdist(precomputed_flat,measured_flat, 'sqeuclidean') #.T
    print "R4\t",time.time()-t0, R.shape, #17.7022118568 (108225, 10500)
    print numpy.max(R) #4176.26290975
    # uses 9 Gb ram
    return R

def getR5(VVm,VVs,HHm,HHs):
    from scipy.spatial.distance import cdist
    t0=time.time()
    precomputed_flat = numpy.column_stack((VVs.flatten(), HHs.flatten()))
    measured_flat = numpy.column_stack((VVm.flatten(), HHm.flatten()))
    R=spdist.cdist(precomputed_flat,measured_flat, 'euclidean') #.T
    print "R5\t",time.time()-t0, R.shape, #15.6070930958 (108225, 10500)
    print numpy.max(R) #64.6240118667
    # uses only 9 Gb ram
    return R

def getR6(VVm,VVs,HHm,HHs):
    from scipy.weave import blitz
    t0=time.time()
    R=VVs.flatten()[numpy.newaxis,:]-VVm.flatten()[:,numpy.newaxis]
    blitz("R=R*R") # R*=R
    R1=HHs.flatten()[numpy.newaxis,:]-HHm.flatten()[:,numpy.newaxis]
    blitz("R1=R1*R1") # R1*=R1
    blitz("R=R+R1") # R+=R1
    del R1
    print "R6\t",time.time()-t0, R.shape, #11.7576191425 (108225, 10500) 
    print numpy.max(R) #4176.26290975
    return R

结果如下：

^{pr2}$

虽然最后一个给出了sqrt（（VVm VVs）^2+（HHm HHs）^2，而其他的给出了（VVm VVs）^2+（HHm HHs）^2，但这并不重要，因为在我的代码中，我对每个i取R[i，：]的最小值，并且sqrt无论如何都不会影响最小值，（如果我对距离感兴趣，我只取sqrt（value），而不是在整个数组上执行sqrt，因此实际上没有时间差。在

问题仍然存在：为什么第一个解决方案是最好的（第二个和第三个方案比较慢的原因是因为Delta=。。。需要5.8秒（这也是为什么这两种方法需要26Gb），为什么sqeuclide比euclidean慢？在

sqeuclidean应该做（VVm VVs）^2+（HHm HHs）^2，而我认为它做了一些不同的事情。有人知道如何找到该方法的源代码（C或其他位于底部的代码）？我认为它是sqrt（（VVm VVs）^2+（HHm HHs）^2）^2（我能想到为什么它会比（VVm VVs）^2+（HHm HHs）^2慢的唯一原因-我知道这是一个愚蠢的原因，有人有更合理的理由吗？）在

既然我对C一无所知，我该如何将它与西皮·韦伯？这些代码是否像python一样可以正常编译？或者我需要安装专门的东西吗？在

编辑：好的，我试过了席皮.韦伯.闪电战，（R6方法），这个速度稍微快一点，但是我假设比我懂更多C的人仍然可以提高这个速度？我只是用a+=b或*=，查找它们在C中的位置，然后把它们放在blitz语句中，但是我想如果我把带有flatten和newaxis语句的行也放在C中，那么它也应该更快，但我不知道我该怎么做（懂C的人也许能解释一下？）。现在，blitz和我的第一个方法之间的区别还不够大，不足以真正由C vs numpy引起？在

我想其他的方法比如deltas=。。。也可以快得多，当我把它放在C语言里的时候？在

Answer 1

每当你有乘法和和运算时，试着使用点积函数或np.einsum。由于要预先分配阵列，而不是为水平坐标和垂直坐标使用不同的阵列，因此请将两者堆叠在一起：

precomputed_flat = np.column_stack((svf.flatten(), shf.flatten()))
measured_flat = np.column_stack((VVmeasured.flatten(), HHmeasured.flatten()))
deltas = precomputed_flat - measured_flat[:, None, :]

从这里，最简单的方法是：

^{pr2}$

你也可以尝试一下：

from numpy.core.umath_tests import inner1d
dist = inner1d(deltas, deltas)

---

当然还有SciPy的空间模块^{}：

from scipy.spatial.distance import cdist
dist = cdist(precomputed_flat, measured_flat, 'euclidean')

---

编辑 我无法在如此大的数据集上运行测试，但这些时间安排非常有启发性：

len_a, len_b = 10000, 1000

a = np.random.rand(2, len_a)
b =  np.random.rand(2, len_b)
c = np.random.rand(len_a, 2)
d = np.random.rand(len_b, 2)

In [3]: %timeit a[:, None, :] - b[..., None]
10 loops, best of 3: 76.7 ms per loop

In [4]: %timeit c[:, None, :] - d
1 loops, best of 3: 221 ms per loop

对于上述较小的数据集，我可以通过在内存中以不同的方式排列数据，scipy.spatial.distance.cdist使其与{}匹配：

precomputed_flat = np.vstack((svf.flatten(), shf.flatten()))
measured_flat = np.vstack((VVmeasured.flatten(), HHmeasured.flatten()))
deltas = precomputed_flat[:, None, :] - measured_flat

import scipy.spatial.distance as spdist
from numpy.core.umath_tests import inner1d

In [13]: %timeit r0 = a[0, None, :] - b[0, :, None]; r1 = a[1, None, :] - b[1, :, None]; r0 *= r0; r1 *= r1; r0 += r1
10 loops, best of 3: 146 ms per loop

In [14]: %timeit deltas = (a[:, None, :] - b[..., None]).T; inner1d(deltas, deltas)
10 loops, best of 3: 145 ms per loop

In [15]: %timeit spdist.cdist(a.T, b.T)
10 loops, best of 3: 124 ms per loop

In [16]: %timeit deltas = a[:, None, :] - b[..., None]; np.einsum('ijk,ijk->jk', deltas, deltas)
10 loops, best of 3: 163 ms per loop