numpy - instalar - ¿Cómo funciona la matriz de distancia condensada?(pdist)

pip install numpy scipy matplotlib pandas

Question

scipy.spatial.distance.pdist devuelve una matriz de distancia condensada. De la documentación :

Devuelve una matriz de distancia condensada Y. Para cada uno y (donde), la métrica dist (u = X [i], v = X [j]) se calcula y almacena en la entrada ij.

Pensé que significaba i*j . Pero creo que podría estar equivocado. Considerar

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

entonces la documentación dice que dist(X[0], X[2]) debe ser dist_matrix[0*2] . Sin embargo, dist_matrix[0*2] es 0 - no 2.8 como debería ser.

¿Cuál es la fórmula que debería usar para acceder a la similitud de dos vectores, dado i y j ?

Answer 1

Matriz de distancia condensada a matriz de distancia completa

Una matriz de distancia condensada como la devuelta por pdist se puede convertir a una matriz de distancia completa usando 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.        ])

Use forma squareform para convertir a matriz completa:

>>> 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 distancia entre el punto i, j se almacena en dist [i, j]:

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

Índices al índice condensado

Se pueden convertir los índices utilizados para acceder a los elementos de la matriz cuadrada al índice en la matriz condensada:

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

Ejemplo:

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

Índice condensado a índices

También la otra dirección es posible sin sqaureform, que es mejor en términos de tiempo de ejecución y consumo de 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

Ejemplo:

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

Comparación en tiempo de ejecución 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

Crear el sqaureform resulta ser muy lento:

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

Si buscamos dos puntos con una distancia máxima, no es sorprendente que la búsqueda en matriz completa sea O (n) mientras que en forma condensada 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

Obtener las oportunidades para los dos puntos casi no lleva tiempo en ambos casos, pero por supuesto hay una sobrecarga para calcular el índice condensado:

>>> 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

Answer 2

El vector de la matriz comprimida corresponde a la región triangular inferior de la matriz cuadrada. Para convertir un punto en esa región triangular, debe calcular el número de puntos a la izquierda en el triángulo y el número de arriba en la columna.

Puede usar la siguiente función para convertir:

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

Comprobar:

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 ) ]

Answer 3

Esta es la versión del triángulo superior ( i <j ), que debe ser interesante para algunos:

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

Esto es muy fácil de entender:

1. con i*n + j vas a la posición en la matriz formada por cuadrados;
2. con - i*(i+1)/2 elimina el triángulo inferior (incluida la diagonal) en todas las líneas antes de i;
3. con - i quita las posiciones en la línea i antes de la diagonal;
4. con - 1 quita las posiciones en la línea i en la diagonal.

Comprobar:

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)]

Answer 4

Puedes verlo de esta manera: supongamos que x es m por n. Los posibles pares de m filas, elegidos dos a la vez, son itertools.combinations(range(m), 2) , por ejemplo, para m=3 :

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

Entonces, si d = pdist(x) , la k tupla en combinations(range(m), 2)) da los índices de las filas de x asociadas con d[k] .

Ejemplo:

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

El primer elemento es dist(x[0], x[1]) , el segundo es dist(x[0], x[2]) y el tercero es dist(x[1], x[2]) .

O puede verlo como los elementos en la parte triangular superior de la matriz de distancia cuadrada, agrupados en una matriz 1D.

P.ej

>>> 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])

Answer 5

Si alguien está buscando una transformación inversa (es decir, dado un índice de elemento idx , averigüe qué elemento (i, j) corresponde), aquí hay una solución de vectores resonables:

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])

Lo usé para calcular los índices de las filas más cercanas en una matriz.

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

Answer 6

Si desea acceder al elemento de pdist correspondiente al elemento (i, j) -th de la matriz de distancia cuadrada, la matemática es la siguiente: Suponga que i < j (de lo contrario, invierta los índices) si i == j , la respuesta es 0.

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

Entonces el elemento (i, j) es dist_matrix [ind] donde

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