#python #numpy #scipy #sparse-matrix #elementwise-operations
#python #numpy #scipy #разреженная матрица #поэлементные операции
Вопрос:
Мне нужен быстрый поэлементный максимум, который сравнивает каждую строку разреженной матрицы scipy размером n на m по элементам с разреженной матрицей размером 1 на m. Это отлично работает в Numpy с использованием np.maximum(mat, vec) широковещательной передачи via Numpy.
Однако у Scipy .maximum() нет широковещательной передачи. Моя матрица большая, поэтому я не могу преобразовать ее в массив numpy.
Мой текущий обходной путь заключается в циклическом переборе множества строк mat mat[row,:].maximum(vec) . Этот большой цикл снижает эффективность моего кода (это нужно делать много раз). Мое медленное решение находится во втором фрагменте кода ниже — есть ли лучшее решение?
# Example
import numpy as np
from scipy import sparse
mat = sparse.csc_matrix(np.arange(12).reshape((4,3)))
vec = sparse.csc_matrix([-1, 5, 100])
# Numpy's np.maximum() gives the **desired result** using broadcasting (but it can't handle sparse matrices):
numpy_result = np.maximum( mat.toarray(), vec.toarray() )
print( numpy_result )
# [[ 0 5 100]
# [ 3 5 100]
# [ 6 7 100]
# [ 9 10 100]]
# Scipy only compares the top row of mat to vec (no broadcasting!):
scipy_result = mat.maximum(vec)
print( scipy_result.toarray() )
# [[ 0 5 100]
# [ 3 4 5]
# [ 6 7 8]
# [ 9 10 11]]
#Reversing the order of mat and vec in the call to vec.maximum(mat) results in a single row output, and also frequently seg faults (!):
Увеличенный пример и текущее решение для тестирования скорости
import numpy as np
from scipy import sparse
import timeit
mat = sparse.csc_matrix( sparse.random(20000, 4000, density=.01, data_rvs=lambda s: np.random.randint(0, 5000, size=s)) )
vec = sparse.csc_matrix( sparse.random(1, 4000, density=.01, data_rvs=lambda s: np.random.randint(0, 5000, size=s)) )
def sparse_elementwise_maximum(mat, vec):
output = sparse.lil_matrix(mat.shape)
for row_idx in range( mat.shape[0] ):
output[row_idx] = mat[row_idx,:].maximum(vec)
return output
# Time it
num_timing_loops = 3.0
starttime = timeit.default_timer()
for _ in range(int(num_timing_loops)):
sparse_elementwise_maximum(mat, vec)
print('time per call is:', (timeit.default_timer() - starttime)/num_timing_loops, 'seconds')
# 15 seconds per call (way too slow!)
Редактировать
Я принимаю ответ Макса, поскольку вопрос касался конкретно высокопроизводительного решения, а решение Max предлагает огромные ускорения в 1000-2500 раз для различных входных данных, которые я пробовал, за счет большего количества строк кода и компиляции Numba. Однако для общего использования однострочник Даниэля Ф. — отличное решение, предлагающее ускорение в 10-50 раз на примерах, которые я пробовал — я, вероятно, буду использовать для многих других вещей.
Комментарии:
1. Хорошо, странно, делать это с
csr_matrixилиcoo_matrixсбой моего ядра.csc_matrixработает, хотя. Кто-нибудь может повторить?2. Да — просто отредактировано, чтобы сделать их csc_matrix, чтобы пример мог работать для других. У меня также были некоторые сбои при использовании csr_matrix с функцией .maximum() (в частности, ошибки seg).
3. Хорошо, возможно, потребуется несколько внутренних мастеров, чтобы выяснить, что происходит в скомпилированном
scipyкоде. Наберитесь терпения, одному из них может потребоваться некоторое время, чтобы увидеть это.4. Итак, насколько я могу судить, каждая форма разреженной матрицы
scipyотличается отcsc_matrixвозврата кcsr_matrixfor doingmaximum, так что, по крайней мере, это только одна основная причина. Я не могу следовать функцииcs_matrix‘s_maxumum_minimum_, или какmaximumона ее вызывает. Похоже, чтоif dense_check(other):этот вызов должен вызывать aValueError, поскольку это массив. Может быть, там что-то закорочено5. @max9111 Re: разреженность, она варьируется, но разумные значения составляют от 0,001 до 0,005 доли ненулевых элементов.
Ответ №1:
Низкоуровневый подход
Как всегда, вы можете подумать о том, как создается правильный формат разреженной матрицы для этой операции, для csr-матриц основными компонентами являются shape, data_arr, indexes и ind_ptr. С этими частями объекта scipy.sparse.csr это довольно просто, но, возможно, требует немного времени для реализации эффективного алгоритма на скомпилированном языке (C, C , Cython, Python-Numba). В своей реализации я использовал Numba, но перенос его на C должен быть легко возможен (изменения синтаксиса) и, возможно, избежать нарезки.
Реализация (первая попытка)
import numpy as np
import numba as nb
# get all needed components of the csr object and create a resulting csr object at the end
def sparse_elementwise_maximum_wrap(mat,vec):
mat_csr=mat.tocsr()
vec_csr=vec.tocsr()
shape_mat=mat_csr.shape
indices_mat=mat_csr.indices
indptr_mat=mat_csr.indptr
data_mat=mat_csr.data
indices_vec=vec_csr.indices
data_vec=vec_csr.data
res=sparse_elementwise_maximum_nb(indices_mat,indptr_mat,data_mat,shape_mat,indices_vec,data_vec)
res=sparse.csr_matrix(res, shape=shape_mat)
return res
@nb.njit(cache=True)
def sparse_elementwise_maximum_nb(indices_mat,indptr_mat,data_mat,shape_mat,vec_row_ind,vec_row_data):
data_res=[]
indices_res=[]
indptr_mat_res=[]
indptr_mat_=0
indptr_mat_res.append(indptr_mat_)
for row_idx in range(shape_mat[0]):
mat_row_ind=indices_mat[indptr_mat[row_idx]:indptr_mat[row_idx 1]]
mat_row_data=data_mat[indptr_mat[row_idx]:indptr_mat[row_idx 1]]
mat_ptr=0
vec_ptr=0
while mat_ptr<mat_row_ind.shape[0] and vec_ptr<vec_row_ind.shape[0]:
ind_mat=mat_row_ind[mat_ptr]
ind_vec=vec_row_ind[vec_ptr]
#value for both matrix and vector is present
if ind_mat==ind_vec:
data_res.append(max(mat_row_data[mat_ptr],vec_row_data[vec_ptr]))
indices_res.append(ind_mat)
mat_ptr =1
vec_ptr =1
indptr_mat_ =1
#only value for the matrix is present vector is assumed 0
elif ind_mat<ind_vec:
if mat_row_data[mat_ptr] >0:
data_res.append(mat_row_data[mat_ptr])
indices_res.append(ind_mat)
indptr_mat_ =1
mat_ptr =1
#only value for the vector is present matrix is assumed 0
else:
if vec_row_data[vec_ptr] >0:
data_res.append(vec_row_data[vec_ptr])
indices_res.append(ind_vec)
indptr_mat_ =1
vec_ptr =1
for i in range(mat_ptr,mat_row_ind.shape[0]):
if mat_row_data[i] >0:
data_res.append(mat_row_data[i])
indices_res.append(mat_row_ind[i])
indptr_mat_ =1
for i in range(vec_ptr,vec_row_ind.shape[0]):
if vec_row_data[i] >0:
data_res.append(vec_row_data[i])
indices_res.append(vec_row_ind[i])
indptr_mat_ =1
indptr_mat_res.append(indptr_mat_)
return np.array(data_res),np.array(indices_res),np.array(indptr_mat_res)
Реализация (оптимизирована)
При таком подходе списки заменяются динамически изменяемым массивом. Я увеличил размер выходных данных с шагом 60 МБ. При создании csr-объекта также не создается копия данных, только ссылки. Если вы хотите избежать перегрузки памяти, вам нужно скопировать массивы в конце.
@nb.njit(cache=True)
def sparse_elementwise_maximum_nb(indices_mat,indptr_mat,data_mat,shape_mat,vec_row_ind,vec_row_data):
mem_step=5_000_000
#preallocate memory for 5M non-zero elements (60 MB in this example)
data_res=np.empty(mem_step,dtype=data_mat.dtype)
indices_res=np.empty(mem_step,dtype=np.int32)
data_res_p=0
indptr_mat_res=np.empty((shape_mat[0] 1),dtype=np.int32)
indptr_mat_res[0]=0
indptr_mat_res_p=1
indptr_mat_=0
for row_idx in range(shape_mat[0]):
mat_row_ind=indices_mat[indptr_mat[row_idx]:indptr_mat[row_idx 1]]
mat_row_data=data_mat[indptr_mat[row_idx]:indptr_mat[row_idx 1]]
#check if resizing is necessary
if data_res.shape[0]<data_res_p shape_mat[1]:
#add at least memory for another mem_step elements
size_to_add=mem_step
if shape_mat[1] >size_to_add:
size_to_add=shape_mat[1]
data_res_2 =np.empty(data_res.shape[0] size_to_add,data_res.dtype)
indices_res_2=np.empty(indices_res.shape[0] size_to_add,indices_res.dtype)
for i in range(data_res_p):
data_res_2[i]=data_res[i]
indices_res_2[i]=indices_res[i]
data_res=data_res_2
indices_res=indices_res_2
mat_ptr=0
vec_ptr=0
while mat_ptr<mat_row_ind.shape[0] and vec_ptr<vec_row_ind.shape[0]:
ind_mat=mat_row_ind[mat_ptr]
ind_vec=vec_row_ind[vec_ptr]
#value for both matrix and vector is present
if ind_mat==ind_vec:
data_res[data_res_p]=max(mat_row_data[mat_ptr],vec_row_data[vec_ptr])
indices_res[data_res_p]=ind_mat
data_res_p =1
mat_ptr =1
vec_ptr =1
indptr_mat_ =1
#only value for the matrix is present vector is assumed 0
elif ind_mat<ind_vec:
if mat_row_data[mat_ptr] >0:
data_res[data_res_p]=mat_row_data[mat_ptr]
indices_res[data_res_p]=ind_mat
data_res_p =1
indptr_mat_ =1
mat_ptr =1
#only value for the vector is present matrix is assumed 0
else:
if vec_row_data[vec_ptr] >0:
data_res[data_res_p]=vec_row_data[vec_ptr]
indices_res[data_res_p]=ind_vec
data_res_p =1
indptr_mat_ =1
vec_ptr =1
for i in range(mat_ptr,mat_row_ind.shape[0]):
if mat_row_data[i] >0:
data_res[data_res_p]=mat_row_data[i]
indices_res[data_res_p]=mat_row_ind[i]
data_res_p =1
indptr_mat_ =1
for i in range(vec_ptr,vec_row_ind.shape[0]):
if vec_row_data[i] >0:
data_res[data_res_p]=vec_row_data[i]
indices_res[data_res_p]=vec_row_ind[i]
data_res_p =1
indptr_mat_ =1
indptr_mat_res[indptr_mat_res_p]=indptr_mat_
indptr_mat_res_p =1
return data_res[:data_res_p],indices_res[:data_res_p],indptr_mat_res
Максимальная память, выделенная в начале
Производительность и удобство использования этого подхода сильно зависят от входных данных. При таком подходе выделяется максимальная память (это может легко вызвать ошибки нехватки памяти).
@nb.njit(cache=True)
def sparse_elementwise_maximum_nb(indices_mat,indptr_mat,data_mat,shape_mat,vec_row_ind,vec_row_data,shrink_to_fit):
max_non_zero=shape_mat[0]*vec_row_data.shape[0] data_mat.shape[0]
data_res=np.empty(max_non_zero,dtype=data_mat.dtype)
indices_res=np.empty(max_non_zero,dtype=np.int32)
data_res_p=0
indptr_mat_res=np.empty((shape_mat[0] 1),dtype=np.int32)
indptr_mat_res[0]=0
indptr_mat_res_p=1
indptr_mat_=0
for row_idx in range(shape_mat[0]):
mat_row_ind=indices_mat[indptr_mat[row_idx]:indptr_mat[row_idx 1]]
mat_row_data=data_mat[indptr_mat[row_idx]:indptr_mat[row_idx 1]]
mat_ptr=0
vec_ptr=0
while mat_ptr<mat_row_ind.shape[0] and vec_ptr<vec_row_ind.shape[0]:
ind_mat=mat_row_ind[mat_ptr]
ind_vec=vec_row_ind[vec_ptr]
#value for both matrix and vector is present
if ind_mat==ind_vec:
data_res[data_res_p]=max(mat_row_data[mat_ptr],vec_row_data[vec_ptr])
indices_res[data_res_p]=ind_mat
data_res_p =1
mat_ptr =1
vec_ptr =1
indptr_mat_ =1
#only value for the matrix is present vector is assumed 0
elif ind_mat<ind_vec:
if mat_row_data[mat_ptr] >0:
data_res[data_res_p]=mat_row_data[mat_ptr]
indices_res[data_res_p]=ind_mat
data_res_p =1
indptr_mat_ =1
mat_ptr =1
#only value for the vector is present matrix is assumed 0
else:
if vec_row_data[vec_ptr] >0:
data_res[data_res_p]=vec_row_data[vec_ptr]
indices_res[data_res_p]=ind_vec
data_res_p =1
indptr_mat_ =1
vec_ptr =1
for i in range(mat_ptr,mat_row_ind.shape[0]):
if mat_row_data[i] >0:
data_res[data_res_p]=mat_row_data[i]
indices_res[data_res_p]=mat_row_ind[i]
data_res_p =1
indptr_mat_ =1
for i in range(vec_ptr,vec_row_ind.shape[0]):
if vec_row_data[i] >0:
data_res[data_res_p]=vec_row_data[i]
indices_res[data_res_p]=vec_row_ind[i]
data_res_p =1
indptr_mat_ =1
indptr_mat_res[indptr_mat_res_p]=indptr_mat_
indptr_mat_res_p =1
if shrink_to_fit==True:
data_res=np.copy(data_res[:data_res_p])
indices_res=np.copy(indices_res[:data_res_p])
else:
data_res=data_res[:data_res_p]
indices_res=indices_res[:data_res_p]
return data_res,indices_res,indptr_mat_res
# get all needed components of the csr object and create a resulting csr object at the end
def sparse_elementwise_maximum_wrap(mat,vec,shrink_to_fit=True):
mat_csr=mat.tocsr()
vec_csr=vec.tocsr()
shape_mat=mat_csr.shape
indices_mat=mat_csr.indices
indptr_mat=mat_csr.indptr
data_mat=mat_csr.data
indices_vec=vec_csr.indices
data_vec=vec_csr.data
res=sparse_elementwise_maximum_nb(indices_mat,indptr_mat,data_mat,shape_mat,indices_vec,data_vec,shrink_to_fit)
res=sparse.csr_matrix(res, shape=shape_mat)
return res
Тайминги
Numba имеет накладные расходы на компиляцию или некоторые накладные расходы для загрузки функции из кэша. Не рассматривайте первый вызов, если вы хотите получить время выполнения, а не компиляцию время выполнения.
import numpy as np
from scipy import sparse
mat = sparse.csr_matrix( sparse.random(20000, 4000, density=.01, data_rvs=lambda s: np.random.randint(0, 5000, size=s)) )
vec = sparse.csr_matrix( sparse.random(1, 4000, density=.01, data_rvs=lambda s: np.random.randint(0, 5000, size=s)) )
%timeit output=sparse_elementwise_maximum(mat, vec)
#for csc input
37.9 s ± 224 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
#for csr input
10.7 s ± 90.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
#Daniel F
%timeit sparse_maximum(mat, vec)
164 ms ± 1.74 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
#low level implementation (first try)
%timeit res=sparse_elementwise_maximum_wrap(mat,vec)
89.7 ms ± 2.51 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
#low level implementation (optimized, csr)
%timeit res=sparse_elementwise_maximum_wrap(mat,vec)
16.5 ms ± 122 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
#low level implementation (preallocation, without copying at the end)
%timeit res=sparse_elementwise_maximum_wrap(mat,vec)
16.5 ms ± 122 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
#low level implementation (preallocation, with copying at the end)
%timeit res=sparse_elementwise_maximum_wrap(mat,vec)
16.5 ms ± 122 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit res=sparse_elementwise_maximum_wrap(mat,vec,shrink_to_fit=False)
14.9 ms ± 110 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit res=sparse_elementwise_maximum_wrap(mat,vec,shrink_to_fit=True)
21.7 ms ± 399 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
#For comparison, copying the result takes
%%timeit
np.copy(res.data)
np.copy(res.indices)
np.copy(res.indptr)
7.8 ms ± 47.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Комментарии:
1. Это потрясающе. На больших входных данных я получаю ускорение ~ 2000x. Вчера я впервые попытался написать это в Numba (никогда раньше не использовал Numba) и в итоге получил что-то намного медленнее, чем оригинал — очевидно, это был вопрос мастерства. Вопрос: основываясь на моих очень ограниченных знаниях C , могу ли я работать еще быстрее, выполняя многопоточность по строкам mat, поскольку каждая строка независима?
2. @cataclysmic Да, это должно быть возможно, но я думаю, что для получения действительно эффективного кода необходимы сложные знания C . Основная проблема реализации Numba, описанной выше, заключается не в отсутствии распараллеливания, а в сборе результатов в Numba-списках. Может быть, я найду немного времени вечером, чтобы оптимизировать это (должно дать еще один коэффициент в два или более)
3. Удивительно, но я думаю, что я, должно быть, делаю что-то не так — я не могу воспроизвести ваши тайминги. После компиляции как исходной, так и оптимизированной версий я получаю 0,0643 секунды для оригинала и 0,0488 секунды для оптимизированной, что лучше, но не 5-кратное ускорение, которое вы получаете для оптимизированной и оригинальной. Я вижу, что ваш код вызывает sparse_elementwise_maximum_wrap(), но я вижу sparse_elementwise_maximum_2_wrap из исходной версии выше — возможно, эта функция изменилась? Еще раз спасибо, кстати, за эти удивительные реализации.
4. Также я должен был добавить в исходный код вызов np.random.seed(12345) непосредственно перед генерацией mat и vec, чтобы они всегда были одинаковыми для согласованности времени.
5. @cataclysmic Шаг памяти был слишком мал. Это значение обычно должно быть адаптировано к реальной проблеме. Оболочка остается неизменной для обеих функций.
Ответ №2:
scipy.sparse матрицы не транслируются. Вообще. Поэтому, если вы не сможете найти какой-либо способ работы с indices and inpts (у меня нет), вы застряли в стекировании. Лучшее, что я могу понять, это просто для vstack ваших vec s, пока они не будут иметь ту же форму, mat что и . Кажется, это дает хорошее ускорение, хотя и не объясняет странность с csr segfault.
#using `mat` and `vec` from the speed test
def sparse_maximum(mat, vec):
vec1 = sparse.vstack([vec for _ in range(mat.shape[0])])
return mat.maximum(vec1)
# Time it
num_timing_loops = 3.0
starttime = timeit.default_timer()
sparse_maximum(mat, vec)
print('time per call is:', (timeit.default_timer() - starttime)/num_timing_loops, 'seconds')
# I was getting 11-12 seconds on your original code
time per call is: 0.514533479333295 seconds
Доказательство того, что он работает с исходными матрицами:
vec = sparse.vstack([vec for _ in range(4)])
print(mat.maximum(vec).todense())
[[ 0 5 100]
[ 3 5 100]
[ 6 7 100]
[ 9 10 100]]
Комментарии:
1. Это дает мне 15-кратное ускорение в примере в моем вопросе и больше похоже на 50-кратное ускорение для большей (50k X 10k) матрицы. Спасибо за публикацию! Я согласен, что построение vstack vecs размером с мат кажется мне «неправильным» (с точки зрения эффективности), но, возможно, обойти это невозможно.
Ответ №3:
Глядя на maximum метод и код, особенно на _binopt метод, /scipy/sparse/compressed.py становится очевидным, что он может работать со скаляром other . Для разреженной other он создает новую разреженную матрицу (того же формата и формы), Используя indptr значения и т. Д. Если другой имеет ту же форму, он работает:
In [55]: mat = sparse.csr_matrix(np.arange(12).reshape((4,3)))
In [64]: mat.maximum(mat)
Out[64]:
<4x3 sparse matrix of type '<class 'numpy.int64'>'
with 11 stored elements in Compressed Sparse Row format>
Сбой заключается в том, что другая — 1d разреженная матрица:
In [65]: mat.maximum(mat[0,:])
Segmentation fault (core dumped)
mat.maximum(mat[:,0]) выполняется без ошибок, хотя я не уверен в значениях. mat[:,0] будет иметь тот же размер indptr .
Я думал mat.maximum(mat[:,0]) , что это даст ту же ошибку, если mat было csc , но это не так.
Давайте будем честными, такого рода операции не являются сильной стороной для разреженных матриц. Ядром его математики является умножение матриц. Это то, для чего они изначально были разработаны — разреженные задачи линейной алгебры, такие как конечная разность и конечный элемент.
Комментарии:
1. Спасибо. Я использую scipy.sparse матрицы для умножения матриц. Очевидно, что эта часть кода выполняется очень быстро (намного быстрее, чем массивы numpy). Однако для одного из точечных произведений необходима матрица, которая является поэлементным максимумом (mat, vec).
2.
vecСущественно ли увеличиваются ненулевые значения к ненулевым значениям результатов? Еслиvecимеет нули в том же месте,matчтоmaximumи это, потенциально намного проще, чем если может превратить 0 вmatненулевое значение. Изменения в разреженности acsr/cscстоят дорого, говоря относительно.3. У меня только что была та же идея! Да, нули в vec также являются нулями в mat (т. Е. 0-столбцы), и у меня уже есть индексы этих нулей из другой операции. Я написал версию numpy, которая намного быстрее, чем моя предыдущая версия numpy. Все еще не уверен, как написать оптимальную версию scipy (я относительно новичок в scipy, поэтому все еще пытаюсь понять, как делать что-то эффективно)
4. Способ numpy: `np.maximum(mat[:,nonzero_cols_booleanvec] , vec [nonzero_cols_boolean_vec] )’ Прямо сейчас я пробую версию, которая выполняет нарезку в scipy, а затем преобразует в массивы numpy (что делает mat и vec намного уже = меньше памяти, чем плотные версии full mat amp;полный vec, потому что они разреженные, и мы сохраняем ненулевые значения только после нарезки), затем использует np.maximum. Вероятно, есть что-то намного лучшее без этого уродливого кастинга, но я этого не знаю.