Come funziona la matrice a distanza condensata? (pdist)

scipy.spatial.distance.pdist restituisce una matrice a distanza condensata. Dalla documentazione :

Restituisce una matrice della distanza condensata Y. Per ogni e (dove), la dist metrica (u = X [i], v = X [j]) viene calcasting e memorizzata nella voce ij.

Pensavo che ij intendesse i*j . Ma penso che potrei sbagliarmi. Prendere in considerazione

 X = array([[1,2], [1,2], [3,4]]) dist_matrix = pdist(X) 

allora la documentazione dice che dist(X[0], X[2]) dovrebbe essere dist_matrix[0*2] . Tuttavia, dist_matrix[0*2] è 0 – non 2.8 come dovrebbe essere.

Qual è la formula che dovrei usare per accedere alla somiglianza di due vettori, dati i e j ?

Puoi guardarlo in questo modo: Supponiamo che x sia m per n. Le possibili coppie di m righe, scelte due alla volta, sono itertools.combinations(range(m), 2) , ad es. Per m=3 :

 >>> import itertools >>> list(combinations(range(3),2)) [(0, 1), (0, 2), (1, 2)] 

Quindi se d = pdist(x) , la k th tuple in combinations(range(m), 2)) fornisce gli indici delle righe di x associate a d[k] .

Esempio:

 >>> x = array([[0,10],[10,10],[20,20]]) >>> pdist(x) array([ 10. , 22.36067977, 14.14213562]) 

Il primo elemento è dist(x[0], x[1]) , il secondo è dist(x[0], x[2]) e il terzo è dist(x[1], x[2]) .

Oppure puoi vederlo come gli elementi nella parte triangular superiore della matrice della distanza quadrata, uniti insieme in una matrice 1D.

Per esempio

 >>> squareform(pdist(x)) array([[ 0. , 10. , 22.361], [ 10. , 0. , 14.142], [ 22.361, 14.142, 0. ]]) >>> y = array([[0,10],[10,10],[20,20],[10,0]]) >>> squareform(pdist(y)) array([[ 0. , 10. , 22.361, 14.142], [ 10. , 0. , 14.142, 10. ], [ 22.361, 14.142, 0. , 22.361], [ 14.142, 10. , 22.361, 0. ]]) >>> pdist(y) array([ 10. , 22.361, 14.142, 14.142, 10. , 22.361]) 

Matrice di distanza condensata a matrice a distanza intera

Una matrice a distanza condensata restituita da pdist può essere convertita in una matrice a distanza intera utilizzando scipy.spatial.distance.squareform :

 >>> import numpy as np >>> from scipy.spatial.distance import pdist, squareform >>> points = np.array([[0,1],[1,1],[3,5], [15, 5]]) >>> dist_condensed = pdist(points) >>> dist_condensed array([ 1. , 5. , 15.5241747 , 4.47213595, 14.56021978, 12. ]) 

Usa la forma squareform per convertire in full matrix:

 >>> dist = squareform(dist_condensed) array([[ 0. , 1. , 5. , 15.5241747 ], [ 1. , 0. , 4.47213595, 14.56021978], [ 5. , 4.47213595, 0. , 12. ], [ 15.5241747 , 14.56021978, 12. , 0. ]]) 

La distanza tra il punto i, j è memorizzata in dist [i, j]:

 >>> dist[2, 0] 5.0 >>> np.linalg.norm(points[2] - points[0]) 5.0 

Indici a indice condensato

Si possono convertire gli indici usati per accedere agli elementi della matrice quadrata all’indice nella matrice condensata:

 def square_to_condensed(i, j, n): assert i != j, "no diagonal elements in condensed matrix" if i < j: i, j = j, i return n*j - j*(j+1)/2 + i - 1 - j 

Esempio:

 >>> square_to_condensed(1, 2, len(points)) 3 >>> dist_condensed[3] 4.4721359549995796 >>> dist[1,2] 4.4721359549995796 

Indice condensato agli indici

Anche l'altra direzione è ansible senza sqaureform, che è migliore in termini di tempo di esecuzione e consumo di memoria:

 import math def calc_row_idx(k, n): return int(math.ceil((1/2.) * (- (-8*k + 4 *n**2 -4*n - 7)**0.5 + 2*n -1) - 1)) def elem_in_i_rows(i, n): return i * (n - 1 - i) + (i*(i + 1))/2 def calc_col_idx(k, i, n): return int(n - elem_in_i_rows(i + 1, n) + k) def condensed_to_square(k, n): i = calc_row_idx(k, n) j = calc_col_idx(k, i, n) return i, j 

Esempio:

 >>> condensed_to_square(3, 4) (1.0, 2.0) 

Confronto di runtime con squareform

 >>> import numpy as np >>> from scipy.spatial.distance import pdist, squareform >>> points = np.random.random((10**4,3)) >>> %timeit dist_condensed = pdist(points) 1 loops, best of 3: 555 ms per loop 

La creazione di sqaureform risulta essere molto lenta:

 >>> dist_condensed = pdist(points) >>> %timeit dist = squareform(dist_condensed) 1 loops, best of 3: 2.25 s per loop 

Se stiamo cercando due punti con la massima distanza non è sorprendente che la ricerca in full matrix sia O (n) mentre in forma condensata solo O (n / 2):

 >>> dist = squareform(dist_condensed) >>> %timeit dist_condensed.argmax() 10 loops, best of 3: 45.2 ms per loop >>> %timeit dist.argmax() 10 loops, best of 3: 93.3 ms per loop 

Ottenere gli elementi necessari per i due punti non richiede quasi tempo in entrambi i casi, ma ovviamente c'è un sovraccarico per il calcolo dell'indice condensato:

 >>> idx_flat = dist.argmax() >>> idx_condensed = dist.argmax() >>> %timeit list(np.unravel_index(idx_flat, dist.shape)) 100000 loops, best of 3: 2.28 µs per loop >>> %timeit condensed_to_square(idx_condensed, len(points)) 100000 loops, best of 3: 14.7 µs per loop 

Il vettore della matrice compressa corrisponde alla regione triangular inferiore della matrice quadrata. Per convertire per un punto in quella regione triangular, è necessario calcolare il numero di punti a sinistra nel triangolo e il numero sopra nella colonna.

È ansible utilizzare la seguente funzione per convertire:

 q = lambda i,j,n: n*j - j*(j+1)/2 + i - 1 - j 

Dai un’occhiata:

 import numpy as np from scipy.spatial.distance import pdist, squareform x = np.random.uniform( size = 100 ).reshape( ( 50, 2 ) ) d = pdist( x ) ds = squareform( d ) for i in xrange( 1, 50 ): for j in xrange( i ): assert ds[ i, j ] == d[ q( i, j, 50 ) ] 

Se vuoi accedere all’elemento di pdist corrispondente pdist (i, j) -th della matrice di distanza quadrata, la matematica è la seguente: Assumi i < j (altrimenti inverti gli indici) se i == j , la risposta è 0.

 X = random((N,m)) dist_matrix = pdist(X) 

Quindi l'elemento (i, j) -th è dist_matrix [ind] dove

 ind = (N - array(range(1,i+1))).sum() + (j - 1 - i). 

Questa è la versione del triangolo superiore ( i ), che deve essere interessante per alcuni:

 condensed_idx = lambda i,j,n: i*n + j - i*(i+1)/2 - i - 1 

Questo è molto facile da capire:

  1. con i*n + j vai nella posizione nella matrice quadrata;
  2. con - i*(i+1)/2 rimuovi il triangolo inferiore (inclusa la diagonale) in tutte le linee prima di i;
  3. con - i rimuovi le posizioni nella riga i prima della diagonale;
  4. con - 1 rimuovi le posizioni nella riga i sulla diagonale.

Dai un’occhiata:

 import scipy from scipy.spatial.distance import pdist, squareform condensed_idx = lambda i,j,n: i*n + j - i*(i+1)/2 - i - 1 n = 50 dim = 2 x = scipy.random.uniform(size = n*dim).reshape((n, dim)) d = pdist(x) ds = squareform(d) for i in xrange(1, n-1): for j in xrange(i+1, n): assert ds[i, j] == d[condensed_idx(i, j, n)] 

Se qualcuno sta cercando una trasformazione inversa (cioè dato un indice di elemento idx , capire a quale elemento (i, j) corrisponde ad esso), ecco una soluzione vettoriale risonabile:

 def actual_indices(idx, n): n_row_elems = np.cumsum(np.arange(1, n)[::-1]) ii = (n_row_elems[:, None] - 1 < idx[None, :]).sum(axis=0) shifts = np.concatenate([[0], n_row_elems]) jj = np.arange(1, n)[ii] + idx - shifts[ii] return ii, jj n = 5 k = 10 idx = np.random.randint(0, n, k) a = pdist(np.random.rand(n, n)) b = squareform(a) ii, jj = actual_indices(idx, n)] assert np.allclose(b[ii, jj, a[idx]) 

L'ho usato per capire gli indici delle righe più vicine in una matrice.

 m = 3 # how many closest lowest_dist_idx = np.argpartition(-a, -m)[-m:] ii, jj = actual_indices(lowest_dist_idx, n) # rows ii[0] and jj[0] are closest