python - 三角/スパースストレージへのnumpy行列乗算?

Question

私は非常に大きな疎行列乗算 (matmul) の問題に取り組んでいます。例として、次のように言いましょう。

A はバイナリ ( 75 x 200,000 ) 行列です。まばらなので、保存には csc を使用しています。次の matmul 操作を行う必要があります。
B = A.transpose() * A
出力は、サイズが 200Kx200K の疎な対称行列になります。

残念ながら、B はラップトップの RAM (または「コア内」) に格納するには大きすぎます。一方、この問題を解決する B のプロパティがいくつかあるので、私は幸運です。

B は対角線とスパースに沿って対称になるため、三角行列 (上/下) を使用して matmul 演算の結果を格納し、スパース行列の格納形式を使用してサイズをさらに縮小できます。

私の質問は... numpyまたはscipyに出力ストレージ要件がどのようになるかを事前に伝えて、numpyを使用してストレージソリューションを選択し、数回後に「マトリックスが大きすぎます」というランタイムエラーを回避できるようにすることはできますか?分（時間）の計算？

言い換えると、行列乗算のストレージ要件は、近似カウントアルゴリズムを使用して 2 つの入力行列の内容を分析することで概算できますか?

https://en.wikipedia.org/wiki/Approximate_counting_algorithm

そうでない場合は、ブルートフォースソリューションを検討しています。次の Web リンクからの map/reduce、アウトオブコアストレージ、または matmul サブディビジョンソリューション (strassens アルゴリズム) に関連するもの:

いくつかの Map/Reduce 問題の細分化ソリューション

アウトオブコア (PyTables) ストレージソリューション

Python と NumPy を使用した非常に大きな行列

matmul サブディビジョンソリューション:

推奨事項、コメント、またはガイダンスをお寄せいただきありがとうございます。

score 3 · Accepted Answer

行列とその転置の積を求めているため、値は基本的に列と元の行列[m, n]のドット積になります。mn

次のマトリックスをおもちゃの例として使用します

a = np.array([[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
              [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
              [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1]])
>>> np.dot(a.T, a)
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
       [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
       [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
       [0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 2]])

これは一定の形(3, 12)をしており、ゼロ以外のエントリが 7 つあります。それとの転置の積はもちろん形状(12, 12)であり、16 個のゼロ以外のエントリがあり、そのうちの 6 個は対角線にあるため、必要な要素は 11 個だけです。

次の 2 つの方法のいずれかを使用して、出力行列のサイズを把握できます。

CSRフォーマット

元のマトリックスにCゼロ以外の列がある場合、新しいマトリックスには最大でもC**2ゼロ以外のエントリがあり、そのエントリCは対角にあり、ゼロにならないことが保証されます。残りのエントリのうち、保持する必要があるのは半分だけです。それはせいぜい(C**2 + C) / 2ゼロ以外の要素です。もちろん、これらの多くもゼロになるため、これはおそらく過大評価です。

マトリックスがcsr形式で格納されている場合indices、対応するオブジェクトの属性には、scipyゼロ以外のすべての要素の列インデックスを含む配列があるため、上記の推定値を次のように簡単に計算できます。

>>> a_csr = scipy.sparse.csr_matrix(a)
>>> a_csr.indices
array([ 2, 11,  1,  7, 10,  4, 11])
>>> np.unique(a_csr.indices).shape[0]
6

したがって、ゼロ以外のエントリを持つ列が 6 つあるため、実際の 16 よりもはるかに多く、最大で 36 の非ゼロエントリが推定されます。

CSCフォーマット

非ゼロ要素の列インデックスの代わりに行インデックスがある場合、実際にはより良い推定を行うことができます。2 つの列の内積が非ゼロになるには、同じ行にゼロ以外の要素が含まれている必要があります。R特定の行にゼロ以外の要素がある場合、それらはR**2製品にゼロ以外の要素を提供します。これをすべての行で合計すると、一部の要素を複数回カウントすることになるため、これは上限でもあります。

行列の非ゼロ要素の行インデックスはindices疎な csc 行列の属性にあるため、この推定は次のように計算できます。

>>> a_csc = scipy.sparse.csc_matrix(a)
>>> a_csc.indices
array([1, 0, 2, 1, 1, 0, 2])
>>> rows, where = np.unique(a_csc.indices, return_inverse=True)
>>> where = np.bincount(where)
>>> rows
array([0, 1, 2])
>>> where
array([2, 3, 2])
>>> np.sum(where**2)
17

これは実際の 16 に非常に近いです。そして、この推定値が実際に次の値と同じであることは、実際には偶然ではありません。

>>> np.sum(np.dot(a.T,a),axis=None)
17

いずれにせよ、次のコードを使用すると、推定がかなり適切であることがわかります。

def estimate(a) :
    a_csc = scipy.sparse.csc_matrix(a)
    _, where = np.unique(a_csc.indices, return_inverse=True)
    where = np.bincount(where)
    return np.sum(where**2)

def test(shape=(10,1000), count=100) :
    a = np.zeros(np.prod(shape), dtype=int)
    a[np.random.randint(np.prod(shape), size=count)] = 1
    print 'a non-zero = {0}'.format(np.sum(a))
    a = a.reshape(shape)
    print 'a.T * a non-zero = {0}'.format(np.flatnonzero(np.dot(a.T,
                                                                a)).shape[0])
    print 'csc estimate = {0}'.format(estimate(a))

>>> test(count=100)
a non-zero = 100
a.T * a non-zero = 1065
csc estimate = 1072
>>> test(count=200)
a non-zero = 199
a.T * a non-zero = 4056
csc estimate = 4079
>>> test(count=50)
a non-zero = 50
a.T * a non-zero = 293
csc estimate = 294

python - 三角/スパースストレージへのnumpy行列乗算?

1 に答える 1

Related

Reference