在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

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]]) 

干杯， 达维德

Answer 2

我觉得这个问题并不完全清楚。我想这个稀疏矩阵给你从原始矩阵的每个系数中取出的样本数量？

查看链接中的代码，可能会加快速度。在转置的矩阵上操作并重写链接的代码以在列上而不是在行上操作。因为这样可以让你的处理器缓存它采样得更好的值，也就是说内存中的跳转较少。

其余的就像我会这样做，使用numpy（并不一定意味着这是最有效的方式）。

如果你需要它更快，你可以尝试在C++中编写函数，并用scipy.weave将它包含到你的python中。在C++中，我会去每一行，并建立一个> 0的位置查找表，在与查找表中的项目数相等的范围内生成整数。我会累积查找表中每个位置的绘制频率，然后将这些数字放回到数组中的正确位置。该代码应该只是简单的几行。