Pythonで3層ニューラルネットの勾配を計算
3層ニューラルネットとは、入力層・中間層・出力層の3層によって構成される多層パーセプトロンで、回帰モデルの1つといえる。ここではその勾配の計算方法とその行列表示をまとめる。行列表示をすることによってPythonのnumpyだけをつかった勾配計算を簡単に書くことができる。いまは3層のみ、忘れない内にメモる。
モデル
\(n\)コの入出力ペア\((x_i,y_i)\in\mathbb{R}^p\times\mathbb{R},i=1,\dots,n\)に対して次の形の関数をフィッティングする。\(\varphi(\cdot):\mathbb{R}\to\mathbb{R}\)は活性化関数。シグモイドやらReLUやらなんでもいい。
\[ f(x)=\sum_{j=1}^d\beta^{(2)}_j\varphi(x^\top\beta^{(1)}_j)\]
ここで、\(\beta^{(1)}_j\in\mathbb{R}^p,\beta^{(2)}_j\in\mathbb{R}\)はモデルのパラメータとする。さらに\(d\)は中間層のユニット数。式形が煩雑になるのでここではバイアス項は考えないとする。バイアス項を考えたい場合は\(p\to p+1\)として\(x\to(x;1)\)と置き換えればよいだけ。（実はより多層にしていくとそんな簡単にいかないんだけども理論的にはそうだよね。）
計画行列とパラメータベクトルによる表示
\(n\)コの説明変数を横にして縦に積んだものを\(X\in\mathbb{R}^{n\times p}\)とする。これを計画行列(design matrix)と言ったりすることがあるらしい。さらに、\(n\)コの出力変数を縦に積んだものを\(Y\in\mathbb{R}^n\)とする。同じようにして、パラメータは縦にして横に積む。つまり\(\beta^{(1)}\in\mathbb{R}^{p\times d},\beta^{(2)}\in\mathbb{R}^d\)とする。これを用いると、モデルに対する関数のフィッティングは
\[ Y=\varphi(X\beta^{(1)})\beta^{(2)}+E\]
という関係式の誤差ベクトル\(E\in\mathbb{R}^n\)最小化問題となる。\(\varphi(\cdot)\)の意味は少し変わっていてelement-wiseに作用すると考える。（いや、普通に要素ごとに作用するって言えよ！なんやねんelement-wiseって！）
考えるのは普通に最小二乗回帰として\(\min\|E\|_2^2/n\)とする(平均二乗誤差、Mean Squared Error : MSE)。最適化を考える場合は\(1/n\)はあってもなくてもいい。考えられるシチュエーションとして、\(n\)を振ったりする場合は\(n\)に対して値がもりもり大きくなっていって比較出来なくなることを踏まえると最初から\(1/n\)しておくのが安定。しかし、Pythonのmeanを使わない今回の行列表示による計算では、グローバルな\(n\)は参照しない、とする場合.shape[0]とかlen()とか気持ち悪いことをして\(n\)を取り出さないといけなくなるので微妙な点もある。統一さえされていれば問題はない。
誤差勾配の計算例：線形回帰の場合
3層ニューラルネットの勾配を計算する前に、例として線形モデル\(y=X\beta\)を考える。まず、同じように計画行列による表示を行うと\(Y=X\beta+E\)となり、
\[ \begin{align} \frac{\partial(\|E\|_2^2/n)}{\partial\beta} &=\frac1n\frac{\partial}{\partial\beta}(Y-X\beta)^2\\ &=-2X^\top(Y-X\beta)/n \end{align}\]
となる。ここでは、3層ニューラルネットの場合にもこんな感じで（総和の記号を使わず）行列形式の勾配を求めたい。いっているように理由はPythonでスマートに勾配を計算したいから、総和とか使いたくないねん。
誤差勾配の計算：出力重み
偏微分なので中間層出力が固定されていると考えるとこれは上記の例と全く同じになる。簡単、結果のみ示す。
\[ \frac{\partial(\|E\|_2^2/n)}{\partial\beta^{(2)}} =-2\varphi\left(X\beta^{(1)}\right)^\top\left(Y-\varphi(X\beta^{(1)})\beta^{(2)}\right)/n\]
つまり線形回帰の場合において\(X\to\varphi(X\beta^{(1)})\)と置き換えただけ。
誤差勾配の計算：中間層重み
中間層重みは\(p\times d\)の行列になっているのでベクトルの微分のように簡単にはいかない。まずは結果を示してPythonで計算するところに移る。よって結果。導出は最後の方に記す。
\[ \frac{\partial(\|E\|_2^2/n)}{\partial\beta^{(1)}} =-2X^\top\left\{\left(Y-\varphi(X\beta^{(1)})\beta^{(2)}\right)\left(\beta^{(2)}\right)^\top *\varphi'(X\beta^{(1)})\right\}/n\]
\(\varphi'\)は活性化関数の微分で、例によって（要素ご、いや、）element-wiseに作用する。少しだけ結果の解説。まず、これは行列による微分なので普通の行列の計算だけでは書き表せないのである。いくら式をこねくり回しても無理。途中の\(*\)はアダマール積といわれる演算で、何かっていうと単純に要素ごとに掛け算する演算となる、つまりnumpyでいうところの*による行列積と同じもの（numpyistならアダマール積よりこのほうがわかりやすいよな）。これを踏まえると他は普通に行列の積として整合性は取れている。
Pythonによる勾配計算：線形回帰の例
ここではscipy.optimize.approx_fprimeとscipy.optimize.check_gradによって勾配計算が正しく行われているかを確認する。のだが、例によってまずは簡単な問題から始める。まずはapprox_fprimeとcheck_gradの動作を確認しよう。こういう知らん関数はベクトルの食わせ方なんかを試しておかないとどこにミスがあるのか特定できなくなる。いちいち解説はせんがoptionalのargsでX,Yをつっこんでいるところがポイントか。グローバル使ってもいいんだがまあどっちでも。（お前の引数の解説を見るよりsite:scipy.orgで検索してドキュメントを見るわ！というあるある）
import numpy as np
import scipy.optimize as opt
n, p = 100, 3
X = np.random.randn(n,p)
Y = np.random.randn(n,1)
def varphi(x):
return 1/(1+np.exp(-x))
def varphi_dash(x):
return np.exp(-x)/(1+np.exp(-x))**2
def example_loss(beta, X, Y):
beta = beta.reshape(p, 1)
E = Y - np.dot(X, beta)
return np.mean(E**2)
def example_grad(beta, X, Y):
beta = beta.reshape(p, 1)
return (-2*np.dot(X.T, Y - np.dot(X, beta))/n).flatten()
beta = np.random.randn(p,1).flatten()
err = opt.check_grad(example_loss, example_grad, beta, X, Y)
g1 = opt.approx_fprime(beta, example_loss, 1e-8, X, Y)
g2 = example_grad(beta, X, Y)
print(err, g1, g2)乱数を振っているので具体的な値は毎回変わるが、このコードの結果の一例は以下。
2.683374514175368e-08 [0.48136775 2.1154666 0.99944373] [0.48136772 2.11546664 0.9994437 ]errは近似計算した勾配と、自分で計算した勾配の差のノルムを表しており小さければ正解。今1e-8のオーダと十分小さくなっているのが確認できるので正しく勾配が計算できていることがわかる。g1とg2はそれぞれ近似計算した勾配と自分で計算した勾配の結果であり、両者は一致していることがわかる。これを踏まえて、同じことを3層ニューラルネットの場合に行う。
Pythonによる勾配計算：3層ニューラルネット
コードは次のようになる。注意すべき点は、betaという変数を出力層重みと中間層重みに分割したり逆にpackingしたりするところ。numpyのndarrayは非常に使いやすいのでこれくらいは数行でできる。
import numpy as np
import scipy.optimize as opt
n, p, d = 100, 3, 2
X = np.random.randn(n,p)
Y = np.random.randn(n,1)
beta1 = np.random.randn(p,d)
beta2 = np.random.randn(d,1)
def varphi(x):
return 1/(1+np.exp(-x))
def varphi_dash(x):
return np.exp(-x)/(1+np.exp(-x))**2
def loss(beta, X, Y):
beta1, beta2 = beta[:p*d].reshape(p,d), beta[p*d:].reshape(d,1)
E = Y - np.dot(varphi(np.dot(X, beta1)), beta2)
return np.mean(E**2)
def loss_grad(beta, X, Y):
beta1, beta2 = beta[:p*d].reshape(p,d), beta[p*d:].reshape(d,1)
grad_beta = np.zeros(p*d+d)
grad_beta[:p*d] = (-2*np.dot(X.T, np.dot(Y-np.dot(varphi(np.dot(X, beta1)), beta2), beta2.T)*varphi_dash(np.dot(X, beta1)))/n).flatten()
grad_beta[p*d:] = (-2*np.dot(varphi(np.dot(X, beta1)).T, Y - np.dot(varphi(np.dot(X, beta1)), beta2))/n).flatten()
return grad_beta
beta = np.zeros(p*d+d)
beta[:p*d] = beta1.flatten()
beta[p*d:] = beta2.flatten()
err = opt.check_grad(loss, loss_grad, beta, X, Y)
g1 = opt.approx_fprime(beta, loss, 1e-8, X, Y)
g2 = loss_grad(beta, X, Y)
print(err, g1, g2)同じく乱数振っているが結果の一例は以下。
3.381798920788484e-08
[ 0.00717659 0.03158676 -0.00301397 0.00627223 -0.00108149 0.00943889
-0.47406048 -0.45371247]
[ 0.0071766 0.03158678 -0.00301397 0.00627223 -0.00108145 0.00943891
-0.47406047 -0.45371244]うむ。正しく計算できている。ここでうまくいかなければscipy.optimize.minimizeなんかに損失関数を喰わせても何にもうまくいかない。間違っている場合はg1、g2の中身を見ると中間層重みと出力重みのどちらが間違っているのか、あるいは\(n\)の割忘れなのか、などがわかるのでとにかく一致するまで数式を打ち込むだけ。行列形式で出しているので本当に勾配計算をするのは各層に付き1行で済んでいることに注目。np.dotの連発によりかなり1行は長いがやっていることは単純に正しく式を打ち込んでいるだけ。varphiの作用させ忘れで結構ハマったことを正直にメモしておく。
ついでに言っておくと同じ計算を何回もやっているので計算効率悪いというツッコミも無しで、分かりやすさ重視であえてこう書いています。実際にminimizeを回すときには気休め程度に人間最適化しておけばいいかと。さらについでに言っておくとキャッシュが効くので人間による最適化など不要という人はご自由に。
中間層重み勾配の行列表示の導出
やることは簡単で、行列による微分といっても列ごとに分解してやれば普通のベクトル微分の計算は通用するので\(k=1,\dots,d\)のそれぞれの列に対して勾配を計算して横に並べるだけ。ただしスマートな行列表示にもっていくのは結構苦労する、というか体系的なやり方がわかんないまま結果にたどり着いた。
\[ \begin{align} \frac{\partial\|E\|_2^2/n}{\partial\beta^{(1)}_k} &=-\frac2n \frac{\partial}{\partial\beta^{(1)}_k}\left(\varphi(X\beta^{(1)})\beta^{(2)}\right) \left(Y-\varphi(X\beta^{(1)})\beta^{(2)}\right)\\ &=-\frac2n \left(X^\top\varphi'(X\beta^{(1)})_k\beta^{(2)}_k\right) \left(Y-\varphi(X\beta^{(1)})\beta^{(2)}\right)\\ &=-2X^\top\beta^{(2)}_k\left\{ \varphi'(X\beta^{(1)})_k *(Y-\varphi(X\beta^{(1)})\beta^{(2)})\right\}/n \end{align}\]
まあわかるだろうけど下付きの\(k\)は列の抜き出しを意味します。例えば\(n\times p\)行列であれば第\(k\)列を抜き出した結果\(n\times 1\)の縦ベクトルになるっていう感じです。ただし\(\beta^{(2)}_k\)だけは第\(k\)行の抜き出しで、もともと縦ベクトルだったのでここではもはやスカラーになります。ゆえに前に出しています。あとはこの式とにらめっこしてこの\(p\times 1\)要素を横に並べていくだけ。だけといっても結構うん？となりますが\(\beta^{(2)}_k\)の処理をいい感じにやれば目的の式にはたどり着けます。
以上