82

scipy.spatial.distance.pdist返回一个压缩的距离矩阵。从文档中

返回压缩距离矩阵 Y。对于每个 and (其中 ),度量 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

4

7 回答 7

105

你可以这样看:假设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元组给出与相关联combinations(range(m), 2))的行的索引。xd[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])

或者您可以将其视为平方距离矩阵的上三角部分中的元素,串在一起形成一维数组。

例如

>>> 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])
于 2012-10-26T02:01:22.267 回答
48

压缩距离矩阵到全距离矩阵

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)

与 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

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

在这两种情况下,获取这两个点的索引几乎都不需要时间,但是计算压缩索引当然会有一些开销:

>>> 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
于 2016-04-26T14:11:58.573 回答
18

压缩矩阵的向量对应于方阵的底部三角形区域。要转换为该三角形区域中的点,您需要计算三角形左侧的点数以及列中上方的点数。

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

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 ) ]
于 2013-10-25T04:42:47.523 回答
8

我有同样的问题。而且我发现使用起来更简单numpy.triu_indices

import numpy as np
from scipy.spatial.distance import pdist, squareform
N = 10

# Calculate distances
X = np.random.random((N,3))
dist_condensed = pdist(X)

# Get indexes: matrix indices of dist_condensed[i] are [a[i],b[i]]
a,b = np.triu_indices(N,k=1)

# Fill distance matrix
dist_matrix = np.zeros((N,N))
for i in range(len(dist_condensed)):
    dist_matrix[a[i],b[i]] = dist_condensed[i]
    dist_matrix[b[i],a[i]] = dist_condensed[i]

# Compare with squareform output
np.all(dist_matrix == squareform(dist_condensed))
于 2019-04-03T11:10:05.300 回答
6

这是上三角形版本(i < j),这对某些人来说一定很有趣:

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

这很容易理解:

  1. i*n + j你一起去方阵中的位置;
  2. 在 i 之前的所有行中- i*(i+1)/2删除下三角形(包括对角线);
  3. - i您一起删除对角线前 i 行中的位置;
  4. - 1您一起删除对角线上第 i 行的位置。

查看:

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)]
于 2014-09-03T16:48:31.657 回答
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).
于 2013-07-26T00:21:27.437 回答
3

如果有人正在寻找逆变换(即给定一个元素索引idx,找出(i, j)对应的元素),这是一个合理的矢量解决方案:

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

我用它来计算矩阵中最近行的索引。

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
于 2018-03-01T21:59:15.710 回答