2016-05-06 24 views
3

私は大きなスパースバイナリ行列を扱っています。私はScipyスパース行列の実装を使用してそれらを凝縮しました。 scipy.spatial.distanceからJaccard distanceの計算は、スパース行列に直接操作をサポートしていないので、次のいずれかgetrow()に対するScipyスパース行列の代替方法

  1. 緻密に全体疎行列を変換した後、空腹メモリ

    あるベクターとして、各行で動作しますまたは

  2. ループを散りばめ、各行をつかみ、getrow()を使用して操作してください。

    または

  3. スパース行列上で動作するように私たちの独自の実装を記述します。

がパースペクティブに物事を置くために、ここではサンプルコードです:

import scipy.spatial.distance as d 
import numpy as np 
from scipy.sparse import csr_matrix 

# benchmark performance 
X = np.random.random((3000, 3000)) 
# binarize 
X[X > 0.3] = 0 
X[X>0] = 1 
mat = csr_matrix(X) 

a = np.zeros(3000) 
a[4] = a[100] = a[22] =1 
a = csr_matrix(a) 

def jaccard_fast(v1,v2): 
    common = v1.dot(v2.T) 
    dis = (v1 != v2).getnnz() 
    if common[0,0]: 
     return 1.0-float(common[0,0])/float(common[0,0]+dis) 
    else: 
     return 0.0 

def benchmark_jaccard_fast(): 
    for i in range(mat.shape[0]): 
     jaccard_fast(mat.getrow(i),a) 

def benchmark_jaccard_internal_todense(): 
    for v1,v2 in zip(mat.todense(),a.todense()): 
     d.jaccard(v1,v2) 

def benchmark_jaccard_internal_getrow(): 
    for i in range(mat.shape[0]): 
     d.jaccard(mat.getrow(i),a) 


print "Jaccard Fast:" 
%time benchmark_jaccard_fast() 
print "Jaccard Scipy (expanding to dense):" 
%time benchmark_jaccard_internal_todense() 
print "Jaccard Scipy (using getrow):" 
%time benchmark_jaccard_internal_getrow() 

jaccard_fastは、私自身の実装です。私の実装は、scipyの疎な行列では内部のものより速いようですが、getrow()は実装が遅くなっているようです。 scipy.spatial.distance.jaccardに対するIベンチマークjaccard_fastように、結果は以下のとおりです。

Jaccard Fast: 
CPU times: user 1.28 s, sys: 0 ns, total: 1.28 s 
Wall time: 1.28 s 
Jaccard Scipy (expanding to dense): 
CPU times: user 28 ms, sys: 8 ms, total: 36 ms 
Wall time: 37.2 ms 
Jaccard Scipy (using getrow): 
CPU times: user 1.82 s, sys: 0 ns, total: 1.82 s 
Wall time: 1.81 s 

getrowボトルネックを回避する方法上の任意の助けいただければ幸いです。私は、メモリ制限のためtodense()を使用して自分の疎な行列を拡張する余裕がありません。

+0

私は思いますあなたの 'benchmark_jaccard_internal_todense'は不正です。 3000回の完全な反復は行っていません。 'a.todense()'は(1,3000)なので、その上でzipし、 'mat'は最小の行数で繰り返します。 1. – hpaulj

+0

'benchmark_jaccard_internal_getrow'は有効なテストではありません。あなた自身の承認で' d.jaccard'がスパースで動作しないからです。 – hpaulj

答えて

2

スパースインデックスは、より遅いことが知られています。 How to read/traverse/slice Scipy sparse matrices (LIL, CSR, COO, DOK) faster?

In [33]: timeit for row in mat: x=row # sparse iteration 
1 loops, best of 3: 510 ms per loop 

In [35]: timeit for row in mat.todense(): x=row # dense iteration 
10 loops, best of 3: 175 ms per loop 

が、私はgetrow要因

In [40]: mrow=mat.getrow(0) 

In [41]: mrowd=mrow.todense() 

In [42]: timeit d.jaccard(mrow, a) # one sparse row 
1000 loops, best of 3: 1.32 ms per loop 

In [43]: timeit d.jaccard(mrow.todense(), a.todense()) # with conversion 
1000 loops, best of 3: 539 µs per loop 

In [44]: timeit d.jaccard(mrowd, ad) # dense 
10000 loops, best of 3: 173 µs per loop 

を排除スパース行列

In [36]: ad=a.todense() 

In [37]: timeit for row in mat.todense(): d.jaccard(row,ad) # all dense 
1 loops, best of 3: 734 ms per loop 

In [38]: timeit for row in mat: d.jaccard(row.todense(),ad) # inner dense 
1 loops, best of 3: 1.69 s per loop 

In [39]: timeit for row in mat: d.jaccard(row,a) # all sparse 
1 loops, best of 3: 4.61 s per loop 

を使用しているとき、あなたのd.jacardも遅くなることがわかり=========== ===========

これらのテストを再実行する必要があるのは、d.jaccardは疎では機能しません(密度が高い場合はjaccard_fastは機能しません)。したがって、疎行の繰り返しの問題をjaccard計算から分離することは、より多くの作業を必要とします。それは密な行のd.jaccard走行に一致する値を返します

def my_jaccard(mat, a): 
    common = mat*a.T # sparse does the large matrix product well 
    dis=np.array([(row!=a).getnnz() for row in mat]) # iterative 
    cA = common.A.ravel() 
    return 1 - cA/(cA + dis) 

私はjaccard_fastを少し手直ししました。 d.jaccardは、commonが0である行に対して1を返します。cAdisが同じスロットで0である可能性がない限り、cAテストは必要ないようです。

In [141]: r=np.array([d.jaccard(row,ad) for row in mat.todense()]) 

In [142]: r1=my_jaccard(mat,a) 

In [143]: np.allclose(r,r1) 
Out[143]: True 

速度は半分です。私が再加工できる場合は、dis calcと同様の速度が必要です。

In [144]: timeit r=np.array([d.jaccard(row,ad) for row in mat.todense()]) 
1 loops, best of 3: 783 ms per loop 

In [145]: timeit r1=my_jaccard(mat,a) 
1 loops, best of 3: 1.29 s per loop 

さらにcalcを微調整します。 commonの値は0にマスクされています。これは2つの目的を持っています。つまり、0で除算する問題が発生しないようにし、disの反復回数を減らしてスピードを少し改善します。

def my_jaccard(mat, a): 
    common=mat*a.T 
    cA = common.A.ravel() 
    mask = cA!=0 
    cA = cA[mask] 
    dis = np.array([(row!=a).getnnz() for row, b in zip(mat,mask) if b]) 
    ret = np.ones(mat.shape[0]) 
    ret[mask] = 1 - cA/(cA+dis) 
    return ret 

これは少し時間がかかります。

In [188]: timeit my_jaccard(mat,a) 
1 loops, best of 3: 1.04 s per loop 

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

私は疎行列を比較を見て、その質問にはPython - Efficient Function with scipy sparse Matrices

の問題で重複があります1行の行列で、行を複製するためにsparse.kronを使用することが、numpyブロードキャストを複製する最も速い方法であることが判明しました。有意に(10倍)を改善するこのタイミングでdisアレイ

def my_jaccard1(mat, a): 
    common = mat*a.T 
    cA = common.A.ravel() 
    aM = sparse.kron(a,np.ones((mat.shape[0],1),int)) 
    dis = (mat!=aM).sum(1) 
    ret = 1-cA/(cA+dis.A1) 
    return ret  

を計算するjaccardにそのアイデアを使用

In [318]: timeit my_jaccard1(mat,a) 
1 loops, best of 3: 97.1 ms per loop 

ゼロによる除算から保護するために、前のように、私はマスキング適用することができます。実際には計算が遅くなります(140msまで)。

def my_jaccard3(mat, a): 
    common = mat*a.T 
    cA = common.A.ravel() 
    mask = cA!=0 
    cA = cA[mask] 
    aM = sparse.kron(a,np.ones((len(cA),1),int)) 
    dis = (mat[mask,:]!=aM).sum(1) 
    ret = np.ones(mat.shape[0]) 
    ret[mask] = 1 - cA/(cA+dis.A1) 
    return ret 

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

編集 - 疑われる例テスト

In [75]: x,y= np.array([1,1,0,0,1,0]), np.array([0,0,1,0,1,0]) 

In [76]: d.jaccard(x,y) 
Out[76]: 0.75 

In [78]: jaccard_fast(sparse.csr_matrix(x),sparse.csr_matrix(y)) 
Out[78]: 0.75 

マイバージョン:

In [79]: my_jaccard(sparse.csr_matrix(x),sparse.csr_matrix(y)) 
Out[79]: array([ 0.75]) 
... 
In [82]: my_jaccard3(sparse.csr_matrix(x),sparse.csr_matrix(y)) 
Out[82]: array([ 0.75]) 

(編集 - 明示的sparse.kronを使用)

+0

これを前にチェックできないという謝罪がありますが、あなたの解決策は 'd.jaccard'と同じ結果をもたらさないと思います。たとえば、ベクトル[1,1,0,0,1,0]、[0,0,1,0,1,0]で試してください。 – Segmented

+0

'jacard_fast'はそのケースを正しく取得しますか?私はそれを複製しようとしました(1 0スイッチを除く)。 – hpaulj

+0

'jacard_fast'はこのケースを正しく取得します。確認しました。お詫び申し上げます、私は明日もまた計算装置にアクセスする際にこれを調べなければなりません。 – Segmented

関連する問題