快速稀疏Numpy/Python中的矩阵

Question

我需要快速对一个矩阵进行稀疏化处理。

稀疏化就是把丰度矩阵转换成均匀的采样深度。

在这个例子中，每一行代表一个样本，而采样深度就是这一行所有数值的总和。我想随机抽取（可以重复抽取）这个矩阵中的样本，抽取的次数是 min(rowsums(matrix))。

假设我有一个矩阵：

>>> m = [ [0, 9, 0],
...       [0, 3, 3],
...       [0, 4, 4] ]

稀疏化的过程是逐行进行的，每一行随机抽取 min(rowsums(matrix)) 次（在这个例子中是6次）。

>>> rf = rarefaction(m)
>>> rf
    [ [0, 6, 0],  # sum = 6
      [0, 3, 3],  # sum = 6
      [0, 3, 3] ] # sum = 6

结果是随机的，但每一行的总和始终是一样的。

>>> rf = rarefaction(m)
>>> rf
    [ [0, 6, 0],   # sum = 6
      [0, 2, 4],   # sum = 6
      [0, 4, 2], ] # sum = 6

PyCogent 有一个函数可以逐行处理这个过程，但在处理大矩阵时速度非常慢。

我觉得Numpy里可能有一个函数可以做到这一点，但我不太确定它叫什么。

Answer 1

我觉得这个问题不是特别清楚。我猜这个稀疏矩阵是用来告诉你从原始矩阵的每个系数中取了多少个样本吧？

看你链接里的代码，可能有提升速度的空间。可以试着对调转的矩阵进行操作，并把你的代码改成处理列而不是行。这样做可以让你的处理器更好地缓存取样的值，也就是说，内存中的跳转会更少。

其他部分我也会这样做，使用numpy（这并不一定是最有效的方法）。

如果你需要更快的速度，可以尝试用C++编写这个函数，然后通过scipy.weave把它包含进你的python代码里。在C++中，我会遍历每一行，建立一个大于0的位置信息表，生成min(rowsums(matrix))个整数，这些整数的范围等于位置信息表中的项目数量。我会统计每个位置在位置信息表中被抽取的次数，然后把这些数字放回数组的正确位置。这个代码应该实际上只需要几行。

Answer 2

import numpy as np
from numpy.random import RandomState

def rarefaction(M, seed=0):
    prng = RandomState(seed) # reproducible results
    noccur = np.sum(M, axis=1) # number of occurrences for each sample
    nvar = M.shape[1] # number of variables
    depth = np.min(noccur) # sampling depth

    Mrarefied = np.empty_like(M)
    for i in range(M.shape[0]): # for each sample
        p = M[i] / float(noccur[i]) # relative frequency / probability
        choice = prng.choice(nvar, depth, p=p)
        Mrarefied[i] = np.bincount(choice, minlength=nvar)

    return Mrarefied

例子：

>>> M = np.array([[0, 9, 0], [0, 3, 3], [0, 4, 4]])
>>> M
array([[0, 9, 0],
       [0, 3, 3],
       [0, 4, 4]])
>>> rarefaction(M)
array([[0, 6, 0],
       [0, 2, 4],
       [0, 3, 3]])
>>> rarefaction(M, seed=1)
array([[0, 6, 0],
       [0, 4, 2],
       [0, 3, 3]])
>>> rarefaction(M, seed=2)
array([[0, 6, 0],
       [0, 3, 3],
       [0, 3, 3]])

谢谢，
达维德