1

我正在尝试编写一个快速算法来计算对数伽玛函数。目前我的实现似乎很幼稚，只是迭代 1000 万次来计算 gamma 函数的日志（我也在使用 numba 来优化代码）。

import numpy as np
from numba import njit
EULER_MAS = 0.577215664901532 # euler mascheroni constant
HARMONC_10MIL = 16.695311365860007 # sum of 1/k from 1 to 10,000,000

@njit(fastmath=True)
def gammaln(z):
"""Compute log of gamma function for some real positive float z"""
    out = -EULER_MAS*z - np.log(z) + z*HARMONC_10MIL
    n = 10000000 # number of iters
    for k in range(1,n+1,4):
        # loop unrolling
        v1 = np.log(1 + z/k)
        v2 = np.log(1 + z/(k+1))
        v3 = np.log(1 + z/(k+2))
        v4 = np.log(1 + z/(k+3))
        out -= v1 + v2 + v3 + v4

    return out

我根据scipy.special.gammaln实现对我的代码进行了计时，而我的代码实际上慢了 100,000 倍。所以我在做一些非常错误或非常幼稚的事情（可能两者兼而有之）。尽管与 scipy 相比，我的答案至少在小数点后 4 位以内是正确的。

我试图阅读实现 scipy 的 gammaln 函数的 _ufunc 代码，但是我不明白 _gammaln 函数所写的 cython 代码。

有没有更快、更优化的方法可以计算对数伽玛函数？我如何理解 scipy 的实现，以便将其与我的结合起来？

4

3 回答 3

3

您的函数的运行时间将随着迭代次数线性扩展（直到一些恒定的开销）。所以减少迭代次数是加速算法的关键。虽然HARMONIC_10MIL预先计算是一个聪明的想法，但当您截断系列时，它实际上会导致更差的准确性；仅计算系列的一部分结果会提供更高的准确性。

下面的代码是上面发布的代码的修改版本（尽管使用cython代替numba）。

from libc.math cimport log, log1p
cimport cython
cdef:
    float EULER_MAS = 0.577215664901532 # euler mascheroni constant

@cython.cdivision(True)
def gammaln(float z, int n=1000):
    """Compute log of gamma function for some real positive float z"""
    cdef:
        float out = -EULER_MAS*z - log(z)
        int k
        float t
    for k in range(1, n):
        t = z / k
        out += t - log1p(t)

    return out

如下图所示，即使经过 100 次近似，它也能得到一个接近的近似值。

在此处输入图像描述

在 100 次迭代中，它的运行时间与以下数量级相同scipy.special.gammaln：

%timeit special.gammaln(5)
# 932 ns ± 19 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit gammaln(5, 100)
# 1.25 µs ± 20.3 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

剩下的问题当然是要使用多少次迭代。该函数log1p(t)可以扩展为小的泰勒级数t（这与大的极限有关k）。尤其是，

log1p(t) = t - t ** 2 / 2 + ...

这样，对于大k， sum 的参数变为

t - log1p(t) = t ** 2 / 2 + ...

因此，总和的自变量在二阶之前为零，如果足够小t，则可以忽略不计。t换言之，迭代次数应至少与 一样大z，最好至少大一个数量级。

但是，如果可能的话，我会坚持使用scipy经过良好测试的实现。

于 2019-02-24T11:00:19.460 回答
0

关于您之前的问题，我想一个将scipy.special函数包装到 Numba 的示例也很有用。

例子

只要只涉及简单的数据类型（int、double、double*、...），包装 Cython cdef 函数就非常容易和可移植。有关如何调用 scipy.special 函数的文档，请查看 this。您实际需要包装函数的函数名称位于scipy.special.cython_special.__pyx_capi__. 可以使用不同数据类型调用的函数名称被破坏，但确定正确的名称非常容易（只需查看数据类型）

#slightly modified version of https://github.com/numba/numba/issues/3086
from numba.extending import get_cython_function_address
from numba import vectorize, njit
import ctypes
import numpy as np

_PTR = ctypes.POINTER
_dble = ctypes.c_double
_ptr_dble = _PTR(_dble)

addr = get_cython_function_address("scipy.special.cython_special", "gammaln")
functype = ctypes.CFUNCTYPE(_dble, _dble)
gammaln_float64 = functype(addr)

@njit
def numba_gammaln(x):
  return gammaln_float64(x)

在 Numba 中的使用

#Numba example with loops
import numba as nb
import numpy as np
@nb.njit()
def Test_func(A):
  out=np.empty(A.shape[0])
  for i in range(A.shape[0]):
    out[i]=numba_gammaln(A[i])
  return out

计时

data=np.random.rand(1_000_000)
Test_func(A): 39.1ms
gammaln(A):   39.1ms

当然，您可以轻松地并行化此函数并优于 scipy 中的单线程 gammaln 实现，并且您可以在任何 Numba 编译函数中有效地调用此函数。

于 2019-02-24T19:40:42.260 回答
0

通过尝试 numba 的并行模式并主要使用矢量化函数，我设法将性能提高了大约 3 倍（遗憾的是，numba 无法理解numpy.substract.reduce）

from functools import reduce
import numpy as np
from numba import njit

@njit(fastmath=True, parallel=True)
def gammaln_vec(z):
    out = -EULER_MAS*z - np.log(z) + z*HARMONC_10MIL
    n = 10000000

    v = np.log(1 + z/np.arange(1, n+1))

    return out-reduce(lambda x1, x2: x1-x2, v, 0)

时间：

#Your function:
%timeit gammaln(1.5)
48.6 ms ± 1.23 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

#My function:
%timeit gammaln_vec(1.5)
15 ms ± 340 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

#scpiy's function
%timeit gammaln_sp(1.5)
1.07 µs ± 18.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

因此，使用 scipy 的功能会更好。没有 C 代码，我不知道如何进一步分解它

于 2019-02-24T11:47:33.893 回答