bardziej skuteczny sposób obliczania odległości w numpy?

Question

mam pytanie, w jaki sposób obliczyć odległości w numpy tak szybko, jak to możliwe,

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 

rezultaty w następujących godzinach:

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 

Chociaż ten ostatni daje sqrt ((wm-VVS)^2 + (HHm-HH)^2), podczas gdy inne dają (VVm-VV)^2 + (HHm-HH)^2, To nie jest naprawdę ważne, ponieważ w przeciwnym razie dalej w moim kodzie biorę minimum R [i ,:] dla każdego i i sqrt nie wpływa na wartość minimalną w każdym razie (i jeśli interesuje mnie odległość, po prostu biorę sqrt (wartość), zamiast robić sqrt w całej tablicy, więc nie naprawdę nie ma czasu Różnica z tego powodu.

Pozostaje pytanie: jak to się stało, że pierwsze rozwiązanie jest najlepsze (powód, dla którego drugi i trzeci są wolniejsze, ponieważ deltas = ... zajmuje 5,8 sekund, (dlatego też te dwie metody pobierają 26 Gb), I dlaczego sqeuclidean wolniej niż euklidesa?

sqeuclidean powinien po prostu zrobić (VVm-VV)^2 + (HHm-HH)^2, podczas gdy myślę, że robi coś innego. Ktoś wie, jak znaleźć kod źródłowy (C lub cokolwiek jest na dole) tej metody? Myślę, że to sqrt ((VVm-VV)^2 + (HHm-HH)^2)^2 (jedyny powód, dla którego mogę myśleć, dlaczego byłby wolniejszy niż (VVm-VV)^2 + (HHm-HH)^2 - Wiem, że to głupi powód, czy ktoś ma bardziej logiczny?)

Ponieważ nic nie wiem o C, w jaki sposób wstawiłbym to za pomocą scipy.weave? i czy ten kod jest kompilowany normalnie tak jak w przypadku pythona? czy potrzebuję specjalnych rzeczy do tego zainstalowanych?

Edycja: ok, próbowałem go z scipy.weave.blitz, (metoda R6), a to trochę szybciej, ale zakładam, że ktoś, kto zna więcej C niż ja, może jeszcze poprawić tę prędkość? Po prostu wziąłem linie, które są w formie a + = b lub * =, i sprawdziłem, jak będą w C, i umieściłem je w oświadczeniu blitz, ale myślę, że jeśli umieściłem linie z instrukcjami spłaszczonymi i newaxis w C również, że powinien też działać szybciej, ale nie wiem, jak to zrobić (ktoś, kto zna C może wyjaśnić?). W tej chwili, różnica między materiałem a blitzem i moją pierwszą metodą nie jest wystarczająco duża, by naprawdę być spowodowana przez C vs numpy?

Sądzę, że inne metody, takie jak deltas = ..., mogą pójść o wiele szybciej, kiedy umieściłbym je w C?

Answer 1

Ilekroć masz mnożenia i sumy, spróbuj użyć jednej z funkcji produktu dot lub np.einsum.Ponieważ jesteś preallocating swoje tablice, zamiast różne tablice dla współrzędnych poziomych i pionowych, układać je obie razem:

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

stąd, najprostsze byłoby:

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

Można również spróbować czegoś jak:

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

---

Istnieje oczywiście również moduł przestrzenny scipy za cdist:

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

---

EDIT nie mogę uruchomić testy na tak dużą zbioru danych, ale te czasy są dość pouczające:

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 

dla wyżej mniejszym zbiorze, mogę dostać nieznaczne przyspieszenie w stosunku do metody z scipy.spatial.distance.cdist i dopasowanie go do inner1d, poprzez ustawienie danych inaczej w pamięci:

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