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

Я копался в коде обеих версий и думаю, что понимаю, что делают обе.

Начните с маленького простого 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