python - 对数 Gamma 函数的快速算法

Question

我正在尝试编写一个快速算法来计算 log gamma function 。目前我的实现看起来很幼稚，只是迭代了 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对我的代码进行了计时实现和我的实际上慢了十万倍。所以我正在做一些非常错误或非常天真的事情(可能两者兼而有之)。尽管与 scipy 相比，我的答案至少正确到小数点后 4 位以内。

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

是否有更快、更优化的方法来计算对数 Gamma 函数？我如何理解 scipy 的实现，以便将其与我的实现合并？

Answer 1

函数的运行时间将随着迭代次数线性扩展(达到一定的恒定开销)。因此，减少迭代次数是加快算法速度的关键。虽然预先计算 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，总和的参数变为

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

因此，在 t 中，总和的参数为零直至二阶，如果 t 足够小，则可以忽略不计。换句话说，迭代次数至少应与 z 一样大，最好至少大一个数量级。

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