fast access of sparse matrix in cython: memoryview vs vector of dictionaries

Viewed 106

I used cython to speed my bottleneck in python. The task is to compute the selective inverse (below S) of a sparse matrix given by its cholesky factorization provided in csc-format (data, indptr, indices). But the the task is not really important, in the end it is a 3 times nested for-loop where I have to access elements of S fast.

When I use a memoryview of a full/huge matrix

double[:,:] Sfull

and access the entries then the algorithm is quite fast and meets my expectations. But it is clear, that this is only possible when the matrix Sfull fits into the memory.

My approach was to use a list/vector of dictionaries/maps such that I can access the elements also relatively fast.

cdef vector[map[int, double]] S

It turned out, that accessing the elements inside the loop with this data structure is around 20 times slower. Is this expected or is there another issue? Do you see any other data structure?

Thank you very much for any comments or help! Best, Manuel

Bellow, the cython code, where the version with the full memoryview is commented out.

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

You can access the arrays independently - e.g.:

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

If that's too inconvenient, you can put them in a C struct as pointers:

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