Как работает уплотненная матрица расстояний? (Pdist)
scipy.spatial.distance.pdist
возвращает сжатую матрицу расстояний. Из документация:
Возвращает матрицу сокращенных расстояний Y. Для каждого и (где) метрика dist (u = X [i], v = X [j]) вычисляется и сохраняется в записи ij.
Я думал, что ij
означает i*j
. Но я думаю, что я ошибаюсь. Рассмотрим
X = array([[1,2], [1,2], [3,4]])
dist_matrix = pdist(X)
то в документации указано, что dist(X[0], X[2])
должен быть dist_matrix[0*2]
. Однако dist_matrix[0*2]
равно 0, а не 2.8, как и должно быть.
Какую формулу я должен использовать для доступа к подобию двух векторов, заданных i
и j
?
Ответы
Ответ 1
Вы можете посмотреть на это следующим образом: предположим, что x
есть m по n. Возможные пары строк m
, выбранных по два за раз, составляют itertools.combinations(range(m), 2)
, например, для m=3
:
>>> import itertools
>>> list(combinations(range(3),2))
[(0, 1), (0, 2), (1, 2)]
Итак, если d = pdist(x)
, k
th кортеж в combinations(range(m), 2))
дает индексы строк x
, связанных с d[k]
.
Пример:
>>> x = array([[0,10],[10,10],[20,20]])
>>> pdist(x)
array([ 10. , 22.36067977, 14.14213562])
Первый элемент dist(x[0], x[1])
, второй - dist(x[0], x[2])
, а третий - dist(x[1], x[2])
.
Или вы можете просмотреть его как элементы в верхней треугольной части квадратной матрицы расстояния, нанизанные в 1D-массив.
например.
>>> 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])
Ответ 2
Вектор сжатой матрицы соответствует нижней треугольной области квадратной матрицы. Чтобы преобразовать точку в этой треугольной области, вам нужно вычислить количество точек слева в треугольнике и число, указанное в столбце.
Вы можете использовать следующую функцию для преобразования:
q = lambda i,j,n: n*j - j*(j+1)/2 + i - 1 - j
Check:
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 ) ]
Ответ 3
Матрица с конденсированным расстоянием до матрицы полного расстояния
Сжатая матрица расстояний, возвращаемая pdist, может быть преобразована в матрицу полного расстояния с помощью 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. ])
Используйте squareform
для преобразования в полную матрицу:
>>> 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. ]])
Расстояние между точкой i, j сохраняется в dist [i, j]:
>>> dist[2, 0]
5.0
>>> np.linalg.norm(points[2] - points[0])
5.0
Индексы к конденсированному индексу
Можно преобразовать индексы, используемые для доступа к элементам квадратной матрицы, к индексу в конденсированной матрице:
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
Пример:
>>> square_to_condensed(1, 2, len(points))
3
>>> dist_condensed[3]
4.4721359549995796
>>> dist[1,2]
4.4721359549995796
Сжатый индекс индексов
Также возможно другое направление без sqaureform, которое лучше с точки зрения времени выполнения и потребления памяти:
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
Пример:
>>> condensed_to_square(3, 4)
(1.0, 2.0)
Сравнение времени выполнения с квадратом
>>> 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
Создание sqaureform оказывается очень медленным:
>>> dist_condensed = pdist(points)
>>> %timeit dist = squareform(dist_condensed)
1 loops, best of 3: 2.25 s per loop
Если мы ищем две точки с максимальным расстоянием, неудивительно, что поиск в полной матрице равен O (n), а в конденсированной форме только 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
Получение inideces для двух точек практически не требует времени в обоих случаях, но, конечно, для вычисления сконденсированного индекса есть некоторые накладные расходы:
>>> 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
Ответ 4
Если вы хотите получить доступ к элементу pdist
, соответствующему (i, j) -тому элементу матрицы квадратного расстояния, математика выглядит следующим образом: Предположим i < j
(иначе индексы переворота), если i == j
, ответ равен 0.
X = random((N,m))
dist_matrix = pdist(X)
Тогда (i, j) -й элемент является dist_matrix [ind], где
ind = (N - array(range(1,i+1))).sum() + (j - 1 - i).
Ответ 5
Это версия верхнего треугольника (я < j), которая должна быть интересна некоторым:
condensed_idx = lambda i,j,n: i*n + j - i*(i+1)/2 - i - 1
Это очень легко понять:
- с
i*n + j
вы переходите к позиции в квадратной матрице;
- с
- i*(i+1)/2
вы удаляете нижний треугольник (включая диагональ) во всех строках до i;
- с
- i
вы удаляете позиции в строке я перед диагональю;
- с
- 1
вы удаляете позиции в строке я по диагонали.
Check:
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)]