Vectores de desplazamiento por pares entre conjunto de puntos

Tengo una matriz de N puntos en d dimensiones (N, d) y me gustaría hacer una nueva matriz de todos los vectores de desplazamiento para cada par (N choose 2, d) . Si solo quiero la magnitud de estos vectores, puedo usar pdist desde scipy.spatial.distance .

Sería genial si pudiera hacer

 pdist(points, lambda u, v: u - v) 

pero la función metric debe devolver un escalar ( ValueError: setting an array element with a sequence. ) ValueError: setting an array element with a sequence.


Mi solución es usar np.triu_indices :

 i, j = np.triu_indices(len(points), 1) displacements = points[i] - points[j] 

Esto es aproximadamente 20-30 veces más lento que usar pdist (lo comparo tomando la magnitud de los displacements , aunque esta no es una parte que consume mucho tiempo, que supongo que en realidad está haciendo el triángulo superior y ejecutando una indexación elegante).

Sencillo sería

 dis_vectors = [l - r for l, r in itertools.combinations(points, 2)] 

Pero dudo que sea rápido. En realidad, %timeit dice:

Por 3 puntos:

 list : 13 us pdist: 24 us 

Pero ya por 27 puntos:

 list : 798 us pdist: 35.2 us 

¿De cuántos puntos estamos hablando aquí?

Otra posibilidad algo como

 import numpy from operator import mul from fractions import Fraction def binomial_coefficient(n,k): # credit to http://stackoverflow.com/users/226086/nas-banov return int( reduce(mul, (Fraction(ni, i+1) for i in range(k)), 1) ) def pairwise_displacements(a): n = a.shape[0] d = a.shape[1] c = binomial_coefficient(n, 2) out = numpy.zeros( (c, d) ) l = 0 r = l + n - 1 for sl in range(1, n): # no point1 - point1! out[l:r] = a[:n-sl] - a[sl:] l = r r += n - (sl + 1) return out 

Esto simplemente “desliza” la matriz contra sí misma en todas las dimensiones y realiza una resta (de difusión) en cada paso. Tenga en cuenta que no se considera la repetición ni pares iguales (por ejemplo, punto1 – punto1).

Esta función aún se desempeña bien en el rango de 1000 puntos con 31.3ms , mientras que pdist es aún más rápido con 20.7 ms y la lista de comprensión ocupa el tercer lugar con 1.23 s .

Si calcula el producto cartesiano completo de las diferencias, aplana la matriz 2D resultante y crea sus propios índices para extraer el triángulo superior, puede hacer que sea “solo” 6 pdist más lento que pdist :

 In [39]: points = np.random.rand(1000, 2) In [40]: %timeit pdist(points) 100 loops, best of 3: 5.81 ms per loop In [41]: %%timeit ...: n = len(points) ...: rng = np.arange(1, n) ...: idx = np.arange(n *(n-1) // 2) + np.repeat(np.cumsum(rng), rng[::-1]) ...: np.take((points[:, None] - points).reshape(-1, 2), idx, axis=0) ...: 10 loops, best of 3: 33.9 ms per loop 

También puede acelerar su solución, crear los índices usted mismo y usar la toma en lugar de la indexación elegante:

 In [75]: %%timeit ...: n = len(points) ...: rng = np.arange(1, n) ...: idx1 = np.repeat(rng - 1, rng[::-1]) ...: idx2 = np.arange(n*(n-1)//2) + np.repeat(n - np.cumsum(rng[::-1]), rng[::-1]) ...: np.take(points, idx1, axis=0) - np.take(points, idx2, axis=0) ...: 10 loops, best of 3: 38.8 ms per loop