利用scipy稀疏阵列上的坐标有效地执行逐元素运算

我试图弄清楚如何在scipy稀疏数组（csc格式）上有效地执行以下操作：
元素伪代码：

try:
    M[i,j] /= (V[i] + V[j])
except ZeroDivisionError:
    pass # leave at old value

注意，这使得零元素为0，因此不增加输入矩阵M的密度。
对于稠密矩阵有一个简单的方法：

# V is a ndarray of shape (N), all values finite and >= 0.0
# M is a ndarray of shape (N,N), all values finite

rows_plus_cols = V[:, np.newaxis] + V[np.newaxis, :]
M = np.divide(M, rows_plus_cols, out=M, where=(rows_plus_cols != 0))

不幸的是，这种方法似乎不适用于稀疏矩阵，因为rows_plus_cols的中间计算会导致大小为(N,N)的密集矩阵。
典型的输入M可以是100 k乘100 k，具有~ 1 M个非零条目。显然，在实现一个完整的密集矩阵是不明智的政权。出于测试目的，以下是“足够接近”典型输入：

n = 100000
d = 10
# a.k.a. "mostly nonzero, but not always"
V = np.maximum((np.random.normal(size=n)+2), 0)

M = scipy.sparse.coo_array((np.random.normal(size=n*d), np.random.randint(n, size=(2, n*d))), shape=(n,n)).tocsc()
M.sum_duplicates()
M.prune()
M.sort_indices()

小问题：scipy稀疏数组似乎不允许使用where子句进行元素级除法。可以解决，但这不是主要问题。我知道您可以使用M.data访问scipy csc数组的底层数据点;我看不出如何在不具体化稠密矩阵的情况下进行作为坐标的函数的元素级计算。
建议？

这里有一种使用coo_array的方法。
在这个例子中，M是一个8x8 coo_array，V是一个长度为8的一维NumPy数组。如果M可能有重复的条目，则应在执行以下操作之前调用M.sum_duplicates()。

In [106]: M
Out[106]: 
<8x8 sparse array of type '<class 'numpy.float64'>'
    with 16 stored elements in COOrdinate format>

In [107]: M.toarray()
Out[107]: 
array([[0.  , 0.  , 0.44, 0.  , 0.93, 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 0.54, 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.2 , 0.21, 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.93, 0.76, 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.32, 0.88, 0.76],
       [0.91, 0.72, 0.  , 0.24, 0.87, 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.74, 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.  , 0.13, 0.  , 0.  , 0.  , 0.  ]])

In [108]: M.data
Out[108]: array([0.44, 0.93, 0.54, 0.2 , 0.21, 0.93, 0.76, 0.32, 0.88, 0.76, 0.91, 0.72, 0.24, 0.87, 0.74, 0.13])

In [109]: M.row
Out[109]: array([0, 0, 1, 2, 2, 3, 3, 4, 4, 4, 5, 5, 5, 5, 6, 7], dtype=int32)

In [110]: M.col
Out[110]: array([2, 4, 6, 4, 5, 5, 6, 5, 6, 7, 0, 1, 3, 4, 4, 3], dtype=int32)

In [111]: V
Out[111]: array([0.5, 1. , 2. , 0. , 0. , 0.5, 1. , 0. ])

首先将分母形成为V[M.row] + V[M.col]，并将任何出现的0替换为1，因此不是除以0，这些位置的值保持不变：

In [112]: denom = V[M.row] + V[M.col]

In [113]: denom[denom == 0] = 1

In [114]: denom
Out[114]: array([2.5, 0.5, 2. , 2. , 2.5, 0.5, 1. , 0.5, 1. , 1. , 1. , 1.5, 0.5, 0.5, 1. , 1. ])

我们将新数组命名为Q。Q_data将是Q的data数组：

In [115]: Q_data = M.data / denom

In [116]: Q_data
Out[116]: array([0.176, 1.86 , 0.27 , 0.1  , 0.084, 1.86 , 0.76 , 0.64 , 0.88 , 0.76 , 0.91 , 0.48 , 0.48 , 1.74 , 0.74 , 0.13 ])

使用与M相同的形状和相同的行和列索引，用Q_data创建一个新的coo_array：

In [117]: Q = sparse.coo_array((Q_data, (M.row, M.col)), shape=M.shape)

In [118]: Q.toarray()
Out[118]: 
array([[0.   , 0.   , 0.176, 0.   , 1.86 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.27 , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.1  , 0.084, 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 1.86 , 0.76 , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.64 , 0.88 , 0.76 ],
       [0.91 , 0.48 , 0.   , 0.48 , 1.74 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.74 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.13 , 0.   , 0.   , 0.   , 0.   ]])

执行密集版本的计算，以验证其是否按预期工作：

In [119]: VV = np.add.outer(V, V)

In [120]: VV[VV == 0] = 1

In [121]: M.toarray() / VV
Out[121]: 
array([[0.   , 0.   , 0.176, 0.   , 1.86 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.   , 0.27 , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.1  , 0.084, 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 1.86 , 0.76 , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.   , 0.64 , 0.88 , 0.76 ],
       [0.91 , 0.48 , 0.   , 0.48 , 1.74 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.   , 0.74 , 0.   , 0.   , 0.   ],
       [0.   , 0.   , 0.   , 0.13 , 0.   , 0.   , 0.   , 0.   ]])

你可以将ipython会话的第112、113、115和117行打包成一个函数。

利用scipy稀疏阵列上的坐标有效地执行逐元素运算

1条答案

相关问题

热门标签

最新问答