Construcción paralela de una matriz de distancia

Question

Trabajo en agrupamiento aglomerativo jerárquico en grandes cantidades de vectores multidimensionales, y noté que el mayor cuello de botella es la construcción de la matriz de distancia. Una implementación ingenua para esta tarea es la siguiente (en este caso en Python):

''' v = an array (N,d), where rows are the observations 
and columns the dimensions''' 
def create_dist_matrix(v): 
    N = v.shape[0] 
    D = np.zeros((N,N)) 
    for i in range(N): 
     for j in range(i+1): 
      D[i,j] = cosine(v[i,:],v[j,:]) # scipy.spatial.distance.cosine() 
    return D 

Me preguntaba cuál es la mejor manera de añadir algo de paralelismo a esta rutina. Una manera fácil sería romper y asignar el bucle externo a una cantidad de trabajos, p. si tiene 10 procesadores, cree 10 trabajos diferentes para diferentes rangos de i y luego concatené los resultados. Sin embargo, esta solución "horizontal" no parece del todo correcta. ¿Hay algún otro algoritmo paralelo (o librerías existentes) para esta tarea? Cualquier ayuda sería muy apreciada.

Answer 1

Dudo que lo obtenga más rápido que pdist en el módulo scipy. Probablemente es por esto que dice

Tenga en cuenta que debe evitar pasar una referencia a una de las funciones de distancia definidas en esta biblioteca. Por ejemplo ,:

dm = pdist(X, sokalsneath) 

sería calcular las distancias por pares entre los vectores en X utilizando la función sokalsneath Python. Esto daría como resultado que sokalsneath se llame n elija 2 veces, que es ineficaz. En cambio, la versión optimizada C es más eficiente, y lo llamamos con la siguiente sintaxis .:

dm = pdist(X, 'sokalsneath') 

Así que no se utiliza ninguna función de Python, si utiliza pdist(X, 'cosine'). Cuando lo ejecuto, para mí, parece que solo usa un núcleo, por lo que si tienes muchos núcleos, puedes obtenerlo más rápido. Pero tenga en cuenta que para lograr esto, su implementación nativa debe ser tan rápida como la de SciPy. Eso no será trivial. Prefieres ser paciente o elegir un método de agrupación diferente, e. gramo. un algoritmo que admite un índice espacial.

Answer 2

Parece que scikit-learn tiene una versión paralela de pdist llamada pairwise_distances

from sklearn.metrics.pairwise import pairwise_distances 

D = pairwise_distances(X = v, metric = 'cosine', n_jobs = -1) 

donde n_jobs = -1 especifica que se utilizarán todas las CPU.

Answer 3

Ver @agartland responder — puede especificar n_jobs en sklearn.metrics.pairwise.pairwise_distances o buscar algoritmo de agrupamiento en sklearn.cluster con n_jobs parámetro. P.ej. sklearn.cluster.KMeans.

Aún así, si se siente aventurero, puede implementar su propio cálculo. Por ejemplo, si necesita matriz de distancia 1D para scipy.cluster.hierarchy.linkage puede utilizar:

#!/usr/bin/env python3 
from multiprocessing import Pool 
import numpy as np 
from time import time as ts 


data = np.zeros((100,10)) # YOUR data: np.array[n_samples x m_features] 
n_processes = 4   # YOUR number of processors 
def metric(a, b):   # YOUR dist function 
    return np.sum(np.abs(a-b)) 


n = data.shape[0] 
k_max = n * (n - 1) // 2 # maximum elements in 1D dist array 
k_step = n ** 2 // 500 # ~500 bulks 
dist = np.zeros(k_max) # resulting 1D dist array 


def proc(start): 
    dist = [] 
    k1 = start 
    k2 = min(start + k_step, k_max) 
    for k in range(k1, k2): 
     # get (i, j) for 2D distance matrix knowing (k) for 1D distance matrix 
     i = int(n - 2 - int(np.sqrt(-8 * k + 4 * n * (n - 1) - 7)/2.0 - 0.5)) 
     j = int(k + i + 1 - n * (n - 1)/2 + (n - i) * ((n - i) - 1)/2) 
     # store distance 
     a = data[i, :] 
     b = data[j, :] 
     d = metric(a, b) 
     dist.append(d) 
    return k1, k2, dist 


ts_start = ts() 
with Pool(n_processes) as pool: 
    for k1, k2, res in pool.imap_unordered(proc, range(0, k_max, k_step)): 
     dist[k1:k2] = res 
     print("{:.0f} minutes, {:,}..{:,} out of {:,}".format(
      (ts() - ts_start)/60, k1, k2, k_max)) 


print("Elapsed %.0f minutes" % ((ts() - ts_start)/60)) 
print("Saving...") 
np.savez("dist.npz", dist=dist) 
print("DONE") 

Para que lo sepas, scipy.cluster.hierarchy.linkage aplicación no es paralelo y su complejidad es al menos O (N * N). No estoy seguro de si scipy tiene una implementación paralela de esta función.

Answer 4

Si decide orquestar el multiprocesamiento usted mismo, es posible que desee dividir el número de cálculos de forma pareja entre las CPU para acortar al máximo los cálculos. Entonces las respuestas a this question on equally splitting the diagonal matrix pueden ser útiles.