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

Question

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

我需要非常有效地计算每对样本之间的余弦距离。我一直在X中使用sklearn pairwise_distances向量子集的函数，这给了我一个密集矩阵 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 的矩阵。

Answer 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的"快速稀疏"代码也可能提供想法。

numpy有一些tri...函数，这些函数是直接的Python代码。 它不会试图通过直接实现 tri 计算来节省时间或空间。 迭代 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