如何有效地缓存 scipy gammaln 并将其应用于 numpy 数组？

Question

我正在寻找一种方法来缓存 gammaln 函数的结果，将其应用于 numpy 数组，然后对结果求和。目前我能够使用简单的字典方法缓存 gammln 结果，但似乎我无法有效地迭代数组。下面是问题的说明性代码片段。

import numpy as np
from scipy.special import gammaln

gammaln_cache_dic = {}
def gammaln_cache(x):
    if x not in gammaln_cache_dic:
        val = gammaln(x)
        gammaln_cache_dic[x] = val
    else:
        val = gammaln_cache_dic[x]
    return val

def cache_gammaln_sum(M):
    res_sum = 0.0
    for x in np.fromiter( [ x for x in M ], int ):
        res_sum += gammaln_cache(x)
    return res_sum

np_array = np.random.randint(5, size=(500))*1000+3
start = time.time()
for i in np_array:
    gammaln(i)
end = time.time()
print("NO cache gammaln %f secs" % (end - start))

start = time.time()
for i in np_array:
    gammaln_cache(i)
end = time.time()
print("cache gammaln %f secs" % (end - start))

start = time.time()
gammaln(np_array).sum()
end = time.time()
print("NO cache gammaln sum %f secs" % (end - start))

start = time.time()
cache_gammaln_sum(np_array)
end = time.time()
print("cache gammaln sum %f secs" % (end - start))

Answer 1

你可以使用np.unique

uniq, counts = np.unique(data, return_counts=True)
result = np.dot(gammaln(uniq), counts)

uniq 是数据中排序后的唯一值，return_counts 关键字指示函数也返回多重性。

第二行对每个不同的元素计算一次gammaln，乘以出现次数和总和。

您也可以保留gammaln(uniq) 并将其用作查找表。

look_up = gammaln(uniq)
indices = np.searchsorted(uniq, new_samples)
hit = uniq[indices] == new_samples
result = look_up[indices[hit]]
miss = new_samples[~hit]

Answer 2

gammaln，就像scipy.special 中的大多数函数一样，是一个 ufunc，即给定一个 numpy 数组，它按元素操作并返回一个 numpy 数组：

In [1]: import numpy as np

In [2]: from scipy.special import gammaln

In [3]: x = np.arange(1, 8)

In [4]: gammaln(x)
Out[4]: 
array([ 0.        ,  0.        ,  0.69314718,  1.79175947,  3.17805383,
        4.78749174,  6.57925121])

In [5]: gammaln(x).sum()
Out[5]: 17.029703434928095

实际上，如果你给它一个数字列表，它会在内部将它转换为一个数组，所以

In [6]: gammaln(list(range(1, 8)))
Out[6]: 
array([ 0.        ,  0.        ,  0.69314718,  1.79175947,  3.17805383,
        4.78749174,  6.57925121])

另请注意，ufunc 迭代发生在编译后的代码中。因此，您不太可能通过任何类型的缓存和 python 迭代等在性能方面做得更好。

Answer 3

如果编写函数以利用 numpy 编译代码（或本身已编译），则很难通过缓存提高其性能。

In [1]: from scipy.special import gammaln
In [2]: np_array = np.random.randint(5, size=(500))*1000+3
In [3]: np_array.shape
Out[3]: (500,)
In [4]: timeit gammaln(np_array).sum()
....
10000 loops, best of 3: 34.9 µs per loop

评估每个元素会慢很多

In [5]: timeit sum([gammaln(i) for i in np_array])
1000 loops, best of 3: 1.8 ms per loop

简单地遍历数组会更慢。

In [6]: timeit sum([i for i in np_array])
10000 loops, best of 3: 127 µs per loop

在这种情况下，值是一小组整数

In [10]: np.unique(np_array)
Out[10]: array([   3, 1003, 2003, 3003, 4003])

但识别这些唯一值所花费的时间几乎与计算所有 500 个函数的时间一样长。

In [11]: timeit np.unique(np_array)
...
In [12]: timeit gammaln(np.unique(np_array))

事实上，如果我还要求可以重建数组的索引，它会更慢：

In [14]: timeit np.unique(np_array, return_inverse=True)    
10000 loops, best of 3: 54 µs per loop

有用于 Python 的缓存或记忆配方，包括 @functools.lru_cache 装饰器。但是我还没有看到太多numpy 数组。如果您的问题适合缓存，我怀疑最好在较低级别的代码中执行此操作，例如使用 Cython，并使用现有的 C 或 C++ 缓存库。

具有唯一性的缓存

我们可以使用 unique 获得类似缓存的行为 - 使用它来获取唯一值和索引，以便我们重新创建原始数组：

In [5]: unq, idx=np.unique(np_array, return_inverse=True)
In [6]: res1= gammaln(np_array)
In [7]: res2= gammaln(unq)
In [8]: res2= res2[idx]
In [9]: np.allclose(res1, res2)
Out[9]: True

但是比较一下时代：

In [10]: timeit res1= gammaln(np_array)
10000 loops, best of 3: 29.1 µs per loop

In [11]: %%timeit
    ...: unq, idx=np.unique(np_array, return_inverse=True)
    ...: res2= gammaln(unq)[idx]
10000 loops, best of 3: 63.3 µs per loop

大部分额外的时间是在计算unique；一旦我们有了映射，应用它就会非常快：

In [12]: timeit res2=gammaln(unq)[idx]
...
100000 loops, best of 3: 6.12 µs per loop

求和的唯一性

返回值的总和而不是实际值会稍微移动相对时间

In [13]: %timeit res1= gammaln(np_array).sum()
10000 loops, best of 3: 35.3 µs per loop
In [14]: %%timeit
    ...: unq, idx=np.unique(np_array, return_inverse=True)
    ...: res2= gammaln(unq)[idx].sum()
10000 loops, best of 3: 71.3 µs per loop

使用 Paul Panzer 的计数方法

In [21]: %%timeit 
    ...: unq, cnt=np.unique(np_array, return_counts=True)
    ...: res=(gammaln(unq)*cnt).sum()
10000 loops, best of 3: 59.5 µs per loop
In [22]: %%timeit 
    ...: unq, cnt=np.unique(np_array, return_counts=True)
    ...: res=np.dot(gammaln(unq),cnt)
10000 loops, best of 3: 52.1 µs per loop