@beam2dさんの記事が面白かったのだが、これを読んでよくわからなかった人もそれなりにいるだろうという想定で、補足説明のつもりでこれを書いている。
いろんなところで話をするのだが、「数学≠数値計算」である。数学の知識と数値計算の知識はちょっと違う。方程式の解が解析的に解けて、やった解けた！と思っても、その値をコンピュータで計算して結果を出すまでにさらに工夫が必要なことがよくある。コンピュータ内部での数値表現の特性を理解した上で、コンピュータフレンドリーな数式にさらに変形しなければいけないことがある。
例1: Sigmoid関数の計算
機械学習でよく出てくるSigmoid関数
$$ \sigma(x) = \frac{1}{1+e^{-x}} $$
であるが、これは\( \lim_{x\to \infty} \sigma(x) =1 \), \( \lim_{x\to -\infty} \sigma(x) =0 \)という性質がある。これは、この式のまま実装すると、\( x \)が小さい時にうまくいかない。
実際にPythonでやってみる。
>>> def sigmoid_(x): ... return 1/(1+np.exp(-x)) ... >>> sigmoid_(10) 0.99995460213129761 >>> sigmoid_(-10) 4.5397868702434395e-05 >>> sigmoid_(-1000) __main__:2: RuntimeWarning: overflow encountered in exp 0.0
これは、np.expのカッコの中の絶対値が大きすぎてオーバーフローしてしまったことが原因である。しかし、expの計算ができなくてもsigmoid関数は0〜1の値なので計算できないとおかしい。
これを解決するにはいろんな方法が考えられるが、最も単純には、\(x\)がある値\(a\)（ただし、aは負で、exp(-\(a\))がオーバーフローしないようなぎりぎりの値）より小さいときは0にしてしまえばよい。つまり次のように定義すればよい。
$$
\sigma(x)=
\left\{
\begin{array}{ll}
0 & (x \lt a)\\
\frac{1}{1+e^{-x}} & (x\geq a)
\end{array}
\right.
$$
（ただし、\(a\)はexp(-\(a\))がオーバーフローしない程度に十分小さい値）
>>> def sigmoid(x): ... return 0.0 if x<-709 else 1/(1+np.exp(-x)) ... >>> sigmoid(10) 0.99995460213129761 >>> sigmoid(-10) 4.5397868702434395e-05 >>> sigmoid(-1000) 0.0
ここで、\(a\)として-709を使った。実際sigmoid(-709)を計算すると結果は1.2167807506234229e-308となり、これはほぼ0である。
例2: Softplus関数の計算
Softplus関数
$$ \tau(x)= \log (1+e^x) $$
の計算をしたいとする。これもこの式のまま計算すると\(x\)が大きい時にexpがオーバーフローするが、一方で関数の真の値はほぼ\(x\)に近い値になる。
ここで@beam2dさんのブログに書かれていたとおり、
\[
\begin{align}
\tau(x)&= \left\{
\begin{array}{ll}
x+\log (1+e^{-x}) & (x > 0) \\
\log (1+e^x) & (x \leq 0)
\end{array}\right.\\
&=\max(0,x)+\log (1+e^{-|x|})
\end{align}
\]
と変形してから実装すればいい。実際、式変形する前と後での動作の違いを見てみるとわかるかと思う。
>>> def softplus_(x): ... return np.log(1+np.exp(x)) ... >>> softplus_(1000) inf >>> def softplus(x): ... return np.maximum(0,x)+np.log(1+np.exp(-np.abs(-x))) ... >>> softplus(1000) 1000.0
ここでsoftplus_は元の定義式をそのまま使ったものだが、これに1000を引数に与えるとinfになってしまう。softplusは式変形したものを利用したもので、1000を与えると1000が返ってくる。数学的にはこっちのほうが正しい。
おまけ: NumPy芸
上記のsoftplusの実装は、すべてユニバーサル関数なので、引数が配列のときも有効である。実際やってみると次のようになる。
>>> softplus(np.array([10,100,1000])) array([ 10.0000454, 100. , 1000. ])
上記のsigmoid関数の方は配列にそのまま適用できないが、numpy.whereを使って書き換えれば簡単にユニバーサル関数化できる（読者の演習問題とする）。
まとめ
* 数学的考察で解析的な解が求まったとしても、そこからコンピュータで数値的に計算できるかまでにはギャップがあることがある。
* 関数値が有限でも、計算過程でとても大きな値が発生するケースが要注意で、その場合に数式の同値変形で回避できることがある。
* Pythonで大量に同じ関数を計算することを考えると、場合分けも含めてできるだけユニバーサル関数で書いておいたほうがいい。なにしろ、PythonではFor文書いたら負けなので。
文献案内
数値計算に興味をもったら、是非これを読むことをお勧めする。
古い本だが内容は今でも十分に役に立つ。桁落ちのワナなどの説明がくわしい。ページ数も少ないし読みやすいと思う。
更新履歴：
2015/8/1 softplusが途中からsoftmaxになってたのを直しました。