三角/疎ストレージへのnumpy行列の乗算？

Question

私は非常に大きな疎な行列乗算（matmul）問題を扱っています。一例としてのように言ってみましょう。

• Aは、バイナリ（75×20万）行列です。それは疎なので、私はストレージのためにcscを使用しています。

• B = A.transpose（）*

• 出力サイズ200Kx200Kのスパース対称行列になるだろう：私は次MATMUL操作を行う必要があります。

残念ながら、Bは私のラップトップ上（または「コア内」）RAMに格納する大へ道になるだろう。一方、私はこの問題を解決すべきいくつかのプロパティがBにあるので私はラッキーです。

Bは対角線と疎で対称になるので、三角行列（上/下）を使用してmatmul演算の結果を格納することができ、疎行列格納形式でサイズをさらに小さくすることができます。

私の質問は... numpyを使用してストレージソリューションを選択できるように出力ストレージの要件がどのように表示されるのかを前もって知ることができます。計算の数分（数時間）後のランタイムエラー？

換言すれば、行列乗算の記憶要件は、近似計数アルゴリズムを用いて2つの入力行列の内容を分析することによって近似できるか？

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

ない場合、私は力ずくのソリューションに探しています。マップ/削減、アウトオブコアストレージ、または以下のWebリンクからMATMUL細分化ソリューション（strassensアルゴリズム）が関与何か：

カップル地図/削減問題の細分化ソリューションを

• http://www.norstad.org/matrix-multiply/index.html
• http://bpgergo.blogspot.com/2011/08/matrix-multiplication-in-python.html

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

• Very large matrices using Python and NumPy

MATMUL細分化ソリューション：

• https://en.wikipedia.org/wiki/Strassen_algorithm
• http://facultyfp.salisbury.edu/taanastasio/COSC490/Fall03/Lectures/FoxMM/example.pdf
• http://eli.thegreenplace.net/2012/01/16/python-parallelizing-cpu-bound-tasks-with-multiprocessing/

任意の推奨事項については、事前のおかげで、コメント、または指導！

Answer 1

転置行列の積の後であるため、[m, n]の値は、基本的に元の行列の列mとnの内積になります。私は

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つの方法のいずれかであることを行っているものの良いアイデアを得ることができます：あなたの元の行列がC非ゼロの列がある場合

CSRのFORMAT

をあなたの新しい行列は最大でC**2の非ゼロエントリーを持ち、そのうちCは対角にあり、ゼロでないことが保証され、残りのエントリーのうち半分だけを保持する必要があるので、最大でも(C**2 + C)/2です。ゼロ要素。もちろん、これらの多くはゼロになるので、おそらくこれは過大な過大評価です。

あなたのマトリックスはcsr形式で格納されている場合、対応するscipyオブジェクトのindices属性は、すべての非ゼロ要素の列インデックスを持つ配列を有しているので、簡単に上記のように推定値を計算することができる：

>>> 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

CSCのFORMATより方法より、せいぜい36非ゼロのエントリ、ためになる

ゼロ以外の要素の列インデックスの代わりに行インデックスがある場合、実際にはより良い見積もりを行うことができます。 2つの列のドット積が非ゼロであるためには、それらは同じ行に非ゼロ要素を持たなければならない。指定された行に0以外の要素がある場合は、R**2以外の要素を製品に提供します。すべての行についてこれを合計すると、いくつかの要素を複数回カウントするようにバインドされているため、これも上限です。あなたの行列の非零要素の

行インデックスは、スパースCSCマトリックスのindices属性であるので、次のように、この推定値を計算することができる。

>>> 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