более эффективный способ вычисления расстояния в numpy?

Question

У меня есть вопрос о том, как вычислить расстояния в NumPy так быстро, как это возможно,

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 

результаты в следующих случаях:

R1 11.7737319469 (108225, 10500) 4909.66881791 
R2 15.1279799938 (108225, 10500) 4909.66881791 
R3 12.7408981323 (108225, 10500) 4909.66881791 
R4 17.3336868286 (10500, 108225) 4909.66881791 
R5 15.7530870438 (10500, 108225) 70.0690289494 
R6 11.670968771 (108225, 10500) 4909.66881791 

В то время как последняя дает SQRT ((VVM-VVS)^2 + (HHm-HHs)^2), а остальные дают (VVm-VVs)^2 + (HHm-HHs)^2. Это не очень важно, так как иначе в моем коде я беру минимум из R [i ,:] для каждого i, и sqrt никак не влияет на минимальное значение (и если меня интересует расстояние, я просто беру sqrt (value) вместо того, чтобы делать sqrt по всему массиву, так что на самом деле нет времени разница из-за этого.

Вопрос остается: как первое решение является лучшим (причина, по которой вторая и третья медленнее, потому что deltas = ... занимает 5.8 сек. (Что также объясняет, почему эти два метода принимают 26 Гб)), И почему sqeuclidean медленнее, чем евклид?

sqeuclidean должен просто делать (VVm-VVs)^2 + (HHm-HHs)^2, в то время как я думаю, что он делает что-то другое. Кто-нибудь знает, как найти исходный код (C или что-то еще внизу) этого метода? Я думаю, что это делает sqrt ((VVm-VVs)^2 + (HHm-HHs)^2)^2 (единственная причина, по которой я могу думать, почему она будет медленнее, чем (VVm-VVs)^2 + (HHm-HHs)^2 - Я знаю, что это глупая причина, кто-нибудь получил более логичный?)

Так как я ничего не знаю о C, как бы это сделать с помощью scipy.weave? и является ли этот код компилируемым, как правило, с помощью python? или мне нужны специальные материалы для этого?

Редактирование: хорошо, я попробовал его с scipy.weave.blitz, (метод R6), и это немного быстрее, но я предполагаю, что кто-то, кто знает больше C, чем я, может еще улучшить эту скорость? Я просто взял строки, которые имеют форму a = = b или * =, и посмотрел, как они будут на C, и помещает их в блиц-предложение, но я думаю, если бы я поместил строки с операторами с flatten и newaxis в C также, что он должен идти быстрее, но я не знаю, как я могу это сделать (кто-то, кто знает C, может объяснить?). Прямо сейчас, разница между материалом с блицом и моим первым методом невелика, чтобы действительно быть вызвана 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, :] 

здесь, проще было бы:

dist = np.einsum('ijk,ijk->ij', deltas, deltas) 

Вы также можете попробовать что-то как:

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

---

EDIT Я не могу запустить тесты на таком большом наборе данных, но эти тайминги весьма поучительно:

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 и сопоставьте его с inner1d, упорядочив данные по-разному в памяти:

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