凝縮されたペアワイズ距離を直接得る方法は？

Question

私は非常に大きなサイピーススパースcsrマトリックスを持っています。これは100,000×2,000,000次元の行列です。それをXとしましょう。各行は、2,000,000次元の空間内のサンプルベクトルです。

サンプルの各ペア間のコサイン距離を非常に効率的に計算する必要があります。私はXのベクトルのサブセットでsklearn pairwise_distances関数を使用しています。これは、稠密な行列Dを与えます。冗長なエントリを含む対の距離の2乗形式です。 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コサインない：

_row_normsが行ドットである

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

からeinsumを使用して：だから

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

これは、最初の場所ですXは配列でなければなりません。生産

（おそらくcythonに）私はcosine_wrapに掘っていないが、何に見える

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配列として同じものを生産します。

Iは、適切なnorm外積とxyドット積を割ることによって同じ正方形アレイを生成することができる：

​​

またはドット積をとる前に、Xを正規化。これはバージョンscikitのように見えます。

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

scikit（sparse.sparseによって提供されるものを超えて、同じcsr形式を使用して）より速く、疎な計算を行うために、いくつかのcythonコードを追加しました：

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には、まっすぐ進むPythonコードであるtri...の機能があります。トライ計算を直接実装することによって時間や空間を節約しようとするものではありません。三角形配列のより複雑な可変長のステップを行うよりも、nd配列の長方形のレイアウト（形状とストライドを持つ）を反復するほうが簡単です。凝縮された形態は、空間と計算ステップを半分だけカットする。

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

ここdot(x[i],y[j])/(norm[i]*norm[j])を計算iと上jを反復処理、c機能pdist_cosine、の主要部分です。あなたが具体的な例を追加する必要があります

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