凝聚距离矩阵是如何工作的？(pdist)

Question

scipy.spatial.distance.pdist返回一个压缩距离矩阵。来自文献资料

为每个和(其中)返回一个凝聚距离矩阵Y，计算度量dist(u=Xi，v=Xj)并将其存储在条目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，我应该使用什么公式来访问两个向量的相似性

Answer 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)，则combinations(range(m), 2))中的k第四元组将给出与d[k]关联的x行的索引。

示例：

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

或者你可以把它看作是方形距离矩阵的上三角部分中的元素，把它们串成一个一维数组。

例如。

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

凝聚距离矩阵到全距离矩阵

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之间的距离存储在disti，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

缩合指数

另外，如果没有memory，就运行时和内存消耗而言，就更好的方面而言，另一个方向是可能的：

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

创建广场改革的过程非常缓慢：

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

Answer 3

压缩矩阵的向量对应于平方矩阵的下三角区域。要转换该三角区域中的一个点，需要计算三角形中左边的点数和列中的上面的点数。

您可以使用以下函数进行转换：

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

检查：

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