具有广播的稀疏 Scipy 矩阵和向量的元素明智最大值答案

【问题标题】：Elementwise maximum of sparse Scipy matrix & vector with broadcasting具有广播的稀疏 Scipy 矩阵和向量的元素明智最大值
【发布时间】：2021-02-28 12:14:24
【问题描述】：

我需要一个快速的元素最大值，它将 n×m scipy 稀疏矩阵的每一行逐元素与稀疏的 1×m 矩阵进行比较。这在 Numpy 中使用 np.maximum(mat, vec) 通过 Numpy 的广播完美运行。

但是，Scipy 的.maximum() 没有广播。我的矩阵很大，所以我无法将其转换为 numpy 数组。

我目前的解决方法是使用mat[row,:].maximum(vec) 循环遍历多行垫子。这个大循环破坏了我的代码效率（必须多次执行）。我的慢解决方案在下面的第二个代码 sn-p 中——有更好的解决方案吗？

# 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 的回答，因为这个问题专门针对高性能解决方案，Max 的解决方案在我尝试的各种输入上提供了 1000x-2500x 的巨大加速，但代价是更多的代码行和 Numba 编译。但是，对于一般用途，Daniel F 的单线是一个很好的解决方案，在我尝试过的示例上提供了 10 到 50 倍的加速——我可能会用于许多其他事情。

【问题讨论】：

好吧，很奇怪，用csr_matrix 或coo_matrix 做这些会使我的内核崩溃。 csc_matrix 虽然有效。任何人都可以复制吗？
是的——刚刚编辑使它们成为 csc_matrix 以便该示例适用于其他人。我也有一些使用 csr_matrix 和 .maximum() 函数的崩溃（特别是段错误）。
好的，可能需要一些后端向导来弄清楚编译后的scipy 代码中发生了什么。请耐心等待，其中一个人可能需要一段时间才能看到这一点。
所以，据我所知，除了csc_matrix 之外，从scipy 到csr_matrix 的每一种稀疏矩阵形式都会因为maximum 而恢复为csr_matrix，所以至少这只是一个根本原因。我无法遵循cs_matrix 的_maxumum_minimum_ 函数，或者maximum 如何调用它。似乎if dense_check(other): 调用应该引发ValueError，因为它是一个数组。也许那里发生了短路
@max9111 回复：稀疏性，它变化但合理的值是非零元素的 0.001 到 0.005 分数。

标签： python numpy scipy sparse-matrix elementwise-operations

【解决方案1】：

scipy.sparse 矩阵不广播。完全没有。所以除非你能想出一些方法来操作indices 和inpts（我没有），否则你会被卡住。我能想到的最好的办法就是vstack 你的vecs，直到它们的形状与mat 相同。它似乎提供了很好的加速效果，尽管它并不能解释 csr 的段错误怪异。

#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]]

【讨论】：

这给了我问题中示例的 15 倍加速，更像是在更大 (50k X 10k) 矩阵上的 50 倍加速。感谢您的发表！我同意构建一个垫子大小的 vecs vstack 对我来说是“错误的”（效率方面），但也许没有办法解决这个问题。

【解决方案2】：

查看maximum 方法和代码，尤其是/scipy/sparse/compressed.py 中的_binopt 方法，很明显它可以与标量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>

失败是另一个是一维稀疏矩阵：

In [65]: mat.maximum(mat[0,:])
Segmentation fault (core dumped)

mat.maximum(mat[:,0]) 运行没有错误，但我不确定这些值。 mat[:,0] 将具有相同的大小 indptr。

如果mat 是csc，我认为mat.maximum(mat[:,0]) 会给出同样的错误，但事实并非如此。

说实话，这种运算对于稀疏矩阵来说并不是一个强项。它的数学核心是矩阵乘法。这就是它们最初开发的目的 - 稀疏线性代数问题，例如有限差分和有限元。

【讨论】：

谢谢。我正在使用 scipy.sparse 矩阵进行矩阵乘法。这部分代码显然非常快（比 numpy 数组快得多）。但是，其中一个点积需要作为 (mat, vec) 的元素最大值的矩阵。
vec 的非零值是否大大增加了结果的非零值？如果vec 在与mat 相同的位置有零，那么maximum 可能比如果可以将mat 中的0 变为非零要简单得多。相对而言，更改 csr/csc 的稀疏性是昂贵的。
我也有同样的想法！是的，vec 中的零也是 mat 中的零（即 0 列），并且我已经从不同的操作中获得了这些零的索引。我写了 numpy 版本，它比我以前的 numpy 版本快得多。仍然不确定如何编写最佳的 scipy 版本（我对 scipy 比较陌生，所以仍在弄清楚如何有效地做事）
numpy 的方式是：`np.maximum( mat[:,nonzero_cols_booleanvec] , vec[nonzero_cols_boolean_vec] )' 现在我正在尝试一个在 scipy 中进行切片的版本，然后转换为 numpy数组（这使得 mat 和 vec 比完整 mat 和完整 vec 的密集版本更窄=更少内存，因为它们是稀疏的，并且我们只在切片后保留非零），然后使用 np.maximum。没有这种丑陋的演员阵容可能会更好，但我不知道。

【解决方案3】：

低级方法

与往常一样，您可以考虑如何为该操作构建适当的稀疏矩阵格式，对于 csr 矩阵，主要组件是 shape、data_arr、indices 和 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 MB 的步长增加了输出的大小。在创建 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)

【讨论】：

这太棒了。在更大的输入上，我得到了约 2000 倍的加速。昨天我第一次尝试在 Numba 中写这个（以前从未使用过 Numba），结果比原来的要慢得多——显然这是一个技巧问题。问题：基于我非常有限的 C++ 知识，我能否通过对 mat 行进行多线程处理来更快，因为每一行都是独立的？
@cataclysmic 是的，这应该是可能的，但我认为复杂的 C++ 知识是必要的，才能获得真正高效的代码。上面 Numba 实现的主要问题不是缺乏并行化，而是在 Numba-Lists 中收集结果。也许我会在晚上找到一些时间来优化这个（应该再给两个或更多的因素）
太棒了，但我想我一定是做错了什么——我无法复制你的时间安排。在编译原始版本和优化版本后，我得到原始版本为 0.0643 秒，优化版本为 0.0488 秒，这更好，但不是优化与原始版本的 5 倍加速。我看到您的代码正在调用 sparse_elementwise_maximum_wrap()，但我从上面的原始版本中看到了 sparse_elementwise_maximum_2_wrap——这个函数可能会改变吗？再次感谢这些惊人的实现。
另外，我应该在生成 mat 和 vec 之前在原始代码中添加对 np.random.seed(12345) 的调用，以便它们在时序一致性方面始终相同。
@cataclysmic memstep 太小了。这个值通常应该适用于现实世界的问题。两个函数的包装器保持不变。