パーミュテーション配列を使ってスパース（Numpy）行列の行を効率的に並べ替えるにはどうすればいいですか？

Question

（高次元の）スパースcsr_matrixを使用してバンド行列を作成するために、Scipy Reverse Cuthill-McKee実装（scipy.sparse.csgraph.reverse_cuthill_mckee）を使用しました。 このメソッドの結果は、私が理解したように私の行列の行をどのように並べ替えるかの指標を私に与える並べ替え配列です。

他のスパース行列（csr、lil_matrixなど）のスパースcsr_matrixでこの順列を実行する効率的なソリューションはありますか？ 私はfor-loopを試しましたが、私の行列は200,000 x 150,000のような次元を持ち、時間がかかりすぎます。

A = csr_matrix((data,(rowind,columnind)), shape=(200000, 150000), dtype=np.uint8) 

permutation_array = csgraph.reverse_cuthill_mckee(A, false) 

result_matrix = lil_matrix((200000, 150000), dtype=np.uint8) 

i=0 
for x in np.nditer(permutation_array): 
    result_matrix[x, :]=A[i, :] 
    i+=1 

reverse_cuthill_mckee呼び出しの結果は、私の順列のためのインデックスを含むタプルに似ている配列です。これは、[199999 54877 54873 ...、12045 9191 0]（サイズ=20万）

：インデックス0を有する 行は、現在のインデックス199999を有する、インデックス1と 行は、現在のインデックス54877を有しているので、この配列は次のようです、インデックス2との 行は現在のインデックス54873、 などを参照していますhttps://en.wikipedia.org/wiki/Permutation#Definition_and_notations （私はリターンを理解できるように）

は、あなたが正しく順列配列を適用している場合、私は疑問に思う

Answer 1

ありがとうございます。

（csr計算はこのDTYPEでは動作しない場合があります、注意してください）ランダム行列（フロート）を作成し、uint8に変換します

In [963]: ran=sparse.random(10,10,.3, format='csr') 
In [964]: A = sparse.csr_matrix((np.ones(ran.data.shape).astype(np.uint8),ran.indices, ran.indptr)) 
In [965]: A.A 
Out[965]: 
array([[1, 1, 0, 0, 0, 0, 1, 0, 0, 0], 
     [0, 1, 1, 1, 1, 1, 1, 0, 1, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [1, 1, 0, 0, 0, 0, 0, 1, 0, 1], 
     [0, 1, 0, 0, 1, 1, 0, 0, 0, 0], 
     [1, 0, 1, 0, 0, 1, 0, 1, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 1, 0, 0, 0, 1], 
     [0, 1, 1, 1, 0, 1, 0, 0, 0, 0], 
     [0, 0, 0, 0, 1, 1, 1, 0, 0, 0]], dtype=uint8) 

（おっと、ここで間違ったマトリックスを使用）：

を

In [994]: permutation_array = csgraph.reverse_cuthill_mckee(A, False) 
In [995]: permutation_array 
Out[995]: array([9, 7, 0, 4, 6, 3, 5, 1, 8, 2], dtype=int32) 

私の最初の傾きが元の行列の単純インデックス行にこのような配列を使用することである。

In [996]: A[permutation_array,:].A 
Out[996]: 
array([[0, 0, 0, 0, 1, 1, 1, 0, 0, 0], 
     [0, 0, 0, 0, 0, 1, 0, 0, 0, 1], 
     [1, 1, 0, 0, 0, 0, 1, 0, 0, 0], 
     [0, 1, 0, 0, 1, 1, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [1, 1, 0, 0, 0, 0, 0, 1, 0, 1], 
     [1, 0, 1, 0, 0, 1, 0, 1, 0, 0], 
     [0, 1, 1, 1, 1, 1, 1, 0, 1, 0], 
     [0, 1, 1, 1, 0, 1, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) 

クラスタリングがあります。多分ランダムな行列から期待できる最高のものです。一方、

あなたはやっているように見える：

In [997]: res = sparse.lil_matrix(A.shape,dtype=A.dtype) 
In [998]: res[permutation_array,:] = A 
In [999]: res.A 
Out[999]: 
array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 0, 0, 0, 0, 1, 0, 0, 0, 1], 
     [0, 0, 0, 0, 1, 1, 1, 0, 0, 0], 
     [1, 0, 1, 0, 0, 1, 0, 1, 0, 0], 
     [1, 1, 0, 0, 0, 0, 0, 1, 0, 1], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [0, 1, 0, 0, 1, 1, 0, 0, 0, 0], 
     [0, 1, 1, 1, 1, 1, 1, 0, 1, 0], 
     [0, 1, 1, 1, 0, 1, 0, 0, 0, 0], 
     [1, 1, 0, 0, 0, 0, 1, 0, 0, 0]], dtype=uint8) 

私はres内の1のクラスタリングにおける任意の改善が表示されません。

---

MATLAB相当用ドキュメントは

R = symrcm（S）を言うSの対称逆Cuthill-McKeeでの並べ替えを返すこれは順列RようにS（R、R）であります対角に近い非ゼロ要素を持つ傾向があります。

In [1019]: I,J=np.ix_(permutation_array,permutation_array) 
In [1020]: A[I,J].A 
Out[1020]: 
array([[0, 0, 0, 1, 1, 0, 1, 0, 0, 0], 
     [1, 0, 0, 0, 0, 0, 1, 0, 0, 0], 
     [0, 0, 1, 0, 1, 0, 0, 1, 0, 0], 
     [0, 0, 0, 1, 0, 0, 1, 1, 0, 0], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0], 
     [1, 1, 1, 0, 0, 0, 0, 1, 0, 0], 
     [0, 1, 1, 0, 0, 0, 1, 0, 0, 1], 
     [0, 0, 0, 1, 1, 1, 1, 1, 1, 1], 
     [0, 0, 0, 0, 0, 1, 1, 1, 0, 1], 
     [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) 

および確か2つのオフ対角のその他0バンドがあります

はnumpy言えば、それが意味します。

そして、MATLABページ上の帯域幅の計算を使用して、https://www.mathworks.com/help/matlab/ref/symrcm.html

In [1028]: i,j=A.nonzero() 
In [1029]: np.max(i-j) 
Out[1029]: 7 
In [1030]: i,j=A[I,J].nonzero() 
In [1031]: np.max(i-j) 
Out[1031]: 5 

MATLABドキュメントはこの順列で、固有値が同じままであることを言います。テスト：

In [1032]: from scipy.sparse import linalg 
In [1048]: linalg.eigs(A.astype('f'))[0] 
Out[1048]: 
array([ 3.14518213+0.j  , -0.96188843+0.j  , 
     -0.58978939+0.62853903j, -0.58978939-0.62853903j, 
     1.09950364+0.54544497j, 1.09950364-0.54544497j], dtype=complex64) 
In [1049]: linalg.eigs(A[I,J].astype('f'))[0] 
Out[1049]: 
array([ 3.14518023+0.j  , 1.09950352+0.54544479j, 
     1.09950352-0.54544479j, -0.58978981+0.62853914j, 
     -0.58978981-0.62853914j, -0.96188819+0.j  ], dtype=complex64) 

固有値は、私たちが先にしようとした行の順列のために同じではありません。

In [1050]: linalg.eigs(A[permutation_array,:].astype('f'))[0] 
Out[1050]: 
array([ 2.95226836+0.j  , -1.60117996+0.52467293j, 
     -1.60117996-0.52467293j, -0.01723826+1.06249797j, 
     -0.01723826-1.06249797j, 0.90314150+0.j  ], dtype=complex64) 
In [1051]: linalg.eigs(res.astype('f'))[0] 
Out[1051]: 
array([-0.05822830-0.97881651j, -0.99999994+0.j  , 
     1.17350495+0.j  , -0.91237622+0.8656373j , 
     -0.91237622-0.8656373j , 2.26292515+0.j  ], dtype=complex64) 

この[I,J]順列は、[スワップのhttp://ciprian-zavoianu.blogspot.com/2009/01/project-bandwidth-reduction.html

In [1058]: B = np.matrix('1 0 0 0 1 0 0 0;0 1 1 0 0 1 0 1;0 1 1 0 1 0 0 0;0 0 0 
     ...: 1 0 0 1 0;1 0 1 0 1 0 0 0; 0 1 0 0 0 1 0 1;0 0 0 1 0 0 1 0;0 1 0 0 0 
     ...: 1 0 1') 
In [1059]: B 
Out[1059]: 
matrix([[1, 0, 0, 0, 1, 0, 0, 0], 
     [0, 1, 1, 0, 0, 1, 0, 1], 
     [0, 1, 1, 0, 1, 0, 0, 0], 
     [0, 0, 0, 1, 0, 0, 1, 0], 
     [1, 0, 1, 0, 1, 0, 0, 0], 
     [0, 1, 0, 0, 0, 1, 0, 1], 
     [0, 0, 0, 1, 0, 0, 1, 0], 
     [0, 1, 0, 0, 0, 1, 0, 1]]) 
In [1060]: Bm=sparse.csr_matrix(B) 
In [1061]: Bm 
Out[1061]: 
<8x8 sparse matrix of type '<class 'numpy.int32'>' 
    with 22 stored elements in Compressed Sparse Row format> 
In [1062]: permB = csgraph.reverse_cuthill_mckee(Bm, False) 
In [1063]: permB 
Out[1063]: array([6, 3, 7, 5, 1, 2, 4, 0], dtype=int32) 
In [1064]: Bm[np.ix_(permB,permB)].A 
Out[1064]: 
array([[1, 1, 0, 0, 0, 0, 0, 0], 
     [1, 1, 0, 0, 0, 0, 0, 0], 
     [0, 0, 1, 1, 1, 0, 0, 0], 
     [0, 0, 1, 1, 1, 0, 0, 0], 
     [0, 0, 1, 1, 1, 1, 0, 0], 
     [0, 0, 0, 0, 1, 1, 1, 0], 
     [0, 0, 0, 0, 0, 1, 1, 1], 
     [0, 0, 0, 0, 0, 0, 1, 1]], dtype=int32)