log(1-exp(x)) を安定に計算する

「ベルヌーイ分布のパラメータ x が与えられたときに，裏が出る確率 1-x を計算したい」「1-x ではなく log(1-x) がほしい」「x ではなく log x が与えられている」という状況で，タイトルのように log(1-exp(x)) を計算したい状況になりました。

機械学習では何かの値の対数を計算することがよくありますが，いくつかの計算については数値計算的に安定な方法が提案されています。有名なのは log Σ_i exp(x_i) を計算する logsumexp です。そのまま exp(x_i) を計算するとオーバーフローする可能性がありますが，最終的には対数をとっているので答えはそこまで大きくならないはずで，exp(x_i) のオーバーフローを回避する方法があります。

log(1-exp(x)) についても類似のテクニックがないかTwitterでお聞きしたところ，数名の方から有用な情報を頂いたので紹介します。

こちらを参考にしてください。https://t.co/0ev1fTVBAp

— 加藤公一 （はむかず） (@hamukazu) 2017年7月11日

Stanでよければlog1m_expという関数があります。そこでGiuHubでstanプロジェクトのmathリポジトリの中から、該当するC++のソースコードを探して読めば、安定に計算する具体的な方法を知れると思います。

— Kentaro Matsuura (@hankagosa) 2017年7月11日

JuliaのStatsFuns.jlにlog1mexpという関数がhttps://t.co/iaYOlkoDKJhttps://t.co/M4awok6u5b

— (｢・ω・)｢ｶﾞｵｰ Julia㌠ (@bicycle1885) 2017年7月11日

今回は Stan のために log(1-exp(x)) が必要だったため，Stan の log1m_exp 関数の中身を調べることにしました。GitHub の stan-dev/math の中で log1m_exp の実装は以下の箇所にあります。

math/log1m_exp.hpp at 8008de4e4d867ae21e6a66d561588a4ebb80e866 · stan-dev/math · GitHub

inline double log1m_exp(double a) {
  using std::log;
  using std::exp;
  if (a >= 0)
    return std::numeric_limits<double>::quiet_NaN();
  else if (a > -0.693147)
    return log(-expm1(a));  // 0.693147 ~= log(2)
  else
    return log1m(exp(a));
}

見たところによると，exp(x) >= 0.5 の場合，すなわち 1-exp(x) が 0 に近い場合は log(-expm1(x)) を計算し，exp(x) < 0.5 の場合，すなわち 1-exp(x) が 1 に近い場合は log1m(exp(x)) を計算しているようです。引き算で生じる桁落ちを防ぎたい気持ちが感じられますね！

ちなみに expm1(x)=exp(x)-1 と log1m(x)=log1p(-x)=log(1-x) の実装はboostのものを使っていて，ソースは以下のページで見られます。

boost/math/special_functions/expm1.hpp - 1.64.0

boost/math/special_functions/log1p.hpp - 1.64.0

数値計算にはいろいろなテクニックがあるので身につけていきたいですね。