Word2Vecモデルをスクラッチで実装してみる ② 基本のNeural Netの実装

このブログは情報系を勉強する女子大生 Advent Calendar 2017 - Qiitaの17日目の記事です。

内容としては Word2Vecモデルをスクラッチで実装してみる ① そもそもWord2Vecって？ - あさりさんの作業ログの記事の続きになります。
前回はざっくりWord2Vecモデルがどんなものか説明したので、いよいよ実装してみようと思います。
基本的にNumpyを使って実装しました。ちなみに諸事情によりPython2系です。
実際にモデルをStanford Sentiment Treebankのデータで学習させた結果を二次元状にプロットすると以下の感じになります。

f:id:akaringo030402:20171231151856p:plain

~~（結果が微妙？？？キニシナイキニシナイ）~~

実装にあたっては前回も紹介したスタンフォード大学の Stanford CS 224N | Natural Language Processing with Deep Learningの授業資料や授業課題を参考にしています。
とても良い授業なのでお正月に暇な方は見てみてください〜 :)

ニューラルネットワーク自体の説明については、自分が拙く行うよりもネットや教科書でとてもわかりやすくまとめてくださっている方がたくさんいるので割愛します。
個人的に日本語でネットで読めるものだと愛媛大学の村上研究室のニューラルネットワークについての記事の第3, 4, 5章がわかりやすいと思いました。

この記事では実装や勾配計算、勾配の確認などをざっくりさらっていきたいと思います('-'*)

基本のNeural Netを実装してみる

今回は隠れ層が一層のみのニューラルネットワークを実装してみます。 f:id:akaringo030402:20171231160255p:plain

活性化関数にはシグモイド関数を、出力層にはソフトマックス関数、コスト関数は交差エントロピー損失を用います。
この記事ではまずソフトマックス関数、シグモイド関数、交差エントロピー関数について実装と勾配を計算した後、順伝播、逆伝播部分の実装を行なって最後に実際に逆伝播の実装に誤りがないか、gradient chekingを行う構成になっています。

ソフトマックス関数を実装する

まずソフトマックス関数を実装してみます。
ソフトマックス関数とは $n$ 次元実数ベクトル $x=(x_1,\ \ldots,\ x_n)$ を受け取って以下を満たす $n$ 次元実数ベクトル $y = (y_1,\ \ldots,\ y_n)$ を返す関数です。

${\displaystyle 　y_k = \frac{\exp{(a_k)}}{\sum_{i=1}^n \exp{(a_i)}} }$

式からもわかるようにソフトマックス関数には

$0 \leq y_i \leq 1$
$\sum_{i=1} ^ n y_i = 1$

という性質があり、分類問題をニューラルネットワークで解く場合に、出力層の活性化関数として用いられます。
例えば手書き文字画像からその文字が実際にどの数字を指しているのか当てるmnist問題では、入力が実際に0から9のどの数字になりそうかの確率を出力し、最も確率の高い数字を予測として出力します。
この場合出力層で最後の隠れ層の出力結果をsoftmax関数に通すことで、「全体の合計が1になるかつそれぞれが0から1の間の数字になる」10次元の実数ベクトルが帰ってくるので、その中で値が最大になる $y_i$ のインデクスが分類予測結果になります。

とりあえず実装してみる

とりあえずnumpyを使ってこのsoftmax関数を実装してみます。

def softmax(x):
　　e_x = np.exp(x)
　　sum_e_x = np.sum(e_x)
　　x = e_x / sum_e_x
　　return x

これみるとちょっとゾワってしません？？ちなみにこのコードで適当にsoftmax(np.array([1, 3, 4, 5, 1000]))と入れて実行してみるとオーバーフローのために出力値が不定値になっていることがわかります。

__main__:2: RuntimeWarning: overflow encountered in exp
__main__:4: RuntimeWarning: invalid value encountered in true_divide
array([  0.,   0.,   0.,   0.,  nan])

ちなみにsoftmax(np.array([1, 3,4, 5, 10]))を実行すると、softmax関数によって合計が1になるような実数値のベクトルが出力されていることがわかります。

array([  1.22157447e-04,   9.02628229e-04,   2.45359791e-03,
         6.66957062e-03,   9.89852046e-01])

理由としてはとてもシンプルで、 $\exp{(x)}$ がxが大きくなると簡単に数値がオーバーフローして不定値になってしまうためです。
試しにnp.exp(1000)を計算してみるとオーバーフローが発生していることが確認できます。

>>> np.exp(1000)
__main__:1: RuntimeWarning: overflow encountered in exp

オーバーフローをしないよう最大値を引く

この対策としてよく行われるのが、あらかじめ入力として与えられた実数値ベクトルにおける最大値を引区という方法で、数式で表すと以下のような計算を行うことになります。
ソフトマックス関数の入力に対して定数オフセットを追加しても（この場合は入力から最大の値を引いても）結果は不変になります。

$\exp{(ab)} = \exp{(a)} \exp{(b)}$ だったりすることを利用して~~工学部なのでざっくり~~と計算してみると、実際に定数 $C$ を追加した場合も出力結果は不変であることがわかります。

${\displaystyle 　y_k' = \frac{\exp{(a_k + C)}}{\sum_{i=1}^n \exp{(a_i + C)}} = \frac{\exp{(C)} \exp{(a_k)}}{\exp{(C)} \sum \exp{(a_i)}} = \frac{\exp{(a_k)}}{\sum \exp{(a_i)}} = y_k }$

実際に入力から最大値を引いたソフトマックス関数をnumpyで実装してみると以下の通りになります。今後のことも考えて、入力が $n$ 次元のベクトルの場合と $n \times m$ のマトリックスである場合両方実装します。

def softmax(x):
    x = x.astype(np.float64)

    if len(x.shape) > 1:
        # Matrix
        e_x = np.exp(x.T - x.max(1))
        x = e_x / e_x.sum(axis=0)
        x = x.T
    else:
        # Vector
        e_x = np.exp(x - np.max(x))
        x = e_x / e_x.sum(axis=0)

    return x

シグモイド関数及び勾配を計算して実装してみる

次に活性化関数のシグモイド関数の実装について考えてみます。
シグモイド関数は以下の式で表され、 $(-\infty ,\infty )\rightarrow (0,1)$ の単調増加連続関数で、1つの変曲点を持つ実関数（つまりどんな大きな入力が与えられても、小さな入力が与えられても出力結果が0から1の間の実数になる）です。

${\displaystyle 　\sigma (x) = \frac{1}{1 + \exp{(-x)}} }$

ソフトマックスの時と同様、numpyを使って実装してみると以下のようになります。

def sigmoid(x):
   x = np.array(x, dtype=np.float128)
   s = 1.0 / (1 + np.exp(-x))
   return s

学習中のオーバーフロー対策としてとりあえずnp.float128でキャスティングしているのですが良い子は真似しないでください。
シグモイド関数でのオーバーフロー対策としては値をクリッピングする（np.float64で表現できる値より大きくなる場合は強制的に値を打ち止める）やscipy.special.explitを使う方法があるみたいです。

このシグモイド関数の勾配についても実装します。
シグモイド関数の勾配は商の微分公式など使って計算でき、シグモイド関数の勾配はシグモイド関数の値と1からその値を引いた値との積で表せることがわかります。

${\displaystyle 　\sigma' (x) = \frac{-(- \exp{(-x)})}{(1 + \exp{(-x)})^ 2} = \frac{\exp{(-x)}}{(1 + \exp{(-x)})^ 2} = \frac{1}{(1 + \exp{(-x)})} \frac{\exp{(-x)}}{1 + \exp{(-x)}} = \sigma(x)(1 - \sigma(x)) }$

このためシグモイド関数の実装はこんな感じでとてもシンプルに書くことができます。

def sigmoid_grad(s):
    ds = s * (1 - s)
    return ds

コスト関数（交差エントロピー損失）の勾配を計算して実装する

コスト関数として用いる交差エントロピー損失についても、実装と勾配の計算を行います。
交差エントロピー損失についての説明及びなぜ交差エントロピー損失を使うべきかなどは本やブログで詳細な説明を行っている方もいるのでここでは割愛します。

$y$ を正解ラベル、 $\hat{y}$ を予測された結果とすると、交差エントロピー関数は以下のように定義されます。

${\displaystyle 　CE (y, \hat{y}) = -\sum_{i} y_i \log{(\hat{y_i})} }$

今回は出力層にsoftmax関数を使っているので $\hat{y} = softmax{(\theta)}$ となります。

def cross_entropy_loss(y, y_hat): 
   # y = labels, y_hat = softmax(theta)
   cross_extropy_loss = -1 * np.sum(y * np.log(y_hat))
   return cross_extropy_loss

この交差エントロピー損失について、入力 $\theta$ についての勾配を計算すると

${\displaystyle 　\frac{\partial CE (y, \hat{y})}{{\partial \theta}} = \hat{{\bf y}} - {\bf y} }$

となるため、単純に正解ラベルと出力結果との差を計算すれば良いことがわかります。

隠れ層一層のニューラルネットワークの勾配を計算する

上で計算・実装したものも組み合わせながら、「隠れ層が1つのニューラルネット」を実装していきます。

f:id:akaringo030402:20171231160255p:plain

まず順伝播部の実装を行います。上の図の $h$ 及び $\hat{y}$ については次の式で表すことができます。

${\displaystyle 　{\bf h} = {\rm sigmoid}({\bf xW_1} + b_1), \ \hat{{\bf y}} = {\rm softmax} ({\bf hW_2} + b_2) }$

${\bf W_i, b_i} (i=1,2)$ については二層のウェイト及びバイアスを示しています。

まず入力データ、初期パラメータ、入力層・隠れ層・出力層のデータを入力とし、隠れ層及び出力層の出力を返す forward_pop()という関数として実装します。

def forward_prop(data, params, dimensions):
   # dataは入力データ(x)
   # dimentionsがそれぞれの層の次元を表している。
   # paramsは各パラメータの初期値。
   # params = np.random.randn((dimensions[0] + 1) * dimensions[1] + (dimensions[1] + 1) * dimensions[2], ) などで初期化
   # Dxが入力層の次元、Hが隠れ層の次元、Dyが出力層の次元
   ofs = 0
   Dx, H, Dy = (dimensions[0], dimensions[1], dimensions[2])

   W1 = np.reshape(params[ofs:ofs+ Dx * H], (Dx, H)) # W1.shape = (Dx, H)
   ofs += Dx * H
   b1 = np.reshape(params[ofs:ofs + H], (1, H)) # b1.shape = (1, H)
   ofs += H
   W2 = np.reshape(params[ofs:ofs + H * Dy], (H, Dy)) # W2.shape = (H, Dy)
   ofs += H * Dy
   b2 = np.reshape(params[ofs:ofs + Dy], (1, Dy)) # b1.shape = (1, Dy)

   h = sigmoid(np.dot(data, W1) + b1) # h.shape = (N, h)
   y_hat = softmax(np.dot(h, W2) + b2) # y_hat.shape = (N, Dy)
   
   return h, y_hat

dataは(N, Dx)の行列で与えられ、paramsは((Dx + 1) * Dx + (H + 1) * Dy, 1)の乱数ベクトルでまとめて与えています。
順伝播を計算する際はこのparamsをそれぞれ適切な形にnp.reshape()でreshapeした後、ドット積・加算などを行い、先ほど実装したsigmoid, softmaxなどの関数を持ちいてy, y_hatを求めます。

これで順伝播の計算ができたので、逆伝播についても計算して実装してみます。

${\bf z_1} = {\bf xW_1} + b_1, {\bf z_2} = {\bf xW_2} + b_2$ とおき、コスト関数の入力xについての微分を順に後ろから計算していきます。
基本的には合成関数の微分を使って

${\displaystyle 　\frac{\partial CE (y, \hat{y})}{{\partial {\bf x}}} = \frac{\partial CE (y, \hat{y})}{{\partial {\bf z_2}}} \frac{\partial {\bf z_2}}{{\partial{\bf h}}} \frac{\partial {\bf h}}{{\partial {\bf z_1}}} \frac{\partial {\bf z_1}}{{\partial {\bf x}}} }$

を計算していくだけなのですが、一気にやると間違えそうなので一つ一つ微分を計算していきます ('-'*)
（一応ベクトルは太字で区別するようにしているのですが、面倒になって太字化し忘れているところがあるかもしれません。年明けに直します…）

まず交差エントロピー損失関数の $z_2$ についての微分を計算すると、これは上で求めた ${\displaystyle 　\frac{\partial CE (y, \hat{y})}{{\partial \theta}} = \hat{{\bf y}} - {\bf y} }$ をそのまま使えば良いので、

${\displaystyle 　{\bf \delta_1} = \frac{\partial CE}{\partial {\bf z_2}} = \hat{{\bf y}} - {\bf y} }$ となります。

を使ってコスト関数の隠れ層の出力 ${\bf h}$ に対する微分を計算すると

${\displaystyle 　{\bf \delta_2} = \frac{\partial CE (y, \hat{y})}{{\partial {\bf h}}} = \delta_1 \frac{\partial z_2}{{\partial {\bf h}}} = \delta_1W_2^{\mathrm{T}} = (\hat{{\bf y}} - {\bf y}) {\bf W_2}^{\mathrm{T}} }$

となります。さらに合成関数の微分を進めていって交差エントロピー関数のxについての微分を計算していきます。

${\displaystyle 　{\bf \delta_3} = \frac{\partial CE}{{\partial {\bf z_1}}} = \delta_2 \frac{\partial {\bf h}}{{\partial {\bf z_1}}} = \delta_2 \circ ( \sigma(x)(1 - \sigma(x)) }$

${\displaystyle 　\frac{\partial CE (y, \hat{y})}{{\partial {\bf x}}} = \delta_3 \frac{\partial {\bf z_1}}{{\partial {\bf x}}} = \delta_3 {\bf W_1}^{\mathrm{T}} }$

これで必要な計算はできたので、逆伝播についても実装してみます。先ほどのforward_prop()とまとめて一つの関数にします。

def forward_backward_prop(data, labels, params, dimensions):
    ofs = 0
    Dx, H, Dy = (dimensions[0], dimensions[1], dimensions[2])
    ### forward propagation
    W1 = np.reshape(params[ofs:ofs+ Dx * H], (Dx, H))
    ofs += Dx * H
    b1 = np.reshape(params[ofs:ofs + H], (1, H))
    ofs += H
    W2 = np.reshape(params[ofs:ofs + H * Dy], (H, Dy))
    ofs += H * Dy
    b2 = np.reshape(params[ofs:ofs + Dy], (1, Dy))

    h = sigmoid(np.dot(data, W1) + b1)
    y_hat = softmax(np.dot(h, W2) + b2)

    ### backward propagation
    delta = y_hat - labels
    gradW2 = h.T.dot(delta)
    gradb2 = np.sum(delta, axis = 0)
    delta = delta.dot(W2.T) * sigmoid_grad(h)
    gradW1 = data.T.dot(delta)
    gradb1 = np.sum(delta, axis = 0)

    ### cost_function
    cross_extropy_loss = -1 * np.sum(labels * np.log(y_hat))
    cost = cross_extropy_loss

    ### Stack gradients
    grad = np.concatenate((gradW1.flatten(), gradb1.flatten(),
        gradW2.flatten(), gradb2.flatten()))

    return cost, grad

一応これで基本のニューラルネットワークの実装ができました！

gradient checkingで実装が間違っていないか確かめる

今回については手計算でできる勾配計算が多いのですが、それでもだんだん量が増えてきてちょっと本当に実装があっているのか不安になってきますよね…
そういう時はちゃんとgradient checkingをして勾配計算の実装にミスがないかどうか確かめるのが大事です。

gradient checkingとはある $i$ 番目に着目して以下の近似式で勾配を求め、その結果とbackpropで求めた勾配がだいたい一致するかを確かめる作業です。

${\displaystyle 　\frac{\partial J(\theta)}{{\partial {\theta_i}}} = \frac{J(\theta_i + h) -J(\theta_i - h) }{2\times h} }$

これが大幅にbackpropの結果とずれてたら、『あっ…(察し)』とどこかで実装がずれていることがわかるので、特に大規模なネットワークを自分で全部書くときは役に立つ~~のかもしれない~~。

実際にこれを実装してみます。

def gradcheck_naive(f, x):
    rndstate = random.getstate()
    random.setstate(rndstate)

    fx, grad = f(x) # Evaluate function value at original point
    h = 1e-4

    # xの全ての次元についてgradient checkingを行う
    it = np.nditer(x, flags=['multi_index'], op_flags=['readwrite'])
    count = 0
    while not it.finished:
        count+=1
        ix = it.multi_index

        # 近似により勾配を求める。
        x[ix] += h
        random.setstate(rndstate)
        fx_plus_h,_ = f(x)
        random.setstate(rndstate)
        x[ix] -= 2*h
        fx_minus_h,_ = f(x)
        x[ix] += h
        numgrad = (fx_plus_h - fx_minus_h) / (2*h)

        # backpropの計算結果と近似した結果の差異を比較する。
        reldiff = abs(numgrad - grad[ix]) / max(1, abs(numgrad), abs(grad[ix]))
        print reldiff

        if reldiff > 1e-5:
            print "Gradient check failed."
            return
        it.iternext() # つぎのdimentionへ

    print "Gradient check passed!"

これで先ほどのforward_backward_prop()に適当なデータを入れて計算がちゃんとできているか確かめてみます。

def sanity_check():
    N = 20
    dimensions = [10, 5, 10]
    data = np.random.randn(N, dimensions[0]) 
    labels = np.zeros((N, dimensions[2]))
    for i in xrange(N):
        labels[i, random.randint(0,dimensions[2]-1)] = 1

    params = np.random.randn((dimensions[0] + 1) * dimensions[1] + (
        dimensions[1] + 1) * dimensions[2], )

    print params.shape
    gradcheck_naive(lambda params:
        forward_backward_prop(data, labels, params, dimensions), params)

実行してみると、

9.49366951697e-11
2.59475394025e-10
...
4.01561339736e-10
5.21833463139e-11
7.50241228166e-11
Gradient check passed!

とほぼ差がないことがわかります。基本のニューラルネットの実装については大丈夫そうですね (◍ ´꒳` ◍)b

とりあえず今回までで基本のニューラルネットの勾配計算だったりをして実装してみるところまでできました。
次回（三が日までにはかけるように頑張ります…）は実際に今回実装したものをベースに、skip-gramやnegative samplingを実装して実際に日本語wikipedia記事で単語ベクトルを作成してみたいと思います。