Plus rapide par paires distance métrique en python

J'ai un 1D tableau de nombres, et que vous voulez calculer toutes les distances euclidiennes. J'ai une méthode de (grâce à) de le faire avec de la radiodiffusion, mais il est inefficace parce qu'il calcule la distance deux fois. Et il n'est pas à l'échelle.

Voici un exemple qui me donne ce que je veux avec un tableau de 1000 numéros.

import numpy as np
import random
r = np.array([random.randrange(1, 1000) for _ in range(0, 1000)])
dists = np.abs(r - r[:, None])

Quelle est la manière la plus rapide de mise en œuvre dans scipy/numpy/scikit-learn que je peux utiliser, pour ce faire, étant donné qu'elle a pour s'adapter à des situations où les 1D tableau a >10k valeurs.

Remarque: la matrice est symétrique, donc je suppose qu'il est possible d'obtenir au moins un 2x speedup par l'aborder, je ne sais juste pas comment.

Il y a une fonction pour ça: scipy.spatial.distance.pdist. Je ne sais pas si c'est l'option la plus rapide, car il a besoin d'avoir des contrôles pour les données multidimensionnelles, non-Euclidienne normes, et d'autres choses, mais il est intégré.
Combien de temps il vous faut pour cela être? Il ne va jamais à l'échelle de mieux que O(n^2), puisque vous avez à remplir n^2 entrées de sortie. Votre solution existante est O(n^2), et il ne semble pas beaucoup de place pour les principales optimisations.
Cela semble échelle >10k valeurs assez bien déjà quand je l'ai essayer. N'oubliez pas que vous devez remplir à 100 millions d'entrées de sortie. C'est presque la moitié d'un gigaoctet de paires distances.
Je ne le pense pas... Si vous suivez le code source, à la fin, ceci est la fonction qui est appelé. Non seulement il n'y a pas de fantaisie, d'optimisation, mais pour 1D vecteurs c'est la quadrature et en prenant la racine carrée de calculer la valeur absolue. Probablement pire que l'OP du code pour son cas d'utilisation particulier.
Si je ne me trompe pas, scipy est toujours compilé avec BLAS, c'est pas une option comme avec numpy.

OriginalL'auteur roblanf | 2013-11-29

17

Ni l'un ni l'autre des réponses tout à fait répondu à la question 1 a été en Cython, on a été plus lent. Mais à la fois très utiles conseils. Leur mise en place suggère que scipy.spatial.distance.pdist est le chemin à parcourir.

Voici un code:
```
import numpy as np
import random
import sklearn.metrics.pairwise
import scipy.spatial.distance

r = np.array([random.randrange(1, 1000) for _ in range(0, 1000)])
c = r[:, None]

def option1(r):
    dists = np.abs(r - r[:, None])

def option2(r):
    dists = scipy.spatial.distance.pdist(r, 'cityblock')

def option3(r):
    dists = sklearn.metrics.pairwise.manhattan_distances(r)
```
De la synchronisation avec IPython:
```
In [36]: timeit option1(r)
100 loops, best of 3: 5.31 ms per loop

In [37]: timeit option2(c)
1000 loops, best of 3: 1.84 ms per loop

In [38]: timeit option3(c)
100 loops, best of 3: 11.5 ms per loop
```
Je n'ai pas essayé la Cython de mise en œuvre (je ne peux pas l'utiliser pour ce projet), mais comparer mes résultats à l'autre réponse qui l'a fait, il ressemble à scipy.spatial.distance.pdist est environ un tiers plus lent que le Cython de mise en œuvre (en prenant en compte les différentes machines par l'analyse comparative sur la np.abs solution).

Je suppose que c'est aussi rapide comme: scikit-learn.org/stable/modules/generated/... ? la version en sklearn?

OriginalL'auteur roblanf

Ici est un Cython de mise en œuvre qui donne plus de 3X amélioration de la vitesse pour cet exemple sur mon ordinateur. Ce calendrier doit être examiné pour les grands tableaux difficile, parce que les routines BLAS pouvez probablement échelle beaucoup mieux que cela, plutôt naïf code.

Je sais que vous avez demandé quelque chose à l'intérieur de scipy/numpy/scikit-learn, mais peut-être cela va ouvrir de nouvelles possibilités pour vous:

Fichier my_cython.pyx:

import numpy as np
cimport numpy as np
import cython

cdef extern from "math.h":
    double abs(double t)

@cython.wraparound(False)
@cython.boundscheck(False)
def pairwise_distance(np.ndarray[np.double_t, ndim=1] r):
    cdef int i, j, c, size
    cdef np.ndarray[np.double_t, ndim=1] ans
    size = sum(range(1, r.shape[0]+1))
    ans = np.empty(size, dtype=r.dtype)
    c = -1
    for i in range(r.shape[0]):
        for j in range(i, r.shape[0]):
            c += 1
            ans[c] = abs(r[i] - r[j])
    return ans

La réponse est un 1-D tableau contenant tous les non-répétée évaluations.

À l'importation dans Python:

import numpy as np
import random

import pyximport; pyximport.install()
from my_cython import pairwise_distance

r = np.array([random.randrange(1, 1000) for _ in range(0, 1000)], dtype=float)

def solOP(r):
    return np.abs(r - r[:, None])

De la synchronisation avec IPython:

In [2]: timeit solOP(r)
100 loops, best of 3: 7.38 ms per loop

In [3]: timeit pairwise_distance(r)
1000 loops, best of 3: 1.77 ms per loop

Sûrement que vous avez voulu dire fabs -- abs est int variante.

OriginalL'auteur Saullo G. P. Castro

3

À l'aide de la moitié de la mémoire, mais 6 fois plus lent que np.abs(r - r[:, None]):
```
triu = np.triu_indices(r.shape[0],1)
dists2 = abs(r[triu[1]]-r[triu[0]])
```
OriginalL'auteur cyborg

Vous devez vous connecter pour publier un commentaire.