modo più efficiente per calcolare la distanza in numpy?

Question

Ho una domanda su come calcolare le distanze in NumPy più velocemente che può,

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 

risultati nei seguenti orari:

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 

mentre l'ultimo dà sqrt ((VVM-VVs)^2 + (HHm-HHs)^2), mentre gli altri danno (VVm-VVs)^2 + (HHm-HHs)^2, Questo non è veramente importante, poiché altrimenti più avanti nel mio codice prendo il minimo di R [i ,:] per ogni i, e sqrt non influenza comunque il valore minimo (e se sono interessato alla distanza, prendo solo sqrt (valore), invece di eseguire sqrt sull'intero array, quindi non è davvero il momento differenza a causa di ciò.

La domanda rimane: come mai la prima soluzione è la migliore, (la ragione per la seconda e la terza sono più lenti è perché delta = ... prende 5,8 secondi, (che è anche il motivo per cui questi due metodi prendono 26GB)), E perché lo sqeuclide è più lento dell'euclideo?

sqeuclidean dovrebbe solo fare (VVm-VVs)^2 + (HHm-HHs)^2, mentre penso che faccia qualcosa di diverso. Qualcuno sa come trovare il codice sorgente (C o qualunque cosa sia in fondo) di quel metodo? Penso che faccia sqrt ((VVm-VVs)^2 + (HHm-HHs)^2)^2 (l'unica ragione per cui posso pensare perché sarebbe più lento di (VVm-VVs)^2 + (HHm-HHs)^2 - So che è uno stupido motivo, qualcuno ne ha avuto uno più logico?)

Dato che non so nulla di C, come potrei metterlo in linea con scipy.weave? ed è quel codice compilabile normalmente come fai con Python? o ho bisogno di cose speciali installate per quello?

Modifica: ok, l'ho provato con scipy.weave.blitz, (metodo R6), e questo è leggermente più veloce, ma presumo che qualcuno che ne sa più C di me possa ancora migliorare questa velocità? Ho appena preso le linee che sono nella forma a + = b o * =, e ho cercato come sarebbero state in C, e le ho inserite nell'istruzione blitz, ma suppongo che se metto le righe con le frasi con appiattire e newax in C, anche questo dovrebbe andare più veloce, ma non so come posso farlo (qualcuno che sa che C può spiegare?). In questo momento, la differenza tra le cose con blitz e il mio primo metodo non sono abbastanza grandi per essere realmente causate da C vs Numpy, immagino?

Immagino che gli altri metodi come con deltas = ... possano andare anche molto più velocemente, quando lo inserirò in C?

Answer 1

Ogni volta che si dispone di moltiplicazioni e somme, provare a utilizzare una delle funzioni del prodotto punto o np.einsum.Dal momento che si preallocare gli array, piuttosto che avere le matrici diverse per le coordinate orizzontali e verticali, impilare tutti e due insieme:

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

Da qui, il più semplice sarebbe:

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

Si potrebbe anche provare qualcosa come:

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

---

V'è, naturalmente, anche il modulo spaziale di SciPy cdist:

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

---

EDIT non riesco a eseguire test su un grande insieme di dati, quali, ma questi tempi sono piuttosto illuminante:

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 

Per l'insieme di dati di cui sopra più piccolo, posso ottenere un leggermente accelerare il metodo con scipy.spatial.distance.cdist e abbinarlo con inner1d, disponendo i dati in modo diverso nella memoria:

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