首页 文章

如何使用置换数组有效地置换稀疏(Numpy)矩阵中的行?

提问于
浏览
0

我使用Scipy Reverse Cuthill-McKee实现(scipy.sparse.csgraph.reverse_cuthill_mckee)来使用(高维)稀疏csr_matrix创建带矩阵 . 这个方法的结果是一个置换数组,它给出了如我所知的如何置换矩阵行的索引 .

现在有没有有效的解决方案在我的稀疏csr_matrix中对任何其他稀疏矩阵(csr,lil_matrix等)进行这种排列?我尝试了一个for循环,但我的矩阵的尺寸为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调用的结果是一个数组,就像包含我的排列索引的tupel一样 . 所以这个数组是这样的:[199999 54877 54873 ...,12045 9191 0](size = 200,000)

这意味着:索引为0的行现在索引为199999,索引为1的行现在索引为54877,索引为2的行现在索引为54873等,请参阅:https://en.wikipedia.org/wiki/Permutation#Definition_and_notations(据我所知返回)

谢谢

1 回答

  • 1

    我想知道你是否正确应用了排列数组 .

    创建一个随机矩阵(float)并将其转换为 uint8 (注意, csr 计算可能不适用于此dtype):

    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 中,我没有看到1s聚类的任何改进 .


    相当于MATLAB的文档说

    r = symrcm(S)返回S的对称反Cuthill-McKee排序 . 这是一个置换r,使得S(r,r)倾向于使其非零元素更接近对角线 .

    numpy 条款中,这意味着:

    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个波段 .

    并使用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)
    

相关问题