紧凑距离矩阵是如何工作的？(pdist)

Question

scipy.spatial.distance.pdist返回一个压缩的距离矩阵。来自文档:

返回一个压缩的距离矩阵Y。对于每对i和j（其中i

我以为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，我应该使用什么公式来计算两个向量的相似度？

- Rafael Almeida

Answer 1

您可以这样看待它：假设x是m乘n的矩阵。从中选择两行的可能组合为itertools.combinations(range(m), 2)，例如当m=3时：

您可以这样看待它：假设x是m乘n的矩阵。从中选择两行的可能组合为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])

- Warren Weckesser

---

我明白了，很有趣。看起来squareform更容易使用。sq_form[i,j]将准确地给出y[i]和y[j]之间的距离。然而，我认为压缩形式在内存方面更好。也许我应该多了解一下squareform的功能。但是，难道没有一个简单的公式可以将i,j转换为dist位置吗？ - Rafael Almeida

这真的是有文档记录的行为吗？它有道理，当然，但API中没有任何东西表明它应该与“combinations(range(m), 2))”相结合，这对应于距离矩阵的下三角形。为什么不是上三角形呢？ - VF1

你为什么说它对应于下三角？例如，list(combinations(range(4), 2))会得到[(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)]。该列表中的每个元组都具有形式（行索引，列索引），因此它对应于上三角。 - Warren Weckesser

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

压缩索引指标

不使用方阵也可以进行另一种方向的操作，这样在运行时间和内存消耗方面更好：

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

生成方阵形式的过程非常缓慢：

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

- lumbric

---

我有一个方形的距离矩阵 - 有没有一种函数可以将其转换为压缩形式？（即相反于squareform）...像linkage这样的函数期望压缩形式... 编辑 squareform函数可以双向转换...太酷了..."... 反之亦然。" 编辑2 我称之为“方形”，应该是“冗余”的 编辑3 并且linkage将与两种形式一起使用... - Nate Anderson

@The Red Pea squareform 是其自身的逆运算（即当在完整距离矩阵上运行时，它将其转换为方阵形式） - moustachio

你是怎么想出 calc_row_idx 和 calc_col_idx 的？ - CMCDragonkai

如果我没有记错的话，@CMCDragonkai 的想法是从 square_to_condensed() 函数中取出公式，并使用二次方程式解决 i 或 j。据我所记，这个计算有点繁琐，但只涉及基本的代数和一些推理，其中某些解是不可能的（因为它是负数或复数）。画一个例子三角形矩阵并清空对角线会有帮助。 - lumbric

在Python 3中，“压缩索引的索引”应使用整数除法“//”，而不是浮点除法“/”，以使输出为整数。 - mic

@mic 谢谢你指出来！我已经更新了答案。我必须承认，在2016年的工作中，我仍在使用Python 2.x。 - lumbric

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

- shaunc

---

请注意，这确实是底部三角形区域，这可能对某些人来说很奇怪。 - shaunc

底部三角形是上部三角形的转置，因为距离矩阵是对称的，即交换j,i -> i,j会得到相同的结果。您的解决方案使用了下三角形的解释，但上三角形版本也没有任何错误（我认为这是人们更常见的思考方式）。 - David Marx

我正在尝试朝相反的方向移动：给定压缩距离矩阵（即平面向量）中的索引，如何在不强制转换为方形形式的情况下获取对应于该值的矩阵索引（i，j）？ - David Marx

squareform按行序列化距离向量。这意味着较小的索引始终是第一个：(0,1),(0,2),(0,3).... - user1261273

Answer 4

我也有同样的问题。我发现使用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))

- Rustam Guliev

Answer 5

这是 上三角 版本 (i < j)，对一些人来说可能很有趣：

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

这很容易理解：

1. 使用 i*n + j 可以定位矩阵中的位置；
2. 使用 - i*(i+1)/2 可以删除第 i 行及其上方（包括对角线）的下三角区域；
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)]

- HongboZhu

Answer 6

如果有人正在寻找反向转换（即给定一个元素索引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

- Ben Usman

Answer 7

如果想访问与平方距离矩阵的(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).

- user1462620

---

请注意，array(range(1,i+1))).sum() == ((i+1)*i)/2 这段代码是关于高斯算法的。 - lumbric

如果我有一个数组 array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])，其中 N = 5，在第 i = 3 行和 j = 2 列的公式 (5 - np.array(range(1,3+1))).sum() + (2 - 1 - 3) 应该给我 5，但实际上它给了我 7。 - Veronica