Wie man die verdichtete Form der paarweisen Abstände direkt erhält?

Question

Ich habe eine sehr große scipy spärliche CSR-Matrix. Es ist eine 100.000x2.000.000 dimensionale Matrix. Nennen wir es X. Jede Zeile ist ein Beispielvektor in einem 2.000.000 dimensionalen Raum.

Ich muss die Kosinusabstände zwischen jedem Probenpaar sehr effizient berechnen. Ich habe sklearn pairwise_distances Funktion mit einer Untermenge von Vektoren in X verwendet, die mir eine dichte Matrix D gibt: die quadratische Form der paarweisen Abstände, die redundante Einträge enthält. Wie kann ich sklearn pairwise_distances verwenden, um die komprimierte Form direkt zu erhalten? Bitte beachten Sie http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.pdist.html, um zu sehen, was die kondensierte Form ist. Es ist die Ausgabe von scipy pdist Funktion.

Ich habe Speicher Einschränkungen und ich kann nicht berechnen die quadratische Form und dann die verdichtete Form. Aufgrund von Speicherbeschränkungen kann ich scipy pdist auch nicht verwenden, da es eine dichte Matrix X benötigt, die nicht wieder in den Speicher passt. Ich dachte daran, verschiedene Blöcke von X zu durchlaufen und die verdichtete Form für jeden Block zu berechnen und sie zu verbinden, um die vollständige kondensierte Form zu erhalten, aber das ist relativ umständlich. Irgendwelche besseren Ideen?

Jede Hilfe wird sehr geschätzt. Danke im Voraus.

Unten ist ein reproduzierbares Beispiel (natürlich zu Demonstrationszwecken X ist viel kleiner):

from scipy.sparse import rand 
from scipy.spatial.distance import pdist 
from sklearn.metrics.pairwise import pairwise_distances 
X = rand(1000, 10000, density=0.01, format='csr') 
dist1 = pairwise_distances(X, metric='cosine') 
dist2 = pdist(X.A, 'cosine') 

Wie Sie sehen dist2 in kondensierter Form und ist ein 499.500-dimensionaler Vektor. Aber dist1 ist in der symmetrischen Quadratform und ist eine 1000x1000 Matrix.

Answer 1

Ich grub in den Code für beide Versionen und denke, ich verstehe, was beide tun.

Beginnen Sie mit einem kleinen einfachen X (dicht):

X = np.arange(9.).reshape(3,3) 

pdist Cosinus tut:

norms = _row_norms(X) 
_distance_wrap.pdist_cosine_wrap(_convert_to_double(X), dm, norms) 

wo _row_norms eine Reihe Punkt ist - mit einsum:

norms = np.sqrt(np.einsum('ij,ij->i', X,X) 

So Dies ist der erste Ort, wo X muss ein Array sein.

Ich habe nicht in die cosine_wrap gegraben, aber es scheint (wahrscheinlich in cython) zu tun

xy = np.dot(X, X.T) 
# or xy = np.einsum('ij,kj',X,X) 

d = np.zeros((3,3),float) # square receiver 
d2 = []      # condensed receiver 
for i in range(3): 
    for j in range(i+1,3): 
     val=1-xy[i,j]/(norms[i]*norms[j]) 
     d2.append(val) 
     d[j,i]=d[i,j]=val 

print('array') 
print(d) 
print('condensed',np.array(d2)) 

from scipy.spatial import distance 
d1=distance.pdist(X,'cosine') 
print(' pdist',d1) 

Herstellung:

array 
[[ 0.   0.11456226 0.1573452 ] 
[ 0.11456226 0.   0.00363075] 
[ 0.1573452 0.00363075 0.  ]] 

condensed [ 0.11456226 0.1573452 0.00363075] 
    pdist [ 0.11456226 0.1573452 0.00363075] 

distance.squareform(d1) produziert das gleiche wie meine d Array.

kann ich die gleiche quadratische Anordnung erzeugen, indem die xy Skalarprodukt mit dem entsprechenden norm äußere Produkt Dividieren:

dd=1-xy/(norms[:,None]*norms) 
dd[range(dd.shape[0]),range(dd.shape[1])]=0 # clean up 0s 

Oder durch X vor der Einnahme Skalarprodukt normalisieren.Dies scheint die Version scikit zu sein.

Xnorm = X/norms[:,None] 
1-np.einsum('ij,kj',Xnorm,Xnorm) 

scikit einige cython Code zu tun schneller spärliche Berechnungen (über diejenigen, die von sparse.sparse, aber mit dem gleichen csr Format) hinzugefügt:

from scipy import sparse 
Xc=sparse.csr_matrix(X) 

# csr_row_norm - pyx of following 
cnorm = Xc.multiply(Xc).sum(axis=1) 
cnorm = np.sqrt(cnorm) 
X1 = Xc.multiply(1/cnorm) # dense matrix 
dd = 1-X1*X1.T 

Um eine schnelle kondensierter Form mit spärlichen Matrizen bekomme ich denke, dass Sie eine schnelle verkürzte Version von X1*X1.T implementieren müssen. Das bedeutet, dass Sie verstehen müssen, wie die Multiplikation der Sparse-Matrix implementiert wird - unter c. Der scikit Cython 'Fast Sparse' Code könnte auch Ideen geben.

numpy hat einige tri... Funktionen, die einfach Python-Code sind. Es wird nicht versucht, Zeit oder Platz zu sparen, indem Tri-Berechnungen direkt implementiert werden. Es ist einfacher, über das rechteckige Layout eines nd-Arrays (mit Form und Schritten) zu iterieren als die komplexeren Schritte mit variabler Länge eines dreieckigen Arrays. Die verdichtete Form schneidet den Raum und die Berechnungsschritte nur um die Hälfte.

============

Hier ist der Hauptteil der cpdist_cosine Funktion, die iteriert über i und die oberen j, dot(x[i],y[j])/(norm[i]*norm[j]) berechnen.

for (i = 0; i < m; i++) { 
    for (j = i + 1; j < m; j++, dm++) { 
     u = X + (n * i); 
     v = X + (n * j); 
     cosine = dot_product(u, v, n)/(norms[i] * norms[j]); 
     if (fabs(cosine) > 1.) { 
      /* Clip to correct rounding error. */ 
      cosine = npy_copysign(1, cosine); 
     } 
     *dm = 1. - cosine; 
    } 
} 

https://github.com/scipy/scipy/blob/master/scipy/spatial/src/distance_impl.h