首页 文章

如何通过广播的密集1d数组元素地乘以scipy.sparse矩阵?

提问于
浏览
41

假设我有一个2d稀疏数组 . 在我的实际用例中,行数和列数都要大得多(比如20000和50000),因此在使用密集表示时它无法适应内存:

>>> import numpy as np
>>> import scipy.sparse as ssp

>>> a = ssp.lil_matrix((5, 3))
>>> a[1, 2] = -1
>>> a[4, 1] = 2
>>> a.todense()
matrix([[ 0.,  0.,  0.],
        [ 0.,  0., -1.],
        [ 0.,  0.,  0.],
        [ 0.,  0.,  0.],
        [ 0.,  2.,  0.]])

现在假设我有一个密集的1d数组,其中包含大小为3的所有非零组件(在我的实际案例中为50000):

>>> d = np.ones(3) * 3
>>> d
array([ 3.,  3.,  3.])

我想使用numpy的常用广播语义来计算a和d的元素乘法 . 然而,scipy中的稀疏矩阵属于np.matrix:'*'运算符被重载使其行为类似于矩阵乘法而不是逐元乘法:

>>> a * d
array([ 0., -3.,  0.,  0.,  6.])

一种解决方案是将'a'切换到'*'运算符的数组语义,这将产生预期的结果:

>>> a.toarray() * d
array([[ 0.,  0.,  0.],
       [ 0.,  0., -3.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  6.,  0.]])

但是我不能这样做,因为对toarray()的调用将实现不适合内存的密集版本'a'(结果也会密集):

>>> ssp.issparse(a.toarray())
False

知道如何构建这个,同时只保留稀疏的数据结构,而不必在'a'的列上进行无效的python循环?

3 回答

  • 1

    我也在scipy.org上回复了,但我想我应该在这里添加一个答案,以防其他人在搜索时找到这个页面 .

    您可以将矢量转换为稀疏对角矩阵,然后使用矩阵乘法(使用*)来执行与广播相同的操作,但效率很高 .

    >>> d = ssp.lil_matrix((3,3))
    >>> d.setdiag(np.ones(3)*3)
    >>> a*d
    <5x3 sparse matrix of type '<type 'numpy.float64'>'
     with 2 stored elements in Compressed Sparse Row format>
    >>> (a*d).todense()
    matrix([[ 0.,  0.,  0.],
            [ 0.,  0., -3.],
            [ 0.,  0.,  0.],
            [ 0.,  0.,  0.],
            [ 0.,  6.,  0.]])
    

    希望有所帮助!

  • 44

    我认为A.multiply(B)应该在scipy稀疏中工作 . 方法乘法进行“逐点”乘法,而不是矩阵乘法 .

    HTH

  • 23

    好吧,这是一个简单的代码,可以做你想要的 . 我不知道它是否像你想的那样高效,所以接受或离开它:

    import scipy.sparse as ssp
    def pointmult(a,b):
        x = a.copy()
        for i in xrange(a.shape[0]):
            if x.data[i]:
                for j in xrange(len(x.data[i])):
                    x.data[i] *= b[x.rows[i]]
        return x
    

    它仅适用于lil矩阵,因此如果您希望它与其他格式一起使用,您将不得不进行一些更改 .

相关问题