按特定顺序将空列插入 scipy 稀疏矩阵

Question

我有一个包含 M 行和 N 列的稀疏矩阵，我想将 K 个额外的 NULL 列连接到该矩阵，这样我的对象现在将具有 M 行和 (N+K) 列。棘手的部分是我还有一个长度为 N 的索引列表，范围从 0 到 N+K，它指示每一列在新矩阵中的位置。

例如，如果 N = 2，K = 1 并且索引列表是 [2, 0]，这意味着我想将 MxN 矩阵的最后一列作为第一列，引入一个空列，然后将我的第一列作为最后一列。

我正在尝试使用以下代码 - 当我已经有 x 但无法在此处上传时。

import numpy as np
from scipy import sparse
M = 5000
N = 10
pad_factor = 1.2
size = int(pad_factor * N)
x = sparse.random(m = M, n = N, density = 0.1, dtype = 'float64')
indeces = np.random.choice(range(size), size=N, replace=False)
null_mat = sparse.lil_matrix((M, size))
null_mat[:, indeces] = x

问题是对于 N = 1,500,000、P = 5,000 和 K = 200，此代码将无法缩放并且会给我一个内存错误。确切的错误是： "return np.zeros(self.shape, dtype = self.dtype, order=order) MemoryError"。

我只想添加一些空列，所以我猜我的切片想法效率低下，尤其是当我的真实数据中有 K << N 时。在某种程度上，我们可以将其视为合并排序问题——我有一个非空数据集和一个空数据集，我想将它们连接起来，但要按特定顺序连接。关于如何让它发挥作用有什么想法吗？

谢谢！

Answer 1

正如我在评论中推断的那样，内存错误是在

null_mat[:, indeces] = x

行因为lil__setitem__方法，做了一个x.toarray()，也就是先把x转成密集数组。将稀疏矩阵直接映射到索引 lil 可能更 space 有效，但需要更多的编码工作。并且 lil 针对迭代分配进行了优化，而不是这种大规模矩阵映射。

sparse.hstack 使用 sparse.bmat 连接稀疏矩阵。这会将所有输入转换为 coo，然后将它们的属性组合成一个新集合，从中构建新矩阵。

直接构建coo矩阵

经过相当多的尝试，我发现以下简单操作有效：

In [479]: z1=sparse.coo_matrix((x.data, (x.row, indeces[x.col])),shape=(M,size))

In [480]: z1
Out[480]: 
<5000x12 sparse matrix of type '<class 'numpy.float64'>'
    with 5000 stored elements in COOrdinate format>

将此与 x 和 null_mat 进行比较：

In [481]: x
Out[481]: 
<5000x10 sparse matrix of type '<class 'numpy.float64'>'
    with 5000 stored elements in COOrdinate format>
In [482]: null_mat
Out[482]: 
<5000x12 sparse matrix of type '<class 'numpy.float64'>'
    with 5000 stored elements in LInked List format>

测试稀疏矩阵的相等性可能很棘手。 coo 值可以以任何顺序出现，例如 x 由 sparse.random.

生成

但是 csr 格式对行进行排序，因此 indptr 属性的这种比较是一个很好的相等性测试：

In [483]: np.allclose(null_mat.tocsr().indptr, z1.tocsr().indptr)
Out[483]: True

时间测试：

In [477]: timeit z1=sparse.coo_matrix((x.data, (x.row, indeces[x.col])),shape=(M,size))
108 µs ± 1.24 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [478]: 
In [478]: timeit null_mat[:, indeces] = x
3.05 ms ± 4.55 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

矩阵乘法方法

csr 列表的格式索引是通过矩阵乘法完成的。它构造一个 extractor 矩阵，并应用它。矩阵乘法是csr_matrix的强项。

我们可以用同样的方式进行重新排序：

In [489]: I = sparse.csr_matrix((np.ones(10),(np.arange(10),indeces)), shape=(10,12))
In [490]: I
Out[490]: 
<10x12 sparse matrix of type '<class 'numpy.float64'>'
    with 10 stored elements in Compressed Sparse Row format>

In [496]: w1=x*I

比较这些矩阵的密集等价物：

In [497]: np.allclose(null_mat.A, z1.A)
Out[497]: True
In [498]: np.allclose(null_mat.A, w1.A)
Out[498]: True


In [499]: %%timeit
     ...: I = sparse.csr_matrix((np.ones(10),(np.arange(10),indeces)),shape=(10,
     ...: 12))
     ...: w1=x*I
1.11 ms ± 5.65 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

这比 lil 索引方法要好，但仍然比直接 coo 矩阵构造慢得多。公平地说，我们应该从 coo 样式输入构建一个 csr 矩阵。该转换需要一些时间：

In [502]: timeit z2=sparse.csr_matrix((x.data, (x.row, indeces[x.col])),shape=(M
     ...: ,size))
639 µs ± 604 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

错误回溯

MemoryError 回溯应该表明错误发生在这个索引赋值中，相关的方法调用是：

Signature: null_mat.__setitem__(index, x)
Source:   
    def __setitem__(self, index, x):
       ....
       if isspmatrix(x):
           x = x.toarray()
       ...

Signature: x.toarray(order=None, out=None)
Source:   
    def toarray(self, order=None, out=None):
        """See the docstring for `spmatrix.toarray`."""
        B = self._process_toarray_args(order, out)
Signature: x._process_toarray_args(order, out)
Source:   
    def _process_toarray_args(self, order, out):
        ...
        return np.zeros(self.shape, dtype=self.dtype, order=order)

我通过在 scipy github 上搜索 np.zeros 调用的代码找到了这个。

按特定顺序将空列插入 scipy 稀疏矩阵

Inserting null columns into a scipy sparse matrix in a specific order

python

scipy

sparse-matrix

直接构建coo矩阵

矩阵乘法方法

错误回溯