沿轴0重复scipy csr稀疏矩阵

Question

我想重复一个scipy的csr稀疏矩阵的行，但当我尝试使用numpy的repeat方法时，它只是把稀疏矩阵当成一个对象来处理，结果只是在一个ndarray中重复这个对象。我查阅了文档，但没找到任何可以重复scipy csr稀疏矩阵行的工具。

于是我写了以下代码，直接操作内部数据，似乎能正常工作

def csr_repeat(csr, repeats):
    if isinstance(repeats, int):
        repeats = np.repeat(repeats, csr.shape[0])
    repeats = np.asarray(repeats)
    rnnz = np.diff(csr.indptr)
    ndata = rnnz.dot(repeats)
    if ndata == 0:
        return sparse.csr_matrix((np.sum(repeats), csr.shape[1]),
                                 dtype=csr.dtype)
    indmap = np.ones(ndata, dtype=np.int)
    indmap[0] = 0
    rnnz_ = np.repeat(rnnz, repeats)
    indptr_ = rnnz_.cumsum()
    mask = indptr_ < ndata
    indmap -= np.int_(np.bincount(indptr_[mask],
                                  weights=rnnz_[mask],
                                  minlength=ndata))
    jumps = (rnnz * repeats).cumsum()
    mask = jumps < ndata
    indmap += np.int_(np.bincount(jumps[mask],
                                  weights=rnnz[mask],
                                  minlength=ndata))
    indmap = indmap.cumsum()
    return sparse.csr_matrix((csr.data[indmap],
                              csr.indices[indmap],
                              np.r_[0, indptr_]),
                             shape=(np.sum(repeats), csr.shape[1]))

而且效率也还不错，但我不太想直接修改这个类。有没有更好的方法来做到这一点呢？

编辑

当我重新审视这个问题时，我在想我当初为什么要发这个帖子。几乎所有我能想到的对重复矩阵的操作，用原始矩阵来做会更简单，然后再进行重复处理。我的想法是，重复之后再处理这个问题，通常会比任何可能的答案都要好。

Answer 1

重复稀疏矩阵的一个非常有效的方法就是OP提到的那种。我修改了 indptr，这样就不会输出全是0的行了。

## original sparse matrix
indptr = np.array([0, 2, 3, 6])
indices = np.array([0, 2, 2, 0, 1, 2])
data = np.array([1, 2, 3, 4, 5, 6])
x = scipy.sparse.csr_matrix((data, indices, indptr), shape=(3, 3))
x.toarray()

array([[1, 0, 2],
       [0, 0, 3],
       [4, 5, 6]])

要重复这个，你需要重复数据和索引，同时还得调整一下 indptr。这虽然不是最优雅的方法，但确实有效。

## repeated sparse matrix
repeat = 5 
new_indptr = indptr
for r in range(1,repeat):
    new_indptr = np.concatenate((new_indptr, new_indptr[-1]+indptr[1:]))
x = scipy.sparse.csr_matrix((np.tile(data,repeat), np.tile(indices,repeat), new_indptr))
x.toarray()

array([[1, 0, 2],
       [0, 0, 3],
       [4, 5, 6],
       [1, 0, 2],
       [0, 0, 3],
       [4, 5, 6],
       [1, 0, 2],
       [0, 0, 3],
       [4, 5, 6],
       [1, 0, 2],
       [0, 0, 3],
       [4, 5, 6],
       [1, 0, 2],
       [0, 0, 3],
       [4, 5, 6]])

Answer 2

在有人给出一个非常聪明的回答后，我又回过头来看我的原始问题，想看看有没有更好的方法。于是我想出了另一种方法，这种方法有一些优缺点。我们可以让scipy重复使用那些重复行的数据，而不是像接受的答案那样重复所有数据，这样就能创建一个类似于原始稀疏数组的视图（就像你可以用broadcast_to那样）。这可以通过简单地调整indptr字段来实现。

repeated = sparse.csr_matrix((orig.data, orig.indices, np.tile(orig.indptr, repeat_num)))

这个技巧会将向量重复repeat_num次，同时只修改indptr。缺点是，由于csr矩阵编码数据的方式，它并不会创建一个尺寸为repeat_num x n的矩阵，而是创建一个尺寸为(2 * repeat_num - 1) x n的矩阵，其中每一行的奇数行都是0。虽然这不算太麻烦，因为每一行都是0，任何操作都会很快，而且之后用类似[::2]的方式切片也很简单，但这并不是最理想的情况。

我认为被标记的答案可能仍然是“最佳”方法。

Answer 3

其实，np.repeat不工作的原因并不奇怪。它把这个操作交给了硬编码的a.repeat方法，如果这个方法不行，它会先把a转成一个数组（如果需要的话，会变成对象）。

在稀疏代码开发的线性代数领域，大部分的组装工作都是在创建稀疏矩阵之前，先在row、col和data数组上完成的。这里的重点是高效的数学运算，而不是添加、删除或索引行和元素。

我没有仔细研究你的代码，但csr格式的矩阵需要这么多工作也不让我感到意外。

我为lil格式写了一个类似的函数（是从lil.copy开始的）：

def lil_repeat(S, repeat):
    # row repeat for lil sparse matrix
    # test for lil type and/or convert
    shape=list(S.shape)
    if isinstance(repeat, int):
        shape[0]=shape[0]*repeat
    else:
        shape[0]=sum(repeat)
    shape = tuple(shape)
    new = sparse.lil_matrix(shape, dtype=S.dtype)
    new.data = S.data.repeat(repeat) # flat repeat
    new.rows = S.rows.repeat(repeat)
    return new

不过也可以通过索引来重复。lil和csr都支持类似于普通numpy数组的索引（至少在较新的版本中是这样）。所以：

S = sparse.lil_matrix([[0,1,2],[0,0,0],[1,0,0]])
print S.A.repeat([1,2,3], axis=0)
print S.A[(0,1,1,2,2,2),:]
print lil_repeat(S,[1,2,3]).A
print S[(0,1,1,2,2,2),:].A

会得到相同的结果。

最棒的是？

print S[np.arange(3).repeat([1,2,3]),:].A

Answer 4

在编程中，有时候我们会遇到一些问题，像是代码运行不正常或者出现错误。这些问题可能是因为我们写的代码有些地方不对，或者是我们没有正确理解某些概念。

比如说，有些人可能会在使用某个功能时，发现它并没有按照预期工作。这时候，我们就需要仔细检查代码，看看是不是哪里出错了。通常，错误信息会给我们一些提示，帮助我们找到问题所在。

另外，了解一些基本的编程概念也很重要。比如变量、循环、条件判断等等，这些都是编程的基础。如果我们能掌握这些基础知识，就能更好地理解代码是如何运行的，也能更容易地解决问题。

总之，编程就像解谜一样，有时候需要耐心和细心去找出答案。遇到问题时，不要着急，慢慢分析，通常就能找到解决办法。

from scipy.sparse import csr_matrix
repeated_row_matrix = csr_matrix(np.ones([repeat_number,1])) * sparse_row