Skip to content

Latest commit

 

History

History
301 lines (227 loc) · 7.81 KB

20230925_api_design_for_gammaln.md

File metadata and controls

301 lines (227 loc) · 7.81 KB
Raw

paddle.gammaln 设计文档

API名称 paddle.gammaln
提交作者 汪昕(GreatV)
提交时间 2023-09-25
版本号 V1.0
依赖飞桨版本 develop
文件名 20230925_api_design_for_gammaln.md

一、概述

1、相关背景

gammaln 函数返回 gamma 函数绝对值的自然对数，定义为：

$$\text{ln}(|\Gamma(x)|)$$

其中 $\Gamma$ 为伽马函数。

2、功能目标

为 Paddle 新增 paddle.gammaln API，提供 gamma 函数绝对值的自然对数计算功能。

3、意义

为 Paddle 新增 paddle.gammaln API，提供 gamma 函数绝对值的自然对数计算功能。

二、飞桨现状

对飞桨框架目前不支持此功能，且飞桨框架暂无 paddle.gamma 函数。但可用其他API组合实现的此功能，代码如下；

import paddle
import numpy as np


def gammaln(x):
    return paddle.log(paddle.abs(paddle.exp(paddle.lgamma(x))))

np_x = np.array([1/5, 1/2, 2/3, 8/7, 3])
x = paddle.to_tensor(np_x, dtype="float64")

out = gammaln(x)
print(out)
# Tensor(shape=[5], dtype=float64, place=Place(cpu), stop_gradient=True,
#        [ 1.52406382,  0.57236494,  0.30315028, -0.06674088,  0.69314718])

需要注意的是，此实现在输入x为大值时，会出现数值溢出的情况，如下所示：

np_x = np.array([1e10, 1e20, 1e40, 1e80])
x = paddle.to_tensor(np_x, dtype="float64")
out = gammaln(x)
print(out)
# Tensor(shape=[4], dtype=float64, place=Place(cpu), stop_gradient=True,
#        [inf., inf., inf., inf.])

这与我们期望的输出有差异

# Tensor(shape=[4], dtype=float64, place=Place(cpu), stop_gradient=True,
#        [220258509288.81057739                                                               ,
#         4505170185988091674624.                                                             ,
#         911034037197618268983181525026762056007680.                                         ,
#         18320680743952364242533014047349914338837359189968031568576067731446620442113081344.])

三、业内方案调研

1. Scipy

在 Scipy 中使用的 API 格式如下：

scipy.special.gammaln(x, out=None) = <ufunc 'gammaln'>

其中，x 为 ndarray 类型, 为实数参数。out 为 ndarray 类型，是 gamma 函数绝对值的自然对数的输出。Scipy 中的实现与 Python 标准库中的 math.lgamma 一致。

scipy 中通过 cephes 库实现，代码如下：

/* Logarithm of Gamma function */
double lgam(double x)
{
    int sign;
    return lgam_sgn(x, &sign);
}

double lgam_sgn(double x, int *sign)
{
    double p, q, u, w, z;
    int i;

    *sign = 1;

    if (!cephes_isfinite(x))
 return x;

    if (x < -34.0) {
 q = -x;
 w = lgam_sgn(q, sign);
 p = floor(q);
 if (p == q) {
   lgsing:
     sf_error("lgam", SF_ERROR_SINGULAR, NULL);
     return (INFINITY);
 }
 i = p;
 if ((i & 1) == 0)
     *sign = -1;
 else
     *sign = 1;
 z = q - p;
 if (z > 0.5) {
     p += 1.0;
     z = p - q;
 }
 z = q * sin(M_PI * z);
 if (z == 0.0)
     goto lgsing;
 /*     z = log(M_PI) - log( z ) - w; */
 z = LOGPI - log(z) - w;
 return (z);
    }

    if (x < 13.0) {
 z = 1.0;
 p = 0.0;
 u = x;
 while (u >= 3.0) {
     p -= 1.0;
     u = x + p;
     z *= u;
 }
 while (u < 2.0) {
     if (u == 0.0)
  goto lgsing;
     z /= u;
     p += 1.0;
     u = x + p;
 }
 if (z < 0.0) {
     *sign = -1;
     z = -z;
 }
 else
     *sign = 1;
 if (u == 2.0)
     return (log(z));
 p -= 2.0;
 x = x + p;
 p = x * polevl(x, B, 5) / p1evl(x, C, 6);
 return (log(z) + p);
    }

    if (x > MAXLGM) {
 return (*sign * INFINITY);
    }

    q = (x - 0.5) * log(x) - x + LS2PI;
    if (x > 1.0e8)
 return (q);

    p = 1.0 / (x * x);
    if (x >= 1000.0)
 q += ((7.9365079365079365079365e-4 * p
        - 2.7777777777777777777778e-3) * p
       + 0.0833333333333333333333) / x;
    else
 q += polevl(p, A, 4) / x;
    return (q);
}

2. Pytorch

在 Pytorch 中使用的 API 格式如下：

torch.special.gammaln(input, *, out=None) → Tensor

其中，input 为 输入Tensor，out 为 输出 Tensor。

Pytorch中直接使用的是C++标准库中的 std::lgamma，代码如下：

  • CPU版本
  Vectorized<T> map(T (*const f)(T)) const {
    Vectorized<T> ret;
    for (int64_t i = 0; i != size(); i++) {
      ret[i] = f(values[i]);
    }
    return ret;
  }

  Vectorized<T> lgamma() const {
    return map(std::lgamma);
  }
  • GPU版本
  const auto lgamma_string = jiterator_stringify(
  template <typename T>
  T lgamma_kernel(T a) {
    return lgamma(a);
  }
); // lgamma_string

四、对比分析

Scipy 基于 cephes 库实现，Pytorch 基于 C++ 标准库实现。Scipy 中的实现与 Python 标准库中的 math.lgamma 一致。Pytorch 中直接使用的是C++标准库中的 std::lgamma。而飞桨中暂无 paddle.gamma 函数，不好通过Python端实现，因此，需要通过C++端实现。

五、设计思路与实现方案

命名与参数设计

API设计为 paddle.gammaln(x, name=None)。其中，x 为 Tensor 类型。Tensor.gammaln() 为 Tensor 的方法版本。

API实现方案

参考 PyTorch 采用C++ 实现，实现位置为 Paddle repo python/paddle/tensor/math.py 目录。并在 python/paddle/tensor/init.py 中，添加 gammaln API，以支持 Tensor.gammaln 的调用方式。头文件放在 paddle/phi/kernels 目录，cc 文件在 paddle/phi/kernels/cpu 目录， cu文件 paddle/phi/kernels/gpu 目录。

其中，Gammaln 的前向计算实现核心内容可为：

template <typename T>
struct GammalnFunctor {
  GammalnFunctor(const T* x, T* output, int64_t numel)
      : x_(x), output_(output), numel_(numel) {}

  HOSTDEVICE void operator()(int64_t idx) const {
    using MT = typename phi::dtype::MPTypeTrait<T>::Type;
    const MT mp_x = static_cast<MT>(x_[idx]);
    output_[idx] = static_cast<T>(std::lgamma(mp_x));
  }

 private:
  const T* x_;
  T* output_;
  int64_t numel_;
};

算子的反向计算实现核心内容可为：

template <typename T>
struct GammalnGradFunctor {
  GammalnGradFunctor(const T* dout, const T* x, T* output, int64_t numel)
      : dout_(dout), x_(x), output_(output), numel_(numel) {}

  HOSTDEVICE void operator()(int64_t idx) const {
    using MT = typename phi::dtype::MPTypeTrait<T>::Type;
    const MT mp_dout = static_cast<MT>(dout_[idx]);
    const MT mp_x = static_cast<MT>(x_[idx]);
    output_[idx] = static_cast<T>(mp_dout * digamma<MT>(mp_x));
  }

 private:
  const T* dout_;
  const T* x_;
  T* output_;
  int64_t numel_;
};

六、测试和验收的考量

可考虑一下场景：

  • 输出数值结果的一致性和数据类型是否正确，使用 scipy 作为参考标准；
  • 对不同 dtype 的输入数据 x 进行计算精度检验 (float32, float64)；
  • 大值输入时，是否出现数值溢出；
  • 输入输出的容错性与错误提示信息；
  • 输出 Dtype 错误或不兼容时抛出异常；
  • 保证调用属性时是可以被正常找到的；
  • 覆盖静态图和动态图测试场景。

七、可行性分析和排期规划

本 API 主要参考 PyTorch 实现，难度适中，工期上能满足要求。

八、影响面

为独立新增API，对其他模块没有影响。

名词解释

附件及参考资料

  1. 伽马函数的详细介绍参见 Gamma Function