如何直接得到成对距离的压缩形式？

Posted 2023-03-12

【中文标题】如何直接得到成对距离的压缩形式？

【英文标题】：How to get the condensed form of pairwise distances directly?

【发布时间】：2016-11-14 23:00:01

【问题描述】：

我有一个非常大的 scipy 稀疏 csr 矩阵。它是一个 100,000x2,000,000 维矩阵。我们称之为X。每行是 2,000,000 维空间中的一个样本向量。

我需要非常有效地计算每对样本之间的余弦距离。我一直在使用sklearn pairwise_distances 函数和X 中的向量子集，这给了我一个密集矩阵D：包含冗余​​条目的成对距离的平方形式。如何使用sklearn pairwise_distances 直接获取压缩形式？请参考http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.pdist.html 以了解压缩形式是什么。是scipy pdist函数的输出。

我有内存限制，我无法计算平方形式，然后得到压缩形式。由于内存限制，我也不能使用scipy pdist，因为它需要一个密集矩阵X，它不再适合内存。我想过循环遍历X的不同块并计算每个块的压缩形式并将它们连接在一起以获得完整的压缩形式，但这相对繁琐。有更好的想法吗？

非常感谢任何帮助。提前致谢。

下面是一个可重现的例子（当然为了演示目的X 要小得多）：

from scipy.sparse import rand
from scipy.spatial.distance import pdist
from sklearn.metrics.pairwise import pairwise_distances
X = rand(1000, 10000, density=0.01, format='csr')
dist1 = pairwise_distances(X, metric='cosine')
dist2 = pdist(X.A, 'cosine')

如您所见，dist2 是压缩形式，是一个 499500 维向量。但是dist1 是对称的正方形，是一个 1000x1000 的矩阵。

【问题讨论】：

你需要添加一个具体的例子；我们可以复制粘贴并运行的东西。显然它不会遇到内存问题。但是除非我们正在处理完全相同的问题，否则很难理解您的口头描述。我很了解稀疏矩阵代码，但没有使用sklearn。所以像“浓缩形式”这样的术语是外来的。

@hpaulj 似乎一切都在***上被问到，最终：***.com/questions/13079563/…

还有关于从值向量中填充上/下三角形（或两者）的问题。

在 scikit-learn 和稀疏距离上搜索会变成 ***.com/q/8956274/901925

@hpaulj：当然。我添加了一个示例，以及对压缩形式的一些参考。压缩形式是线性代数中使用的通用术语。对于处理具有特定结构的大型矩阵的算法，使用代数运算以压缩形式表示矩阵通常效率更高。有多种压缩形式表示，一些使用特征值/特征向量产生块对角矩阵。在这里，因为成对距离矩阵是对称的，所以最简单的压缩形式仅由其上（或下）三角形元素组成。

【参考方案1】：

我研究了这两个版本的代码，并认为我理解了这两个版本的作用。

从一个简单的小X（密集）开始：

X = np.arange(9.).reshape(3,3)

pdist 余弦：

norms = _row_norms(X)
_distance_wrap.pdist_cosine_wrap(_convert_to_double(X), dm, norms)

_row_norms 是一个行点 - 使用 einsum:

norms = np.sqrt(np.einsum('ij,ij->i', X,X)

所以这是X 必须是数组的第一个地方。

我还没有深入研究 cosine_wrap，但它似乎可以（可能在 cython 中）

xy = np.dot(X, X.T)
# or xy = np.einsum('ij,kj',X,X)

d = np.zeros((3,3),float)   # square receiver
d2 = []                     # condensed receiver
for i in range(3):
    for j in range(i+1,3):
         val=1-xy[i,j]/(norms[i]*norms[j])
         d2.append(val)
         d[j,i]=d[i,j]=val

print('array')
print(d)
print('condensed',np.array(d2))

from scipy.spatial import distance
d1=distance.pdist(X,'cosine')
print('    pdist',d1)

制作：

array
[[ 0.          0.11456226  0.1573452 ]
 [ 0.11456226  0.          0.00363075]
 [ 0.1573452   0.00363075  0.        ]]

condensed [ 0.11456226  0.1573452   0.00363075]
    pdist [ 0.11456226  0.1573452   0.00363075]

distance.squareform(d1) 产生的结果与我的 d 数组相同。

我可以通过将xy 点积除以适当的norm 外积来生成相同的方阵：

dd=1-xy/(norms[:,None]*norms)
dd[range(dd.shape[0]),range(dd.shape[1])]=0 # clean up 0s

或者在使用点积之前对X 进行标准化。这似乎是scikit 版本所做的。

Xnorm = X/norms[:,None]
1-np.einsum('ij,kj',Xnorm,Xnorm)

scikit 添加了一些 cython 代码来进行更快的稀疏计算（超出 sparse.sparse 提供的那些，但使用相同的 csr 格式）：

from scipy import sparse
Xc=sparse.csr_matrix(X)

# csr_row_norm - pyx of following
cnorm = Xc.multiply(Xc).sum(axis=1)
cnorm = np.sqrt(cnorm)
X1 = Xc.multiply(1/cnorm)  # dense matrix
dd = 1-X1*X1.T

要获得具有稀疏矩阵的快速压缩形式，我认为您需要实现 X1*X1.T 的快速压缩版本。这意味着您需要了解稀疏矩阵乘法是如何实现的——在c 代码中。 scikit cython 'fast sparse' 代码也可能提供一些想法。

numpy 有一些 tri... 函数，它们是直接的 Python 代码。它不会试图通过直接实施三计算来节省时间或空间。迭代 nd 数组的矩形布局（具有形状和步幅）比执行更复杂的三角形数组的可变长度步骤更容易。精简的形式只是将空间和计算步骤减少了一半。

============

这是c函数pdist_cosine的主要部分，它迭代i和上面的j，计算dot(x[i],y[j])/(norm[i]*norm[j])。

for (i = 0; i < m; i++) 
    for (j = i + 1; j < m; j++, dm++) 
        u = X + (n * i);
        v = X + (n * j);
        cosine = dot_product(u, v, n) / (norms[i] * norms[j]);
        if (fabs(cosine) > 1.) 
            /* Clip to correct rounding error. */
            cosine = npy_copysign(1, cosine);
        
        *dm = 1. - cosine;
    

https://github.com/scipy/scipy/blob/master/scipy/spatial/src/distance_impl.h

【讨论】：

感谢您如此全面的回复。我必须尝试理解 cython 代码！让我们看看...