随机矩阵所有行的快速随机加权选择

Question

numpy.random.choice 允许从向量中进行加权选择，即

arr = numpy.array([1, 2, 3])
weights = numpy.array([0.2, 0.5, 0.3])
choice = numpy.random.choice(arr, p=weights)

选择概率为0.2的1，概率为0.5的2和概率为0.3的3。

如果我们想对每个行都是概率向量的2D数组（矩阵）以向量化的方式快速进行操作，该怎么办？也就是说，我们想要一个来自随机矩阵的选择向量吗？这是超级慢的方式：

import numpy as np

m = 10
n = 100 # Or some very large number

items = np.arange(m)
prob_weights = np.random.rand(m, n)
prob_matrix = prob_weights / prob_weights.sum(axis=0, keepdims=True)

choices = np.zeros((n,))
# This is slow, because of the loop in Python
for i in range(n):
    choices[i] = np.random.choice(items, p=prob_matrix[:,i])

print(choices)：

array([ 4.,  7.,  8.,  1.,  0.,  4.,  3.,  7.,  1.,  5.,  7.,  5.,  3.,
        1.,  9.,  1.,  1.,  5.,  9.,  8.,  2.,  3.,  2.,  6.,  4.,  3.,
        8.,  4.,  1.,  1.,  4.,  0.,  1.,  8.,  5.,  3.,  9.,  9.,  6.,
        5.,  4.,  8.,  4.,  2.,  4.,  0.,  3.,  1.,  2.,  5.,  9.,  3.,
        9.,  9.,  7.,  9.,  3.,  9.,  4.,  8.,  8.,  7.,  6.,  4.,  6.,
        7.,  9.,  5.,  0.,  6.,  1.,  3.,  3.,  2.,  4.,  7.,  0.,  6.,
        3.,  5.,  8.,  0.,  8.,  3.,  4.,  5.,  2.,  2.,  1.,  1.,  9.,
        9.,  4.,  3.,  3.,  2.,  8.,  0.,  6.,  1.])

这篇文章表明，cumsum并且bisect可能是一种潜在的方法，而且很快。但是虽然numpy.cumsum(arr, axis=1)可以沿numpy数组的一个轴执行此操作，但是该bisect.bisect函数一次只能在单个数组上运行。同样，也numpy.searchsorted仅适用于一维阵列。

是否有一种仅使用矢量化操作来完成此操作的快速方法？

Answer 1

这是一个非常快速的完全矢量化版本：

def vectorized(prob_matrix, items):
    s = prob_matrix.cumsum(axis=0)
    r = np.random.rand(prob_matrix.shape[1])
    k = (s < r).sum(axis=0)
    return items[k]

从理论上讲 ，searchsorted是用于在累积求和的概率中查找随机值的正确函数，但m由于相对较小，因此k = (s < r).sum(axis=0)最终速度要快得多。它的时间复杂度为O（m），而searchsorted方法为O（log（m）），但这仅对大得多才有意义m。
另外
，cumsum也是O（m），所以vectorized和@perimosocordiaeimproved都是O（m）。（m实际上，如果您的计算机更大，则必须运行一些测试以查看m此方法变慢之前可以达到的大小。）

这是我使用m = 10和n = 10000（使用函数original和improved@perimosocordiae的答案）得到的时间：

In [115]: %timeit original(prob_matrix, items)
1 loops, best of 3: 270 ms per loop

In [116]: %timeit improved(prob_matrix, items)
10 loops, best of 3: 24.9 ms per loop

In [117]: %timeit vectorized(prob_matrix, items)
1000 loops, best of 3: 1 ms per loop

定义功能的完整脚本为：

import numpy as np


def improved(prob_matrix, items):
    # transpose here for better data locality later
    cdf = np.cumsum(prob_matrix.T, axis=1)
    # random numbers are expensive, so we'll get all of them at once
    ridx = np.random.random(size=n)
    # the one loop we can't avoid, made as simple as possible
    idx = np.zeros(n, dtype=int)
    for i, r in enumerate(ridx):
        idx[i] = np.searchsorted(cdf[i], r)
    # fancy indexing all at once is faster than indexing in a loop
    return items[idx]


def original(prob_matrix, items):
    choices = np.zeros((n,))
    # This is slow, because of the loop in Python
    for i in range(n):
        choices[i] = np.random.choice(items, p=prob_matrix[:,i])
    return choices


def vectorized(prob_matrix, items):
    s = prob_matrix.cumsum(axis=0)
    r = np.random.rand(prob_matrix.shape[1])
    k = (s < r).sum(axis=0)
    return items[k]


m = 10
n = 10000 # Or some very large number

items = np.arange(m)
prob_weights = np.random.rand(m, n)
prob_matrix = prob_weights / prob_weights.sum(axis=0, keepdims=True)