Как получить сжатую форму парных расстояний напрямую?

Question

У меня очень большая скудная разреженная матрица 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.

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. Он не пытается сэкономить время или пространство путем непосредственной реализации трех вычислений. Проще выполнить итерацию по прямоугольному расположению массива 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