在 cython 中快速访问稀疏矩阵：memoryview 与字典向量

Posted 2023-02-23

【中文标题】在 cython 中快速访问稀疏矩阵：memoryview 与字典向量

【英文标题】：fast access of sparse matrix in cython: memoryview vs vector of dictionaries

【发布时间】：2021-09-24 05:46:49

【问题描述】：

我使用 cython 来加速我在 python 中的瓶颈。任务是计算稀疏矩阵的选择性逆（低于 S），该稀疏矩阵由其以 csc 格式（数据、indptr、索引）提供的cholesky 分解给出。但任务并不是很重要，最后它是一个 3 次嵌套的 for 循环，我必须快速访问 S 的元素。

当我使用完整/巨大矩阵的内存视图时

double[:,:] Sfull

然后访问条目然后算法非常快并且符合我的期望。但很明显，这只有在矩阵 Sfull 适合内存时才有可能。

我的方法是使用字典/地图的列表/向量，这样我也可以相对快速地访问元素。

cdef vector[map[int, double]] S

事实证明，使用这种数据结构访问循环内的元素要慢大约 20 倍。这是预期的还是有其他问题？你看到其他数据结构了吗？

非常感谢您的任何 cmets 或帮助！ 最好的， 曼努埃尔

下面是 cython 代码，其中带有完整内存视图的版本被注释掉了。

cdef int invTakC12( double[:] id_diag, double[:] data, int len_i, int[:] indptr, int[:] indices, double[:, :] Sfull):


cdef vector[map[int, double]] S = testDictC(len_i-1) #list of empty dicts
cdef int i, j, j_i, lc
cdef double q 




for i in range(len_i-2, -1, -1):

    for j_i in range(indptr[i+1]-1, indptr[i]-1, -1):
        j = indices[j_i]

        q = 0

        for lc in range(indptr[i+1] -1, indptr[i], -1):

            q += data[lc] * S[j][ indices[lc] ]
            #q += data[lc] * Sfull[ indices[lc], j ]
  
        S[j][i] = -q
        #Sfull[i,j] = -q
        

        if i==j:
            S[j][i] +=  id_diag[i]
            #Sfull[i,j] += id_diag[i]
        else:
            S[i][j] -=  q
            #Sfull[j,i] -= id_diag[i]
    

return 0

【问题讨论】：

使用unordered_map 可能比使用map 更好。很多操作都更快

【参考方案1】：

您可以独立访问数组 - 例如：

cdef double[:] S_data = S.data
cdef np.int32_t[:] S_ind = S.indices
cdef np.int32_t[:] S_indptr = S.indptr

如果太不方便，可以将它们作为指针放在 C 结构中：

cdef struct mapped_csc:
    double *data
    np.int32_t *indices
    np.int32_t *indptr