如何加快 Numpy 中转换矩阵的创建速度？答案

【问题标题】：How can I speed up transition matrix creation in Numpy?如何加快 Numpy 中转换矩阵的创建速度？
【发布时间】：2012-11-05 03:34:23
【问题描述】：

以下是我所知道的对马尔可夫链中的转换进行计数并使用它来填充转换矩阵的最基本方法：

def increment_counts_in_matrix_from_chain(markov_chain, transition_counts_matrix):
    for i in xrange(1, len(markov_chain)):
        old_state = markov_chain[i - 1]
        new_state = markov_chain[i]
        transition_counts_matrix[old_state, new_state] += 1

我尝试了 3 种不同的方式来加速它：

1) 使用基于此 Matlab 代码的稀疏矩阵单线：

transition_matrix = full(sparse(markov_chain(1:end-1), markov_chain(2:end), 1))

在 Numpy/SciPy 中，如下所示：

def get_sparse_counts_matrix(markov_chain, number_of_states):
    return coo_matrix(([1]*(len(markov_chain) - 1), (markov_chain[0:-1], markov_chain[1:])), shape=(number_of_states, number_of_states))

我还尝试了一些 Python 调整，比如使用 zip()：

for old_state, new_state in zip(markov_chain[0:-1], markov_chain[1:]):
    transition_counts_matrix[old_state, new_state] += 1

还有队列：

old_and_new_states_holder = Queue(maxsize=2)
old_and_new_states_holder.put(markov_chain[0])
for new_state in markov_chain[1:]:
    old_and_new_states_holder.put(new_state)
    old_state = old_and_new_states_holder.get()
    transition_counts_matrix[old_state, new_state] += 1

但是这 3 种方法都没有加快速度。事实上，除了 zip() 解决方案之外的所有解决方案都至少比我原来的解决方案慢 10 倍。

还有其他值得研究的解决方案吗？

从大量链构建转移矩阵的改进解决方案
上述问题的最佳答案是帝斯曼。但是，对于任何想要根据数百万个马尔可夫链的列表填充转换矩阵的人来说，最快的方法是：

def fast_increment_transition_counts_from_chain(markov_chain, transition_counts_matrix):
    flat_coords = numpy.ravel_multi_index((markov_chain[:-1], markov_chain[1:]), transition_counts_matrix.shape)
    transition_counts_matrix.flat += numpy.bincount(flat_coords, minlength=transition_counts_matrix.size)

def get_fake_transitions(markov_chains):
    fake_transitions = []
    for i in xrange(1,len(markov_chains)):
        old_chain = markov_chains[i - 1]
        new_chain = markov_chains[i]
        end_of_old = old_chain[-1]
        beginning_of_new = new_chain[0]
        fake_transitions.append((end_of_old, beginning_of_new))
    return fake_transitions

def decrement_fake_transitions(fake_transitions, counts_matrix):
    for old_state, new_state in fake_transitions:
        counts_matrix[old_state, new_state] -= 1

def fast_get_transition_counts_matrix(markov_chains, number_of_states):
    """50% faster than original, but must store 2 additional slice copies of all markov chains in memory at once.
    You might need to break up the chains into manageable chunks that don't exceed your memory.
    """
    transition_counts_matrix = numpy.zeros([number_of_states, number_of_states])
    fake_transitions = get_fake_transitions(markov_chains)
    markov_chains = list(itertools.chain(*markov_chains))
    fast_increment_transition_counts_from_chain(markov_chains, transition_counts_matrix)
    decrement_fake_transitions(fake_transitions, transition_counts_matrix)
    return transition_counts_matrix

【问题讨论】：

标签： python numpy scipy

【解决方案1】：

只是为了好玩，因为我一直想尝试一下，所以我申请了Numba 来解决您的问题。在代码中，只需要添加一个装饰器（虽然我已经直接调用，所以我可以测试 numba 在这里提供的 jit 变体）：

import numpy as np
import numba

def increment_counts_in_matrix_from_chain(markov_chain, transition_counts_matrix):
    for i in xrange(1, len(markov_chain)):
        old_state = markov_chain[i - 1]
        new_state = markov_chain[i]
        transition_counts_matrix[old_state, new_state] += 1

autojit_func = numba.autojit()(increment_counts_in_matrix_from_chain)
jit_func = numba.jit(argtypes=[numba.int64[:,::1],numba.double[:,::1]])(increment_counts_in_matrix_from_chain)

t = np.random.randint(0,50, 500)
m1 = np.zeros((50,50))
m2 = np.zeros((50,50))
m3 = np.zeros((50,50))

然后是时间安排：

In [10]: %timeit increment_counts_in_matrix_from_chain(t,m1)
100 loops, best of 3: 2.38 ms per loop

In [11]: %timeit autojit_func(t,m2)                         

10000 loops, best of 3: 67.5 us per loop

In [12]: %timeit jit_func(t,m3)
100000 loops, best of 3: 4.93 us per loop

autojit 方法根据运行时输入进行一些猜测，jit 函数具有指定的类型。您必须小心一点，因为在这些早期阶段的 numba 不会传达jit 存在错误，如果您为输入传递了错误的类型。它只会吐出一个不正确的答案。

尽管如此，在不更改任何代码的情况下获得 35 倍和 485 倍的加速，并且只添加对 numba 的调用（也可以称为装饰器）在我的书中是相当令人印象深刻的。使用 cython 可能会得到类似的结果，但它需要更多样板文件并编写 setup.py 文件。

我也喜欢这个解决方案，因为代码仍然可读，并且您可以按照您最初考虑实现算法的方式编写它。

【讨论】：

整洁！启动成本是多少？
@DSM 不确定这是否是计时的最佳方式，但%timeit autojit_func = numba.autojit()(increment_counts_in_matrix_from_chain); autojit_func(t,m2) 给了 81 我们。当我为普通的jit 做类似的事情时，我会收到一堆垃圾收集警告，我认为这些警告会搞砸时间。

【解决方案2】：

利用np.bincount 这样的事情怎么样？不是超级强大，但功能强大。 [感谢@Warren Weckesser 的设置。]

import numpy as np
from collections import Counter

def increment_counts_in_matrix_from_chain(markov_chain, transition_counts_matrix):
    for i in xrange(1, len(markov_chain)):
        old_state = markov_chain[i - 1]
        new_state = markov_chain[i]
        transition_counts_matrix[old_state, new_state] += 1

def using_counter(chain, counts_matrix):
    counts = Counter(zip(chain[:-1], chain[1:]))
    from_, to = zip(*counts.keys())
    counts_matrix[from_, to] = counts.values()

def using_bincount(chain, counts_matrix):
    flat_coords = np.ravel_multi_index((chain[:-1], chain[1:]), counts_matrix.shape)
    counts_matrix.flat = np.bincount(flat_coords, minlength=counts_matrix.size)

def using_bincount_reshape(chain, counts_matrix):
    flat_coords = np.ravel_multi_index((chain[:-1], chain[1:]), counts_matrix.shape)
    return np.bincount(flat_coords, minlength=counts_matrix.size).reshape(counts_matrix.shape)

给出：

In [373]: t = np.random.randint(0,50, 500)
In [374]: m1 = np.zeros((50,50))
In [375]: m2 = m1.copy()
In [376]: m3 = m1.copy()

In [377]: timeit increment_counts_in_matrix_from_chain(t, m1)
100 loops, best of 3: 2.79 ms per loop

In [378]: timeit using_counter(t, m2)
1000 loops, best of 3: 924 us per loop

In [379]: timeit using_bincount(t, m3)
10000 loops, best of 3: 57.1 us per loop

[编辑]

避免flat（以不在原地工作为代价）可以为小型矩阵节省一些时间：

In [80]: timeit using_bincount_reshape(t, m3)
10000 loops, best of 3: 22.3 us per loop

【讨论】：

我会接受这个答案，但我想跟进一个额外的问题。当我反复使用 bincount 来填充基于数千个马尔可夫链的转换计数矩阵时，我的原始代码仍然更快。我认为这是因为 counts_matrix.flat += numpy.bincount(flat_coords, minlength=counts_matrix.size) 在更新 counts_matrix 时比我的原始代码慢。对此有什么想法？
对此更新：我找到的基于大量马尔可夫链填充转换矩阵的最快解决方案是将链一个接一个地合并在一起，使用 bincounts，然后获取假转换（来自一个链的末端到下一个链的开头），然后减少每个假转换的计数。该解决方案比我原来的解决方案快 25%。
@some-guy：随意采用您为您的用例找到的最佳解决方案，将其作为答案发布并接受。

【解决方案3】：

这是一种更快的方法。这个想法是计算每个转换的出现次数，并在矩阵的矢量化更新中使用这些计数。（我假设相同的转换可以在 markov_chain 中多次出现。）collections 库中的 Counter 类用于计算每个转换的出现次数。

from collections import Counter

def update_matrix(chain, counts_matrix):
    counts = Counter(zip(chain[:-1], chain[1:]))
    from_, to = zip(*counts.keys())
    counts_matrix[from_, to] += counts.values()

时序示例，在 ipython 中：

In [64]: t = np.random.randint(0,50, 500)

In [65]: m1 = zeros((50,50))

In [66]: m2 = zeros((50,50))

In [67]: %timeit increment_counts_in_matrix_from_chain(t, m1)
1000 loops, best of 3: 895 us per loop

In [68]: %timeit update_matrix(t, m2)
1000 loops, best of 3: 504 us per loop

它更快，但不是快几个数量级。为了真正加快速度，您可以考虑在 Cython 中实现它。

【讨论】：

【解决方案4】：

好的，很少有可以篡改的想法，稍有改进（以人类不理解为代价）

让我们从长度为 3000 的 0 到 9 之间整数的随机向量开始：

L = 3000
N = 10
states = array(randint(N),size=L)
transitions = np.zeros((N,N))

在我的机器上，您的方法的 timeit 性能为 11.4 毫秒。

首先要稍微改进一下，避免两次读取数据，将其存储在一个临时变量中：

old = states[0]
for i in range(1,len(states)):
    new = states[i]
    transitions[new,old]+=1
    old=new

这会给您带来约 10% 的改进，并将时间缩短到 10.9 毫秒。

一种更内卷的方法使用跨步：

def rolling(a, window):
    shape = (a.size - window + 1, window)
    strides = (a.itemsize, a.itemsize)
    return np.lib.stride_tricks.as_strided(a, shape=shape, strides=strides)

state_2 = rolling(states, 2)
for i in range(len(state_2)):
    l,m = state_2[i,0],state_2[i,1]
    transitions[m,l]+=1

strides 允许您读取数组的连续数字，从而诱使数组认为行以不同的方式开始（好吧，它没有很好地描述，但是如果您花一些时间阅读有关 strides 的信息，您会明白的) 这种方法会损失性能，达到 12.2 毫秒，但它是更多欺骗系统的通道。将转换矩阵和跨步数组都展平为一维数组，可以进一步提高性能：

transitions = np.zeros(N*N)
state_2 = rolling(states, 2)
state_flat = np.sum(state_2 * array([1,10]),axis=1)
for i in state_flat:
    transitions[i]+=1
transitions.reshape((N,N))

这下降到 7.75 毫秒。这不是一个数量级，但无论如何都要好 30% :)

【讨论】：