Непосредственно используйте библиотеку Intel mkl на Scipy разреженной матрице для вычисления точки A.T с меньшей памятью

Question

Я хочу позвонить mkl.mkl_scsrmultcsr из python. Цель состоит в том, чтобы вычислить разреженную матрицу C в формате compressed sparse row. Разреженная матрица C является матричным продуктом между A и транспонированием A, где A также является разреженной матрицей в формате csr. При вычислении C = A dot (AT) с scipy, scipy кажется (?), Чтобы выделить новую память для переноса транспонирования A (AT) и определенно выделяет память для новой матрицы C (это означает, что я не могу использовать существующий C матрица). Поэтому я хочу попытаться использовать функцию mkl c напрямую, чтобы уменьшить использование памяти.

Here - это ответ, который работает для другой функции mkl. В этом ответе функция mkl была быстрее в 4 раза.

Следующая версия, наконец, работает после работы над ней в течение 4 дней. Я хочу знать, потерял ли я память. Будет ли ctype делать любую копию массива numpy? Является ли это необходимым преобразование csr-> csc в тестовую функцию? Функция intel c может вычислять (A.T) точку A или точку A, но не точку (A.T).

Еще раз спасибо.

from ctypes import * 
import scipy.sparse as spsp 
import numpy as np 
import multiprocessing as mp 
# June 2nd 2016 version. 

# Load the share library 
mkl = cdll.LoadLibrary("libmkl_rt.so") 

def get_csr_handle(A,clear=False): 
    if clear == True: 
     A.indptr[:] = 0 
     A.indices[:] = 0 
     A.data[:] = 0 

    a_pointer = A.data.ctypes.data_as(POINTER(c_float)) 
    # Array containing non-zero elements of the matrix A. 
    # This corresponds to data array of csr_matrix 
    # Its length is equal to #non zero elements in A 
    # (Can this be longer than actual #non-zero elements?) 
    assert A.data.ctypes.data % 16 == 0 # Check alignment 

    ja_pointer = A.indices.ctypes.data_as(POINTER(c_int)) 
    # Array of column indices of all non-zero elements of A. 
    # This corresponds to the indices array of csr_matrix 
    assert A.indices.ctypes.data % 16 == 0 # Check alignment 

    ia_pointer = A.indptr.ctypes.data_as(POINTER(c_int)) 
    # Array of length m+1. 
    # a[ia[i]:ia[i+1]] is the value of nonzero entries of 
    # the ith row of A. 
    # ja[ia[i]:ia[i+1]] is the column indices of nonzero 
    # entries of the ith row of A 
    # This corresponds to the indptr array of csr_matrix 
    assert A.indptr.ctypes.data % 16 == 0 # Check alignment 

    A_data_size = A.data.size 
    A_indices_size = A.indices.size 
    A_indptr_size = A.indptr.size 

    return (a_pointer, ja_pointer, ia_pointer, A) 

def csr_dot_csr_t(A_handle, C_handle, nz=None): 
    # Calculate (A.T).dot(A) and put result into C 
    # 
    # This uses one-based indexing 
    # 
    # Both C.data and A.data must be in np.float32 type. 
    # 
    # Number of nonzero elements in C must be greater than 
    #  or equal to the size of C.data 
    # 
    # size of C.indptr must be greater than or equal to 
    #  1 + (num rows of A). 
    # 
    # C_data = np.zeros((nz), dtype=np.single) 
    # C_indices = np.zeros((nz), dtype=np.int32) 
    # C_indptr = np.zeros((m+1),dtype=np.int32) 

    #assert len(c_pointer._obj) >= 1 + A_shape[0] 


    (a_pointer, ja_pointer, ia_pointer, A) = A_handle 
    (c_pointer, jc_pointer, ic_pointer, C) = C_handle 
    #print "CCC",C 
    #assert type(C.data[0]) == np.float32 
    #assert type(A.data[0]) == np.float32 
    #assert C.indptr.size >= A.shape[0] + 1 

    #CC = A.dot(A.T) 
    #assert C.data.size >= nz 
    #assert C.indices.size >= nz 

    trans_pointer = byref(c_char('T')) 
    sort_pointer = byref(c_int(0)) 

    (m, n)   = A.shape 
    sort_pointer  = byref(c_int(0)) 
    m_pointer   = byref(c_int(m))  # Number of rows of matrix A 
    n_pointer   = byref(c_int(n))  # Number of columns of matrix A 
    k_pointer   = byref(c_int(n))  # Number of columns of matrix B 
               # should be n when trans='T' 
          # Otherwise, I guess should be m 
    ### 
    b_pointer = a_pointer 
    jb_pointer = ja_pointer 
    ib_pointer = ia_pointer 
    ### 
    if nz == None: 
     nz = n*n #*n # m*m # Number of nonzero elements expected 
      # probably can use lower value for sparse 
      # matrices. 
    nzmax_pointer = byref(c_int(nz)) 
    # length of arrays c and jc. (which are data and 
    # indices of csr_matrix). So this is the number of 
    # nonzero elements of matrix C 
    # 
    # This parameter is used only if request=0. 
    # The routine stops calculation if the number of 
    # elements in the result matrix C exceeds the 
    # specified value of nzmax. 

    info = c_int(-3) 
    info_pointer = byref(info) 
    request_pointer_list = [byref(c_int(0)), byref(c_int(1)), byref(c_int(2))] 
    return_list = [] 
    for ii in [0]: 
     request_pointer = request_pointer_list[ii] 
     ret = mkl.mkl_scsrmultcsr(trans_pointer, request_pointer, sort_pointer, 
        m_pointer, n_pointer, k_pointer, 
        a_pointer, ja_pointer, ia_pointer, 
        b_pointer, jb_pointer, ib_pointer, 
        c_pointer, jc_pointer, ic_pointer, 
        nzmax_pointer, info_pointer) 
     info_val = info.value 
     return_list += [ (ret,info_val) ] 
    return return_list 


def show_csr_internal(A, indent=4): 
    # Print data, indptr, and indices 
    # of a scipy csr_matrix A 
    name = ['data', 'indptr', 'indices'] 
    mat = [A.data, A.indptr, A.indices] 
    for i in range(3): 
     str_print = ' '*indent+name[i]+':\n%s'%mat[i] 
     str_print = str_print.replace('\n', '\n'+' '*indent*2) 
     print(str_print) 

def fix_for_scipy(C,A): 
    n = A.shape[1] 
    print "fix n", n 
    nz = C.indptr[n] - 1 # -1 as this is still one based indexing. 
    print "fix nz", nz 
    data = C.data[:nz] 

    C.indptr[:n+1] -= 1 
    indptr = C.indptr[:n+1] 
    C.indices[:nz] -= 1 
    indices = C.indices[:nz] 
    return spsp.csr_matrix((data, indices, indptr), shape=(n,n)) 

def test(): 
    AA= [[1,0,0,1], 
     [1,0,1,0], 
     [0,0,1,0]] 
    AA = np.random.choice([0,1], size=(3,750000), replace=True, p=[0.99,0.01]) 
    A_original = spsp.csr_matrix(AA) 
    #A = spsp.csr_matrix(A_original, dtype=np.float32) 
    A = A_original.astype(np.float32).tocsc() 
    #A_original = A.todense() 
    A = spsp.csr_matrix((A.data, A.indices, A.indptr)) 
    print "A:" 
    show_csr_internal(A) 
    print A.todense() 
    A.indptr += 1 # convert to 1-based indexing 
    A.indices += 1 # convert to 1-based indexing 
    A_ptrs = get_csr_handle(A) 

    C = spsp.csr_matrix(np.ones((3,3)), dtype=np.float32) 
    #C.data = C.data[:16].view() 
    #C.indptr = C.indptr 
    C_ptrs = get_csr_handle(C, clear=True) 

    print "C:" 
    show_csr_internal(C)  
    print "=call mkl function=" 
    return_list= csr_dot_csr_t(A_ptrs, C_ptrs) 
    print "(ret, info):", return_list 
    print "C after calling mkl:" 
    show_csr_internal(C)  

    C_fix = fix_for_scipy(C,A) 
    print "C_fix for scipy:" 
    show_csr_internal(C_fix) 
    print C_fix.todense() 

    print "Original C after fixing:" 
    show_csr_internal(C) 

    print "scipy's (A).dot(A.T)" 
    scipy_ans = (A_original).dot(A_original.T) 
    #scipy_ans = spsp.csr_matrix(scipy_ans) 
    show_csr_internal(scipy_ans) 
    print scipy_ans.todense() 

if __name__ == "__main__": 
    test() 

Результат:

A: 
    data: 
     [ 1. 1. 1. ..., 1. 1. 1.] 
    indptr: 
     [ 0  0  0 ..., 22673 22673 22673] 
    indices: 
     [1 0 2 ..., 2 1 2] 
[[ 0. 0. 0.] 
[ 0. 0. 0.] 
[ 0. 0. 0.] 
..., 
[ 0. 0. 0.] 
[ 0. 0. 0.] 
[ 0. 0. 0.]] 
C: 
    data: 
     [ 0. 0. 0. 0. 0. 0. 0. 0. 0.] 
    indptr: 
     [0 0 0 0] 
    indices: 
     [0 0 0 0 0 0 0 0 0] 
=call mkl function= 
(ret, info): [(2, 0)] 
C after calling mkl: 
    data: 
     [ 7576. 77. 83. 77. 7607. 104. 83. 104. 7490.] 
    indptr: 
     [ 1 4 7 10] 
    indices: 
     [1 2 3 1 2 3 1 2 3] 
fix n 3 
fix nz 9 
C_fix for scipy: 
    data: 
     [ 7576. 77. 83. 77. 7607. 104. 83. 104. 7490.] 
    indptr: 
     [0 3 6 9] 
    indices: 
     [0 1 2 0 1 2 0 1 2] 
[[ 7576. 77. 83.] 
[ 77. 7607. 104.] 
[ 83. 104. 7490.]] 
Original C after fixing: 
    data: 
     [ 7576. 77. 83. 77. 7607. 104. 83. 104. 7490.] 
    indptr: 
     [0 3 6 9] 
    indices: 
     [0 1 2 0 1 2 0 1 2] 
scipy's (A.T).dot(A) 
    data: 
     [ 83 77 7576 104 77 7607 83 104 7490] 
    indptr: 
     [0 3 6 9] 
    indices: 
     [2 1 0 2 0 1 0 1 2] 
[[7576 77 83] 
[ 77 7607 104] 
[ 83 104 7490]] 

вещи узнал:

1. индекс начинается с 1 для матриц A, B и C.
2. Я перепутал c_indptr и c_indices в оригинале код. Должно быть ia = indptr scipy csr_matrix. ja = индексы scipy csr_matrix.

3. с кодом here. Все передается в mkl_? Csrmultcsr как указатель. mkl_scsrmultcsr(&ta, &r[1], &sort, &m, &m, &m, a, ja, ia, a, ja, ia, c, jc, ic, &nzmax, &info);

4. Я хочу функцию mkl, которая может работать с индексом на основе нуля. Функция mkl.mkl_scsrmultcsr может работать только с однонаправленной индексацией. (Или я один на основе индексации для всего. Это означает, что с помощью Intel C функцию вместо SciPy/NumPy для большинства линейных шагов алгебры.)

Answer 1

Посмотрите на коде Python для SciPy разреженного продукта. Обратите внимание, что он вызывает скомпилированный код за 2 прохода.

Похоже, что код MKL делает то же самое

https://software.intel.com/en-us/node/468640

Если запрос = 1, подпрограмма вычисляет только значения массива гс длины т + 1, память для этого массива должны быть выделены заранее. При выходе значение ic (m + 1) - 1 является действительным числом элементов в массивах c и jc.

Если запрос = 2, процедура была вызвана ранее с запросом параметра = 1, массивы вывода jc и c распределены в вызывающей программе и они имеют длину ic (m + 1) - 1, по крайней мере, ,

Вы первый выделил ic на основе количества строк C (вы знаете, что от входов), и вызвать код МКЛ с request=1.

Для request=2 вам необходимо выделить c и jc массивы, исходя из их размера в ic(m+1) - 1. Это не то же самое, что число nnz во входных массивах.

Вы используете request1 = c_int(0), который требует, чтобы c массивы быть правильный размер - что вы не знаете, без фактического выполнения dot (или request 1).

==================

File:  /usr/lib/python3/dist-packages/scipy/sparse/compressed.py 
Definition: A._mul_sparse_matrix(self, other) 

проход 1 выделяет indptr (размер примечание), и передает указатели (данные не имеет значения в этом передать)

indptr = np.empty(major_axis + 1, dtype=idx_dtype) 

    fn = getattr(_sparsetools, self.format + '_matmat_pass1') 
    fn(M, N, 
     np.asarray(self.indptr, dtype=idx_dtype), 
     np.asarray(self.indices, dtype=idx_dtype), 
     np.asarray(other.indptr, dtype=idx_dtype), 
     np.asarray(other.indices, dtype=idx_dtype), 
     indptr) 

    nnz = indptr[-1] 

проход 2 выделяет indptr (различного размера), и на основе nnzindices и data.

indptr = np.asarray(indptr, dtype=idx_dtype) 
    indices = np.empty(nnz, dtype=idx_dtype) 
    data = np.empty(nnz, dtype=upcast(self.dtype, other.dtype)) 

    fn = getattr(_sparsetools, self.format + '_matmat_pass2') 
    fn(M, N, np.asarray(self.indptr, dtype=idx_dtype), 
     np.asarray(self.indices, dtype=idx_dtype), 
     self.data, 
     np.asarray(other.indptr, dtype=idx_dtype), 
     np.asarray(other.indices, dtype=idx_dtype), 
     other.data, 
     indptr, indices, data) 

Последний изготовить новый массив, используя эти массивы.

return self.__class__((data,indices,indptr),shape=(M,N)) 

mkl библиотека должна использоваться таким же образом.

===================

https://github.com/scipy/scipy/blob/master/scipy/sparse/sparsetools/csr.h

имеет код С для csr_matmat_pass1 и csr_matmat_pass2

======= =============

В случае, если это помогает, вот чистая реализация этих проходов на Python. Дословный перевод без какой-либо попытки использовать любые операции с массивами.

def pass1(A, B): 
    nrow,ncol=A.shape 
    Aptr=A.indptr 
    Bptr=B.indptr 
    Cp=np.zeros(nrow+1,int) 
    mask=np.zeros(ncol,int)-1 
    nnz=0 
    for i in range(nrow): 
     row_nnz=0 
     for jj in range(Aptr[i],Aptr[i+1]): 
      j=A.indices[jj] 
      for kk in range(Bptr[j],Bptr[j+1]): 
       k=B.indices[kk] 
       if mask[k]!=i: 
        mask[k]=i 
        row_nnz += 1 
     nnz += row_nnz 
     Cp[i+1]=nnz 
    return Cp 

def pass2(A,B,Cnnz): 
    nrow,ncol=A.shape 
    Ap,Aj,Ax=A.indptr,A.indices,A.data 
    Bp,Bj,Bx=B.indptr,B.indices,B.data 

    next = np.zeros(ncol,int)-1 
    sums = np.zeros(ncol,A.dtype) 

    Cp=np.zeros(nrow+1,int) 
    Cj=np.zeros(Cnnz,int) 
    Cx=np.zeros(Cnnz,A.dtype) 
    nnz = 0 
    for i in range(nrow): 
     head = -2 
     length = 0 
     for jj in range(Ap[i], Ap[i+1]): 
      j, v = Aj[jj], Ax[jj] 
      for kk in range(Bp[j], Bp[j+1]): 
       k = Bj[kk] 
       sums[k] += v*Bx[kk] 
       if (next[k]==-1): 
        next[k], head=head, k 
        length += 1 
     print(i,sums, next) 
     for _ in range(length): 
      if sums[head] !=0: 
       Cj[nnz], Cx[nnz] = head, sums[head] 
       nnz += 1 
      temp = head 
      head = next[head] 
      next[temp], sums[temp] = -1, 0 
     Cp[i+1] = nnz    
    return Cp, Cj, Cx     

def pass0(A,B): 
    Cp = pass1(A,B) 
    nnz = Cp[-1] 
    Cp,Cj,Cx=pass2(A,B,nnz) 
    N,M=A.shape[0], B.shape[1] 
    C=sparse.csr_matrix((Cx, Cj, Cp), shape=(N,M)) 
    return C 

Оба A и B должны быть csr. Он использует транспонирование, для которого требуется преобразование, например. B = A.T.tocsr().