如何通过索引将 scipy.sparse 矩阵分配给 NumPy 数组？

Posted 2023-02-16

【中文标题】如何通过索引将 scipy.sparse 矩阵分配给 NumPy 数组？

【英文标题】：How to assign scipy.sparse matrix to NumPy array via indexing?

【发布时间】：2014-11-14 16:18:24

【问题描述】：

当我尝试将 scipy.sparse 矩阵 s（任何可用的稀疏类型）分配给 NumPy 数组 a 时，如下所示：

a[:] = s

我收到了TypeError：

TypeError: float() 参数必须是字符串或数字

有没有办法解决这个问题？

我知道todense() 和toarray() 方法，但我真的很想避免不必要的复制，而且我更愿意对NumPy 数组和SciPy 稀疏矩阵使用相同的代码。 目前，我并不担心从稀疏矩阵中获取值效率低下。

是否有某种围绕稀疏矩阵的包装器可以与 NumPy 索引分配一起使用？

如果没有，有什么建议我可以自己构建这样的东西吗？

在这种情况下是否有不同的稀疏数组库与 NumPy 配合？

更新：

我在 NumPy 源中四处寻找，搜索错误消息字符串，我想我在函数 @TYPE@_setitem() 中找到了在 numpy/core/src/multiarray/arraytypes.c.src around line 187 中发生索引分配的部分。

我仍然不明白，但在某些时候，float() 函数似乎被调用（如果a 是一个浮点数组）。因此，我尝试对其中一个 SciPy 稀疏矩阵类进行猴子修补，以允许调用此函数：

import scipy
s = scipy.sparse.dok_matrix((5, 1))
def myfloat(self):
    assert self.shape == (1, 1)
    return self[0, 0]
scipy.sparse.dok.dok_matrix.__float__ = myfloat
a[:] = s

遗憾的是，这不起作用，因为 float() 是在整个稀疏矩阵上调用的，而不是在其中的单个项目上。

所以我想我的新问题是：如何进一步更改稀疏矩阵类以使 NumPy 遍历所有项目并在每个项目上调用 float()？

另一个更新：

我在 Github (https://github.com/FRidh/sparse) 上找到了一个稀疏数组模块，它允许分配给 NumPy 数组。遗憾的是，该模块的功能非常有限（例如，切片还没有真正起作用），但它可能有助于理解如何分配给 NumPy 数组。 我会进一步调查...

还有一个更新：

我做了更多的挖掘，发现一个更有趣的源文件可能是numpy/core/src/multiarray/ctors.c。 我怀疑函数 PySequence_Check() (docs/code) 在分配期间的某个时候被调用。 https://github.com/FRidh/sparse 中的简单稀疏数组类通过了测试，但看起来 SciPy 中的稀疏矩阵类没有通过测试（尽管在我看来它们是序列）。

检查了__array_struct__、__array_interface__ 和__array__，然后以某种方式确定它们不是序列。 属性__getitem__ 和__len__（所有稀疏数组类都有！）未被检查。

这让我想到了另一个问题：我如何以传递PySequence_Check() 的方式操作稀疏矩阵类（或其对象）？

我认为一旦它们被识别为序列，分配就应该起作用，因为 __getitem__() 和 __len__() 应该足够了。

【问题讨论】：

什么是s？ a 是如何创建的？什么稀疏格式？数据类型？

例如，像这样：s = scipy.sparse.dok_matrix((4, 5))，所以它是 Dictionary Of Keys 格式，dtype 是 float64。我还尝试了其他稀疏格式，但它们也不起作用。如果我在a 中使用整数dtype，则错误消息略有不同：“TypeError: long() argument must be a string or a number, not 'dok_matrix'”。

如果我使用lil_matrix，错误信息是“ValueError: setting an array element with a sequence.”。

@hpaulj：是的，a[:]=s.A 做了我想要的，但它创建了一个新的临时数组，我想避免这种情况。它的另一个问题是s.A 现在仅 可以与稀疏数组一起使用，而不再与“普通”NumPy 数组一起使用，所以我会有更多的代码路径，我需要某种@987654361 @检查决定走哪条路。是的，我还想插入所有的零。

假设稀疏矩阵是一个“序列”，它在迭代时应该产生什么值？只是非零，还是全部（好像它很密集）？

【参考方案1】：

正如对我的问题的评论中提到的，序列接口不适用于稀疏矩阵，因为它们在使用单个数字索引时不会丢失维度。 无论如何，为了尝试它，我在纯 Python 中创建了一个非常有限的快速和肮脏的稀疏 array 类，当用单个数字索引时，它返回一个“行”类（它持有一个视图原始数据），它再次可以用单个数字索引以在该索引处产生实际值。使用我班级的实例 s，分配给 NumPy 数组 a 完全按照要求工作：

a[:] = s

我预计这会有些低效，但它真的、真的、真的、非常慢。分配一个 500.000 x 100 的稀疏数组需要几分钟！ 不过，好消息是在分配期间没有创建完整大小的临时数组。在分配期间内存使用量保持不变（当其中一个 CPU 达到最大值时）。

所以这基本上是原始问题的一种解决方案。

为了使分配更高效并且仍然不使用密集数组数据的临时副本，NumPy 必须在内部执行类似于

s.toarray(out=a)

据我所知，目前没有办法让 NumPy 做到这一点。

但是，有一种方法可以做一些非常相似的事情，方法是提供一个返回 NumPy 数组的__array__() 方法。顺便说一句，SciPy 稀疏矩阵已经有了这样的方法，只是名称不同：toarray()。所以我只是重命名了它：

scipy.sparse.dok_matrix.__array__ = scipy.sparse.dok_matrix.toarray
a[:] = s

这就像一个魅力（也适用于其他稀疏矩阵类）并且非常快！

根据我对这种情况的有限理解，这应该创建一个与a 大小相同的临时 NumPy 数组，该数组包含来自s 的所有值（以及许多零），然后分配给a . 但奇怪的是，即使我使用一个非常大的 a 几乎占据了我所有的可用 RAM，分配仍然很快发生，并且没有使用额外的 RAM。

所以我想这是我原来问题的另一个更好的解决方案。

这留下了另一个问题：为什么没有临时数组也能工作？

【参考方案2】：

如何使用nonzero 来识别哪些元素不为零？

x = np.ones((3,4))
s = sparse.csr_matrix((3,4))
s[0,0] = 2
s[1,2] = 3
I,J = s.nonzero()
x[:] = 0  # omit if just changing nonzero values
x[I,J] = s.data
x

nonzero 对密集数组和csr 数组的功能相同。其他格式我没试过。

---

对于 csr（和 coo）稀疏矩阵，值存储在 s.data 数组中。在这个例子中，它看起来像：

array([ 2.,  4.])

x 值位于data 缓冲区x.data 中。在这种情况下，它是 12 个连续的浮点数。

x.ravel()
# array([ 2.,  0.,  0.,  0.,  0.,  0.,  4.,  0.,  0.,  0.,  0.,  0.])

s 的这 2 个值无法在不复制的情况下映射到 x 的 12 个值。稀疏数据值通常不会与其密集等效值中的连续值块匹配。

您担心I 和J 数组的大小。如果稀疏矩阵是coo 格式，它的row 和col 属性可以以相同的方式使用：

sc=s.tocoo()
x[sc.row, sc.col]=sc.data

但是从一种稀疏格式转换为另一种需要复制数据。而且 coo 数组是不可下标的。

---

x = np.zeros((3,4))
x[:]=['123','321','0','1']

生产

array([[ 123.,  321.,    0.,    1.],
       [ 123.,  321.,    0.,    1.],
       [ 123.,  321.,    0.,    1.]])

它确实将float 应用于右侧的每个项目，然后“广播”它以适应x 的大小。

[] 转换为对-_setitem__ 的调用。

x.__setitem__((1,2),3)  # same as x[1,2]=3
x.__setitem__((None,2),'3') # sets the last row

它似乎在任何可迭代的每个项目上调用float（需要仔细检查）。但如果该值是其他对象，我们会收到类似于您原来的错误：

class Foo(): pass
x.__setitem__((1,2), Foo())
# TypeError: float() argument must be a string or a number, not 'Foo'

稀疏的coo 和dok 格式会产生类似的错误，而csr 和lil 会产生一个

ValueError: setting an array element with a sequence.

我还没弄清楚这里使用的是稀疏矩阵的哪个方法或属性。

看看np.broadcast。我认为这复制了这些作业中使用的迭代类型。

b = np.broadcast(x[:], [1,2,3,4])
list(b)

我们可以从一个包含 dtype 对象的数组开始来消除浮点转换的复杂性，它可以容纳任何东西：

xa=np.zeros((3,4), dtype=object)
xa[:]=s

但现在s 出现在xa 的每个元素中。它没有将s 的值分配给xa。

我猜当s 不是np.array 时，numpy 在做作业时首先将其包装起来，例如：

x[:] = np.array(s)

当s 是一个标量或列表时，这会生成一个可以广播以适应x 的数组。但是当它是一个对象时（稀疏数组不是 numpy 数组），这个包装只是一个 dtype=object 的 0d 数组。您需要通过一个函数将s 传递给将其转换为可以广播的可迭代对象。最明显的是toarray()。

【讨论】：

这将避免将整个矩阵复制到一个新数组中，但仍需要创建新的中间数组（I 和J）。

另外，稀疏情况需要额外的代码。我真的很想保持a[:] = s 部分原样，并通过以某种方式操纵 s 事先解决问题。

“操作 s”是指填充零个元素，以便无需特殊索引即可将它们复制到a？

如何操作s 是我不知道的部分......我只想对s 做一些事情，这使得下面的a[:] = s 工作。更具体地说（我在最初的问题中忽略了这部分），我想切入s，它应该仍然可以工作：a[:] = s[start:stop]。所以基本上我需要一些存储稀疏数据的东西，并允许在不复制数据的情况下将结果切片和分配给 NumPy 数组。

我添加了一些关于稀疏数组中数据布局的注释。