ロジスティック回帰からソフトマックス関数と交差エントロピーのまでの導出を簡潔かつ分かりやすくメモ。

ロジスティック回帰のモデル

2クラス分類に用いる回帰の方法。特徴\(x\)がラベル\(0\)である確率\(p_0\)とラベル\(1\)である確率\(p_1\)

\[ \begin{align} p_0&\equiv P(Y=0\;|\;X=x)\\ p_1&\equiv P(Y=1\;|\;X=x) \end{align}\]

に関して、2クラス分類なので\(p_0+p_1=1\)が成立している。このとき次のモデルを考える。

\[ \begin{align} \log p_0&\propto\beta^{(0)}_0+x^\top\beta^{(0)}\\ \log p_1&\propto\beta^{(1)}_0+x^\top\beta^{(1)} \end{align}\]

これは\(p_0+p_1=1\)を用いて\(p_0,p_1\)について解くと

\[ \begin{align} p_0&=\frac{\exp(\beta^{(0)}_0+x^\top\beta^{(0)})} {\exp(\beta^{(0)}_0+x^\top\beta^{(0)})+\exp(\beta^{(1)}_0+x^\top\beta^{(1)})}\\ p_1&=\frac{\exp(\beta^{(1)}_0+x^\top\beta^{(1)})} {\exp(\beta^{(0)}_0+x^\top\beta^{(0)})+\exp(\beta^{(1)}_0+x^\top\beta^{(1)})} \end{align}\]

\(p_1\)に注目する。

\[ \begin{align} p_1&=\frac{\exp(\beta^{(1)}_0+x^\top\beta^{(1)})} {\exp(\beta^{(0)}_0+x^\top\beta^{(0)})+\exp(\beta^{(1)}_0+x^\top\beta^{(1)})}\\ &=\frac{1} {1+\exp(\beta^{(0)}_0+x^\top\beta^{(0)})\exp(-\beta^{(1)}_0-x^\top\beta^{(1)})}\\ &=\frac{1} {1+\exp(\beta^{(0)}_0-\beta^{(1)}_0+x^\top(\beta^{(0)}-\beta^{(1)}))} \end{align}\]

ここで\(\beta_0\equiv\beta^{(0)}_0-\beta^{(1)}_0\)、\(\beta\equiv\beta^{(0)}-\beta^{(1)}\)と置き直す。

\[ p_1=\frac{1} {1+\exp(\beta_0+x^\top\beta)}\]

これがいわゆるロジスティック回帰のモデル。\(p_0=1-p_1\)から必然的にもう片方は決まる。もしくは対称性からどっちに注目しても本質的に同じ結果が得られる。

ロジスティック回帰のモデル微分

ロジスティック回帰のモデルはシグモイド関数と同じ形をしている。

\[ \sigma(x)=\frac{1}{1+\exp(-ax)}\]

この微分に関して次が成り立つ。

\[ \begin{align} \sigma'(x)&=\frac{-a\exp(-ax)}{(1+\exp(-ax))^2}\\ &=a\cdot\frac{1}{1+\exp(-ax)}\cdot \frac{-\exp(-ax)}{1+\exp(-ax)}\\ &=a\cdot\frac{1}{1+\exp(-ax)}\cdot \frac{1-(1+\exp(-ax))}{1+\exp(-ax)}\\ &=a\cdot\frac{1}{1+\exp(-ax)}\cdot \left(\frac{1}{1+\exp(-ax)}-1\right)\\ &=a\sigma(x)(\sigma(x)-1) \end{align}\]

これを踏まえると以下がそれぞれ成立。

\[ \begin{align} \frac{\partial p_1}{\partial\beta_0}&=p_1(1-p_1)\\ \frac{\partial p_1}{\partial\beta}&=p_1(1-p_1)x \end{align}\]

このように入力と出力だけの簡単な計算からモデルの勾配が求まる。

ロジスティック回帰のパラメータ推定

\(y_i,x_i,\;(i=1,\dots,n)\)というデータセットが与えられている。\(y_i\in\{0,1\}\)で\(x\in\mathbb{R}^d\)とする。このデータセットからパラメータ\((\beta_0,\beta)\in\mathbb{R}\times\mathbb{R}^d\)を求めたいとする。確率を考えているのだから最尤推定量を用いるのが妥当だろう。尤度はデータが生成される同時確率\(L(\beta_0,\beta)\)、対数尤度\(\ell(\beta_0,\beta)=\log L(\beta_0,\beta)\)は、

\[ \ell(\beta_0,\beta)=\sum_{i=1}^n\log P(Y=y_i\;|\;X=x_i)\]

特に\(y_i=\{0,1\}\)を利用するとこれは\(p_i\equiv p_1(\beta_0,\beta,x_i)\)を用いて、

\[ \ell(\beta_0,\beta)=\sum_{i=1}^n(y_i\log p_i+(1-y_i)\log (1-p_i))\]

これこれを最大化するような\(\beta_0,\beta\)によってモデルのパラメータを決定する。これは符号を逆にして最小化するのと同じ。対数尤度はエントロピー関数と同じ形をしていてエントロピーが凹関数（concave）であることとシグモイド関数の単調性から符号を逆にした対数尤度関数は凸関数になり、勾配降下法による最適化が使える。

\[ \begin{align} \frac{\partial\ell(\beta_0,\beta)}{\partial\beta_0} &=\sum_{i=1}^n\frac{\partial}{\partial p_i}(y_i\log p_i+(1-y_i)\log (1-p_i))\cdot \frac{\partial p_i}{\partial \beta_0}\\ &=\sum_{i=1}^n\left(\frac{y_i}{p_i}-\frac{1-y_i}{1-p_i}\right)\cdot p_i(1-p_i)\\ &=\sum_{i=1}^n\frac{y_i(1-p_i)-p_i(1-y_i)}{p_i(1-p_i)}\cdot p_i(1-p_i)\\ &=\sum_{i=1}^n(y_i-p_i)\\ \frac{\partial\ell(\beta_0,\beta)}{\partial\beta} &=\sum_{i=1}^nx_i(y_i-p_i) \end{align}\]

という比較的簡単な形の勾配を用いて勾配降下法を使えばロジスティック回帰のモデルパラメータは決定できる。

多クラス分類問題のロジスティック回帰

ここまでの説明は多クラス分類問題の場合にも容易に拡張できる。\(K\)クラス分類を考えて同じように以下を考える。

\[ \begin{align} p_1&\equiv P(Y=1\;|\;X=x)\\ p_2&\equiv P(Y=2\;|\;X=x)\\ &\vdots\\ p_K&\equiv P(Y=K\;|\;X=x)\\ \end{align}\]

\(K-1\)がうるさいので\(1\)スタートにした。次に

\[ \begin{align} \log p_1&\propto \beta^{(1)}_0+x^\top\beta^{(1)}\\ \log p_2&\propto \beta^{(2)}_0+x^\top\beta^{(2)}\\ &\vdots\\ \log p_K&\propto \beta^{(K)}_0+x^\top\beta^{(K)}\\ \end{align}\]

これをそれぞれについて解く。並べるのは面倒なので\(k=1,\dots,K\)に対して一気にいく。

\[ p_k=\frac{\exp(\beta^{(k)}_0+x^\top\beta^{(k)})}{\sum_{j=1}^K\exp(\beta^{(j)}_0+x^\top\beta^{(j)})}\]

これはソフトマックス関数といわれているものと同じ。多クラスのロジスティック回帰からも導かれることがわかる。多クラスの場合はある1つに着目して分母分子割っても見通しが良くならないのでこのままいく。

ソフトマックスの微分

シグモイドを微分したのと同じようにソフトマックス関数の微分を計算しておく。パラメータの組が\(k=1,\dots,K\)だけあり、\(\partial p_k/\partial \beta^{(\ell)}\)は\(k=\ell\)と\(k\neq\ell\)の場合に分けて考えることが出来る。

\(k=\ell\)のとき、

\[ \begin{align}\frac{\partial p_k}{\partial\beta^{(\ell)}}&=\frac{x\exp(\beta^{(k)}_0+x^\top\beta^{(k)})}{\sum_{j=1}^K\exp(\beta^{(j)}_0+x^\top\beta^{(j)})}-\frac{x\exp(\beta^{(k)}_0+x^\top\beta^{(k)})^2}{\left(\sum_{j=1}^K\exp(\beta^{(j)}_0+x^\top\beta^{(j)})\right)^2}\\&=x\cdot\frac{\exp(\beta^{(k)}_0+x^\top\beta^{(k)})}{\sum_{j=1}^K\exp(\beta^{(j)}_0+x^\top\beta^{(j)})}\left(1-\frac{\exp(\beta^{(k)}_0+x^\top\beta^{(k)})}{\sum_{j=1}^K\exp(\beta^{(j)}_0+x^\top\beta^{(j)})}\right)\\&=xp_k(1-p_k)\\\frac{\partial p_k}{\partial\beta^{(\ell)}_0}&=p_k(1-p_k)\end{align}\]

\(k\neq\ell\)のとき

\[ \begin{align}\frac{\partial p_k}{\partial\beta^{(\ell)}}&=-\frac{x\exp(\beta^{(k)}_0+x^\top\beta^{(k)})\exp(\beta^{(\ell)}_0+x^\top\beta^{(\ell)})}{\left(\sum_{j=1}^K\exp(\beta^{(j)}_0+x^\top\beta^{(j)})\right)^2}\\&=-x\cdot\frac{\exp(\beta^{(k)}_0+x^\top\beta^{(k)})}{\sum_{j=1}^K\exp(\beta^{(j)}_0+x^\top\beta^{(j)})}\cdot\frac{\exp(\beta^{(\ell)}_0+x^\top\beta^{(\ell)})}{\sum_{j=1}^K\exp(\beta^{(j)}_0+x^\top\beta^{(j)})}\\&=-xp_kp_\ell\\\frac{\partial p_k}{\partial\beta^{(\ell)}}&=-p_kp_\ell\end{align}\]

これを\(\delta_{k\ell}\)をつかってまとめると

\[ \begin{align}\frac{\partial p_k}{\partial \beta^{(\ell)}}&=xp_k(\delta_{k\ell}-p_\ell)\\\frac{\partial p_k}{\partial \beta^{(\ell)}_0}&=p_k(\delta_{k\ell}-p_\ell)\\&\delta_{k\ell}=\left\{\begin{matrix}1 & (k=\ell)\\0 & (k\neq\ell)\end{matrix}\right.\end{align}\]

と書ける。シグモイド関数と同じく入力と出力を用いて簡素な形で書ける。

交差エントロピーとパラメータ推定

ソフトマックスときたら次は交差エントロピーを導出する。同じようにして対数尤度を最大化することを考える。

\[ \ell(\beta^{(1)}_0,\dots,\beta^{(K)}_0,\beta^{(1)},\dots,\beta^{(K)})=\sum_{i=1}^n\log P(Y=y_i\;|\;X=x_i)\]

この場合は\(y_i=1,\dots,K\)だとうまいこと表せないので次のone-hot表現を使う。

\[ t_{ik}=\left\{\begin{matrix}1 & (y_i=k)\\0 & (y_i\neq k)\end{matrix}\right.\]

これと\(p_{ik}\equiv p_k(\beta^{(1)}_0,\dots,\beta^{(K)}_0,\beta^{(1)},\dots,\beta^{(K)},x_i)\)使うと

\[ \ell(\beta^{(1)}_0,\dots,\beta^{(K)}_0,\beta^{(1)},\dots,\beta^{(K)})=\sum_{i=1}^n\sum_{k=1}^K t_{ik}\log p_{ik}\]

となる。この符号を逆にしたものが交差エントロピー損失（cross-entropy loss）に他ならない。これもソフトマックスの微分を用いて損失関数の勾配が求められる。

\[ \begin{align}\frac{\partial \ell}{\partial\beta^{(\ell)}}&=\sum_{i=1}^n\sum_{k=1}^K\frac{\partial}{\partial p_{ik}}(t_{ik}\log p_{ik})\cdot\frac{\partial p_{ik}}{\partial \beta^{(\ell)}}\\&=\sum_{i=1}^n\sum_{k=1}^K\frac{t_{ik}}{p_{ik}}\cdot x_ip_{ik}(\delta_{k\ell}-p_{i\ell})\\&=\sum_{i=1}^nx_i\sum_{k=1}^K t_{ik}(\delta_{k\ell}-p_{i\ell})\\&=\sum_{i=1}^nx_i\left(\sum_{k=1}^K t_{ik}\delta_{k\ell}-p_{i\ell}\sum_{k=1}^K t_{ik}\right)\\&=\sum_{i=1}^nx_i(t_{i\ell}\delta_{\ell\ell}-p_{i\ell})\\&=\sum_{i=1}^nx_i(t_{i\ell}-p_{i\ell})\\\frac{\partial \ell}{\partial\beta^{(\ell)}_0}&=\sum_{i=1}^n(t_{i\ell}-p_{i\ell})\end{align}\]

for \(\ell=1,\dots,K\)、というように求められる。このようにしてもとめた勾配を用いて勾配降下法をやれば多クラスロジスティック回帰の解が求まる。

数値計算上のTips

ソフトマックスと交差エントロピー損失の計算は素朴にやるとオーバーフローが発生する。そこで、ソフトマックスの計算には次を使う。

\[ \begin{gather}p_k=\frac{\exp(\beta^{(k)}_0+x^\top\beta^{(k)}-C)}{\sum_{j=1}^K\exp(\beta^{(j)}_0+x^\top\beta^{(j)}-C)}\\C=\max_{j=1,\dots,K}(\beta^{(k)}_0+x^\top\beta^{(k)})\end{gather}\]

このようにすると各項の最大値は\(1\)となるのでオーバーフローは発生しない。分母分子を\(\exp(-C)\)で割れば値としては元のソフトマックスと変わらない。

次に交差エントロピー損失の計算について、上記のようにしてソフトマックスを計算すると計算機上\(p_k=0\)となることがありえる。これは指数部があまりに違い過ぎるもの同士の演算で小さい側の値が結果に一切反映されないという情報落ち（正確には指数部のオーバーフローも）が発生するためだ。

このとき\(\log p_k\)は負の無限大に発散するのでこの部分の計算結果は-infとなり、これと\(t_k=0\)を掛け算すると不定形となり結果はnanとなる。これを避けるために十分小さい\(\varepsilon\)を用いて以下を用いる（\(10^{-12}\)など）。

\[ t_k\log(\varepsilon+(1-\varepsilon)p_k)\]

くわえて、\(t_k=0\)がわかっているならその部分は最終的な損失に影響しないから対応する\(\log p_k\)の計算を行わないようにする。\(\varepsilon\)と内分しておくのはあくまで保険ととらえる。

これで計算機上でも安定してクロスエントロピー誤差を計算できる。もしクロスエントロピー周りでオーバーフローやnanが発生したらばこれを試すとよい。

参考文献

Trevor Hastie, Robert Tibshirani, Jerome Friedman, The Elements of Statistical Learning, Data mining, Inference, and Prediction, Second Edition, 2017, p. 119, 4.4 Logistic Regression.