Ускоренный способ построения матрицы ожидаемых частот из разреженной матрицы отсчетов

Question

У меня есть матрица сжатого разреженного ряда, содержащая подсчеты. Я хочу построить матрицу, содержащую ожидаемые частоты для этих отсчетов. Вот код, который я использую сейчас:

from scipy.sparse import coo_matrix 

#m is a csr_matrix 

col_total = m.sum(axis=0) 
row_total = m.sum(axis=1) 
n = int(col_total.sum(axis=1)) 
A = coo_matrix(m) 

for i,j in zip(A.row,A.col): 
    m[i,j]= col_total.item(j)*row_total.item(i)/n 

Это прекрасно работает на небольшой матрице. На большей матрице (> 1 Гб) цикл for занимает несколько дней. Есть ли способ сделать это быстрее?

Answer 1

Чтобы немного расширить на @hpaulj «s ответа, вы можете избавиться от петли for пути создания выходной матрицы непосредственно от ожидаемых частот и строк/индексы столбцов ненулевых элементов в m:

from scipy import sparse 
import numpy as np 

def fast_efreqs(m): 

    col_total = np.array(m.sum(axis=0)).ravel() 
    row_total = np.array(m.sum(axis=1)).ravel() 

    # I'm casting this to an int for consistency with your version, but it's 
    # not clear to me why you would want to do this... 
    grand_total = int(col_total.sum()) 

    ridx, cidx = m.nonzero()   # indices of non-zero elements in m 
    efreqs = row_total[ridx] * col_total[cidx]/grand_total 

    return sparse.coo_matrix((efreqs, (ridx, cidx))) 

Для сравнения, вот ваш исходный код в функции:

def orig_efreqs(m): 

    col_total = m.sum(axis=0) 
    row_total = m.sum(axis=1) 
    n = int(col_total.sum(axis=1)) 

    A = sparse.coo_matrix(m) 
    for i,j in zip(A.row,A.col): 
     m[i,j]= col_total.item(j)*row_total.item(i)/n 

    return m 

Тест эквивалентность на маленькую матрице:

m = sparse.rand(100, 100, density=0.1, format='csr') 
print((orig_efreqs(m.copy()) != fast_efreqs(m)).nnz == 0) 
# True 

производительность Benchmark на большую матрицу:

In [1]: %%timeit m = sparse.rand(10000, 10000, density=0.01, format='csr') 
    .....: orig_efreqs(m) 
    .....: 
1 loops, best of 3: 2min 25s per loop 

In [2]: %%timeit m = sparse.rand(10000, 10000, density=0.01, format='csr') 
    .....: fast_efreqs(m) 
    .....: 
10 loops, best of 3: 38.3 ms per loop 

Answer 2

m.data = (col_total[:,A.col].A*(row_total[A.row,:].T.A)/n)[0] - полностью векторизованный способ расчета m.data. Вероятно, его можно немного почистить. col_total - matrix, поэтому для выполнения умножения элементов по элементу требуется некоторый дополнительный синтаксис.

Я продемонстрирую:

In [37]: m=sparse.rand(10,10,.1,'csr') 
In [38]: col_total=m.sum(axis=0) 
In [39]: row_total=m.sum(axis=1) 
In [40]: n=int(col_total.sum(axis=1)) 

In [42]: A=m.tocoo() 

In [46]: for i,j in zip(A.row,A.col): 
    ....:   m[i,j]= col_total.item(j)*row_total.item(i)/n 
    ....:  

In [49]: m.data 
Out[49]: 
array([ 0.39490171, 0.64246488, 0.19310878, 0.13847277, 0.2018023 , 
     0.008504 , 0.04387622, 0.10903026, 0.37976005, 0.11414632]) 

In [51]: col_total[:,A.col].A*(row_total[A.row,:].T.A)/n 
Out[51]: 
array([[ 0.39490171, 0.64246488, 0.19310878, 0.13847277, 0.2018023 , 
     0.008504 , 0.04387622, 0.10903026, 0.37976005, 0.11414632]]) 

In [53]: (col_total[:,A.col].A*(row_total[A.row,:].T.A)/n)[0] 
Out[53]: 
array([ 0.39490171, 0.64246488, 0.19310878, 0.13847277, 0.2018023 , 
     0.008504 , 0.04387622, 0.10903026, 0.37976005, 0.11414632])