使用带有稀疏 scipy 矩阵的广播答案

【问题标题】：Using broadcasting with sparse scipy matrices使用带有稀疏 scipy 矩阵的广播
【发布时间】：2017-03-31 19:54:17
【问题描述】：

我有一个numpy 数组Z，形状为(k,N)，第二个数组X，形状为(N,n)。

使用numpy广播，我可以轻松获得一个新数组H，其形状为(n,k,N)，其切片为数组Z，其行已乘以X的列：

H = Z.reshape((1, k, N)) * X.T.reshape((n, 1, N))

这很好用，而且速度惊人。现在，X 非常稀疏，我想使用稀疏矩阵运算进一步加快此运算。

但是，如果我执行以下操作：

import scipy.sparse as sprs
spX = sprs.csr_matrix(X)
H = (Z.reshape((1,k,N))*spX.T.reshape((n,1,N))).dot(Z.T)

我收到以下错误：

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:\Python27\lib\site-packages\scipy\sparse\base.py", line 126, in reshape
    self.__class__.__name__)
NotImplementedError: Reshaping not implemented for csc_matrix.

有没有办法使用带有稀疏scipy 矩阵的广播？

【问题讨论】：

除了二维限制之外，稀疏数学仅适用于其他稀疏矩阵。使用密集数组时，稀疏矩阵将转换为密集矩阵。与spX 一起玩，看看各种数学运算会发生什么。

标签： python arrays numpy sparse-matrix

【解决方案1】：

Scipy 稀疏矩阵仅限于二维形状。但是您可以以“稀疏”的方式使用 Numpy：

H = np.zeros((n,k,N), np.result_type(Z, X))
I, J = np.nonzero(X)
Z_ = np.broadcast_to(Z, H.shape)
H[J,:,I] = Z_[J,:,I] * X[I,J,None]

不幸的是，H 的结果仍然是一个密集数组。

注意使用None 进行索引是在所需轴上添加单位长度维度的便捷方法。将高级索引与切片相结合时的结果顺序为explained in the docs。

【讨论】：