如何用广播的密集1d数组逐元素乘以scipy.sparse矩阵？

Question

假设我有一个二维稀疏数组。在我的实际应用中，行和列的数量都非常大（比如说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.]])

现在假设我有一个一维的密集数组，里面全是非零的元素，大小是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循环呢？

Answer 1

好吧，这里有一段简单的代码，可以实现你想要的功能。不过我不确定它的效率是否符合你的要求，所以你可以选择用或者不使用：

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

这段代码只适用于小矩阵，如果你想让它适用于其他格式，就需要做一些修改。

Answer 2

我觉得在scipy的稀疏矩阵中，A.multiply(B)应该是可以用的。这个multiply方法是进行“逐点”相乘，而不是矩阵相乘。

希望对你有帮助！

Answer 3

我在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.]])

希望这对你有帮助！