Эффективно вычислять столбцовую сумму разреженного массива, где каждый ненулевой элемент равен 1

Question

У меня есть куча данных в формате SciPy compressed sparse row (CSR). Конечно, большинство элементов равно нулю, и я также знаю, что все ненулевые элементы имеют значение 1. Я хочу вычислить суммы по различным подмножествам строк моей матрицы. На данный момент я делаю следующее:

import numpy as np 
import scipy as sp 
import scipy.sparse 

# create some data with sparsely distributed ones 
data = np.random.choice((0, 1), size=(1000, 2000), p=(0.95, 0.05)) 
data = sp.sparse.csr_matrix(data, dtype='int8') 

# generate column-wise sums over random subsets of rows 
nrand = 1000 
for k in range(nrand): 
    inds = np.random.choice(data.shape[0], size=100, replace=False) 

    # 60% of time is spent here 
    extracted_rows = data[inds] 

    # 20% of time is spent here 
    row_sum = extracted_rows.sum(axis=0) 

Последние нескольких строк есть узкое место в большем вычислительном трубопроводе. Как я отметил в коде, 60% времени потрачено нарезание данных из случайных индексов, а 20% потрачено на вычисление фактической суммы.

Мне кажется, что я должен использовать свои знания о данных в массиве (т. Е. Любое ненулевое значение в разреженной матрице будет равно 1, других значений нет) для более эффективного вычисления этих сумм. К сожалению, я не могу понять, как это сделать. Возможно, только с data.indices? Я пробовал другие структуры разреженности (например, CSC-матрицу), а также сначала преобразовывался в плотный массив, но эти подходы были все медленнее, чем этот матричный подход CSR.

Answer 1

Хорошо известно, что индексирование разреженных матриц относительно медленное. И есть вопросы о том, как обойти это, обратившись непосредственно к атрибутам данных.

Но сначала некоторые тайминги. Использование data и ind как вы показать мне получить

In [23]: datad=data.A # times at 3.76 ms per loop 

In [24]: timeit row_sumd=datad[inds].sum(axis=0) 
1000 loops, best of 3: 529 µs per loop 

In [25]: timeit row_sum=data[inds].sum(axis=0) 
1000 loops, best of 3: 890 µs per loop 

In [26]: timeit d=datad[inds] 
10000 loops, best of 3: 55.9 µs per loop 

In [27]: timeit d=data[inds] 
1000 loops, best of 3: 617 µs per loop 

Редкой версия работает медленнее, чем густые один, но не на много. Редкое индексирование происходит гораздо медленнее, но его сумма несколько выше.

Редкая сумма делается с матричным продуктом

def sparse.spmatrix.sum 
    .... 
    return np.asmatrix(np.ones((1, m), dtype=res_dtype)) * self 

Это наводит на мысль, что более быстрый способ - превратить inds в соответствующий массив 1s и размножаться.

In [49]: %%timeit 
    ....: b=np.zeros((1,data.shape[0]),'int8') 
    ....: b[:,inds]=1 
    ....: rowmul=b*data 
    ....: 
1000 loops, best of 3: 587 µs per loop 

Это делает разреженную операцию примерно такой же плотной, как и ее эквивалентная плотность. (Но преобразование в густых гораздо медленнее)

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

Тест последний раз отсутствует в np.asmatrix, который присутствует в разреженной sum , Но времена одинаковы, и результаты те же, и результаты те же самые, и результаты те же самые, и результаты те же самые, и результаты те же самые

In [232]: timeit b=np.zeros((1,data.shape[0]),'int8'); b[:,inds]=1; x1=np.asmatrix(b)*data 
1000 loops, best of 3: 661 µs per loop 

In [233]: timeit b=np.zeros((1,data.shape[0]),'int8'); b[:,inds]=1; x2=b*data 
1000 loops, best of 3: 605 µs per loop 

Производится матрица, другая - массив. Но оба делают матричный продукт, 2-й тусклый B против 1-го из data. Несмотря на то, что b является массивом, задача фактически делегирована data, а ее матричный продукт - не столь прозрачным способом.

In [234]: x1 
Out[234]: matrix([[9, 9, 5, ..., 9, 5, 3]], dtype=int8) 

In [235]: x2 
Out[235]: array([[9, 9, 5, ..., 9, 5, 3]], dtype=int8) 

b*data.A - умножение элементов и вызывает ошибку; np.dot(b,data.A) работает, но медленнее.

Новый numpy/python имеет оператора matmul.Я вижу то же самое время картина:

In [280]: timeit b@dataA   # dense product 
100 loops, best of 3: 2.64 ms per loop 

In [281]: timeit b@data.A   # slower due to `.A` conversion 
100 loops, best of 3: 6.44 ms per loop 

In [282]: timeit b@data    # sparse product 
1000 loops, best of 3: 571 µs per loop 

np.dot может также делегировать действия по sparse, хотя вы должны быть осторожны. Я просто повесил свою машину с помощью np.dot(csr_matrix(b),data.A).

Answer 2

Вот Векторизованный подход после преобразования data в плотный массив, а также получать все те inds в векторизованном способе с использованием argpartition-based method -

# Number of selections as a parameter 
n = 100 

# Get inds across all iterations in a vectorized manner as a 2D array. 
inds2D = np.random.rand(nrand,data.shape[0]).argpartition(n)[:,:n] 

# Index into data with those 2D array indices. Then, convert to dense NumPy array, 
# reshape and sum reduce to get the final output 
out = np.array(data.todense())[inds2D.ravel()].reshape(nrand,n,-1).sum(1) 

Продолжительности испытание -

определение 1) Функции:

def org_app(nrand,n): 
    out = np.zeros((nrand,data.shape[1]),dtype=int) 
    for k in range(nrand): 
     inds = np.random.choice(data.shape[0], size=n, replace=False)   
     extracted_rows = data[inds]  
     out[k] = extracted_rows.sum(axis=0)  
    return out 


def vectorized_app(nrand,n): 
    inds2D = np.random.rand(nrand,data.shape[0]).argpartition(n)[:,:n] 
    return np.array(data.todense())[inds2D.ravel()].reshape(nrand,n,-1).sum(1) 

Сроки:

In [205]: # create some data with sparsely distributed ones 
    ...: data = np.random.choice((0, 1), size=(1000, 2000), p=(0.95, 0.05)) 
    ...: data = sp.sparse.csr_matrix(data, dtype='int8') 
    ...: 
    ...: # generate column-wise sums over random subsets of rows 
    ...: nrand = 1000 
    ...: n = 100 
    ...: 

In [206]: %timeit org_app(nrand,n) 
1 loops, best of 3: 1.38 s per loop 

In [207]: %timeit vectorized_app(nrand,n) 
1 loops, best of 3: 826 ms per loop