scipy.sparse矩阵：对非零元素减去行均值

Question

我有一个以csr_matrix格式表示的稀疏矩阵。对于每一行，我需要从非零元素中减去行均值。这些均值必须在行的非零元素数量上计算（而不是行的长度）。 我找到了一种快速计算行均值的方法，代码如下：

# M is a csr_matrix
sums = np.squeeze(np.asarray(M.sum(1)))    # sum of the nonzero elements, for each row
counts = np.diff(M.tocsr().indptr)         # count of the nonzero elements, for each row


# for the i-th row the mean is just sums[i] / float(counts[i]) 

问题在于更新部分。我需要一种快速的方法来完成这个任务。 实际上，我正在将M转换为lil_matrix，并以这种方式执行更新。

M = M.tolil()

for i in xrange(len(sums)):
    for j in M.getrow(i).nonzero()[1]:
        M[i, j] -= sums[i] / float(counts[i])

速度很慢。有没有更快的解决方案建议？

- revy

---

我会尝试使用np.repeat和counts来复制行均值，并直接从M.data数组中进行替换。 - hpaulj

Answer 1

这个有些棘手。我想我明白了。基本思路是我们尝试获得一个带有均值的对角矩阵，以及一个在M中非零数据位置处具有1的矩阵。然后我们将它们相乘并从M中减去积。以下是操作步骤...

>>> import numpy as np
>>> import scipy.sparse as sp
>>> a = sp.csr_matrix([[1., 0., 2.], [1.,2.,3.]])
>>> a.todense()
matrix([[ 1.,  0.,  2.],
        [ 1.,  2.,  3.]])
>>> tot = np.array(a.sum(axis=1).squeeze())[0]
>>> tot
array([ 3.,  6.])
>>> cts = np.diff(a.indptr)
>>> cts
array([2, 3], dtype=int32)
>>> mu = tot/cts
>>> mu
array([ 1.5,  2. ])
>>> d = sp.diags(mu, 0)
>>> d.todense()
matrix([[ 1.5,  0. ],
        [ 0. ,  2. ]])
>>> b = a.copy()
>>> b.data = np.ones_like(b.data)
>>> b.todense()
matrix([[ 1.,  0.,  1.],
        [ 1.,  1.,  1.]])
>>> (d * b).todense()
matrix([[ 1.5,  0. ,  1.5],
        [ 2. ,  2. ,  2. ]])
>>> (a - d*b).todense()
matrix([[-0.5,  0. ,  0.5],
        [-1. ,  0. ,  1. ]])

祝你好运！希望这可以帮到你。

- Dthal

---

是的，我也在考虑同样的事情。谢谢！ - revy

Answer 2

从 @Dthal 的示例开始：

In [92]: a = sparse.csr_matrix([[1.,0,2],[1,2,3]])
In [93]: a.A
Out[93]: 
array([[ 1.,  0.,  2.],
       [ 1.,  2.,  3.]])

In [94]: sums=np.squeeze(a.sum(1).A)
# sums=a.sum(1).A1   # shortcut
In [95]: counts=np.diff(a.tocsr().indptr)
In [96]: means=sums/counts
In [97]: sums
Out[97]: array([ 3.,  6.])
In [98]: counts
Out[98]: array([2, 3], dtype=int32)
In [99]: means
Out[99]: array([ 1.5,  2. ])

repeat命令可以复制means，生成一个与矩阵data大小相同的数组。

In [100]: mc = np.repeat(means, counts)
In [101]: mc
Out[101]: array([ 1.5,  1.5,  2. ,  2. ,  2. ])

这个mc和@Dthal的(b*d).data是相同的。

现在只需从data中减去它。

In [102]: a.data -= mc
In [103]: a.A
Out[103]: 
array([[-0.5,  0. ,  0.5],
       [-1. ,  0. ,  1. ]])

- hpaulj