OCamlとSLAPで作る型安全ニューラルネット(と深層学習)

  • 59
    いいね
  • 0
    コメント
この記事は最終更新日から1年以上が経過しています。

この投稿は Machine Learning Advent CalendarML Advent Calendar の18日目の記事です.

今日は関数型プログラミング言語 OCaml と線形代数演算ライブラリ SLAP を使った型安全なニューラルネットワークの実装について書きたいと思います.最近,深層学習とかいうニューラルネットの応用が流行っていますし,一方で,関数型プログラミング言語とかいうのも流行っているので,2つの流行に(むりやり同時に)のってみました.

私はOCaml を使って,次元の合わない行列演算をコンパイル時に検出する機能を持った変な線形代数演算ライブラリ Sized Linear Algebra Package (SLAP) を作っています.世の中には便利な線形代数ライブラリ(BLAS とか LAPACK とか)や数値計算言語(MatLab とか R とか)が沢山ありますが,そのほとんどは「2次元ベクトルと3次元ベクトルの足し算」のような次元の合わない変な計算を実行時に検査しているので,例外やメモリ破壊などの実行時エラーが起きます.多くの場合,プログラマは泣きながら,どこで次元が合わなくなったのか調べることになりますが,SLAP を使うとこの手のバグをコンパイル時に検出できます.つまり,2次元ベクトルと3次元ベクトルを加算するプログラムを書くと,OCaml コンパイラが「お前の頭はちょっとおかしい」と教えてくれます.もちろん,ファイルからベクトルを読み込むなど,実行時に初めてサイズが決まるような場合も,ちゃんと検査できます.SLAP 自体の説明は以前書いた記事がわかりやすいと思います.今日はSLAPを使って,実際に機械学習で使う2層ニューラルネットをプログラミングしてみます.

静的サイズ検査の簡単なデモは http://akabe.github.io/slap/#demo に載っています.この記事の最後で,今回作ったニューラルネットのプログラムが静的サイズ検査されることを確認しますが,説明だけでピンとこない人は,先にデモをやってみると良いかもしれません.

インストール

インストールには時間がかかるはずなので,その間に次の節を読んでおくと丁度いいと思います.

まず,SLAP の前に

をインストールします.そして,次のコマンドで SLAP と utop(OCaml の高級対話環境)をインストールします(Emacs の M-x run-ocaml でも対話環境が使えます).

$ opam install slap utop

ニューラルネットワークについて

ロジスティック回帰(1層ニューラルネットワーク)

機械学習において,ロジスティック回帰はクラス分類に使われるモデルです.例えば,音楽データを J-POP とかクラシックに分類する,文書(テキスト)をニュースや小説に分類する,などといったことに使います.このモデルは,分類対象(音楽データやテキスト)を表現する特徴ベクトル(フィルタの応答とか単語の出現頻度など)を入力すると,あるクラス(J-POP とかニュース)に属する確率を返してくれます.つまり,音楽データ(の特徴ベクトル)を入れると,「80% の確率で J-POP です」とか答えてくれます(50% 未満だと J-POP ではない可能性が高い).なかなか,便利ですね!しかも,計算も簡単です.入力となる特徴量 $\boldsymbol{x} = (x_1, x_2, \dots, x_L)^\top$ に対して,出力となる確率 $y$ は

y=\mathrm{sigm}\left(\sum^L_{i=1} w_i x_i + b\right)=\mathrm{sigm}\left(\boldsymbol{w}^\top \boldsymbol{x} + b\right)\quad\mbox{where}\quad\mathrm{sigm}(a)=\frac{1}{1+\exp(-a)}

で与えられます.ここで,$\boldsymbol{w} = (w_1, w_2, \dots, w_L)^\top$ は重み(特徴 $x_i$ の重要度)であり,$\mathrm{sigm}$ はロジスティック・シグモイド関数と呼ばれます.この式は,よく次のような模式図で表されます(真ん中の丸は ユニット と言います).

ロジスティック回帰

1層ニューラルネットだと,ここで紹介したロジスティック回帰や,パーセプトロンが有名です.どちらも非常に簡単に計算できるので,初心者向けのプログラミングの練習問題には丁度いいかもしれませんが,簡単すぎてつまらないので,この記事では触れません(一応,以下のURLにサンプルコードがあります).

2層ニューラルネットワーク

この記事では次の図のような2層ニューラルネットを扱います.ロジスティック回帰モデルを横に並べたもの()を2段組にしただけです.(入力信号も層と考えて,以下のニューラルネットを3層と呼ぶ文献もありますが,ここでは PRML にあわせて2層と表記します.)

two-layer-nn.png

第1層「隠れ層」は入力 $\boldsymbol{x} = (x_1,x_2,\dots,x_L)^\top$ を受け取り,

\boldsymbol{y}=\mathrm{sigmv}(\boldsymbol{W}_1\boldsymbol{x}+\boldsymbol{b}_1)

という式で出力 $\boldsymbol{y} = (y_1,y_2,\dots,y_M)^\top$ を計算します.ただし,$\mathrm{sigmv}$ はベクトルの各要素にシグモイド関数を適用する関数で,$\boldsymbol{W}_1$ は重み行列,$\boldsymbol{b}_1$ はバイアスです.この層の出力 $\boldsymbol{y}$ は外部から観測できないので,隠れ層と呼ばれます.一番上の「出力層」は,この $\boldsymbol{y}$ を入力として,ニューラルネットワーク全体の出力

\boldsymbol{z}=\mathrm{sigmv}(\boldsymbol{W}_2\boldsymbol{y}+\boldsymbol{b}_2)

を計算します.この $\boldsymbol{z} = (z_1,z_2,\dots,z_N)^\top$ は外部から観測することができます.出力信号が複数ありますが,これは多クラス分類とかに使います.例えば,音楽データがJ-POP,クラシック,・・・みたいな複数のカテゴリに属するか一気にチェックできます(クラス分類目的なら誤差関数を交差エントロピーにして,出力をソフトマックスにすべきですが,今日はそのへん割愛します).今回はパラメータ $\boldsymbol{W}_1$, $\boldsymbol{W}_2$, $\boldsymbol{b}_1$, $\boldsymbol{b}_2$ を学習することになります.

多層ニューラルネットワークと深層学習

2層ニューラルネットにさらに層を追加して構造を深くすることできます.昔は深い構造を持つニューラルネットは過学習を起こしたり,誤差逆伝搬法(ニューラルネットの学習アルゴリズム)の勾配が下層部で消失したりして,うまく学習できませんでしたが,最近は色々なテクニックでうまいこと学習できるようになりました.こういった深い構造のネットワークを使った手法を深層学習と言います.深層学習の何がうれしいかは,分かりやすい資料が沢山あるので,そちらを参考にして下さい.ちなみに,この記事の最後に,深層学習の手法の1つである Dropout と多層ニューラルネットを実装したサンプルプログラムを紹介します.

2層ニューラルネットをプログラミングしてみる

では,OCaml と SLAP で2層ニューラルネットをプログラミングしてみましょう.文法が C 言語や MatLab 等とだいぶ違うので,文法を中心に大雑把に説明します.OCaml 自体について,ちゃんと勉強したい人は http://ocaml.jp/ とかを読んで下さい.

準備

まず,utop を立ちあげて,次の「おまじない」で SLAP を読み込んでおきます(行先頭の # はプロンプトなので入力しないでね)

# #use "topfind";;
# #require "slap.top";;
# open Slap.Size;;
# open Slap.Common;;
# open Slap.D;;

各層のユニット数の定義

次に,各層のユニット数 L (= input_dim), M (= hidden_dim), N (= output_dim) を決めます.OCaml では,トップレベルの変数は let 変数名 = 式 で定義するので,

# let input_dim = two;;
# let hidden_dim = five;;
# let output_dim = one;;

とします.ちなみに,twofiveone は「次元」を表す特殊な値で,整数とは違うものです(let n = 42 とかすれば普通の整数を変数に束縛できます).

シグモイド関数を作る

関数は let 関数名 引数1 引数2 ... 引数n = 式 のように定義します.例えば,受け取った値をそのまま返すだけの関数(恒等関数)だと,let f x = x とします.シグモイド関数は

# let sigm a = 1.0 /. (1.0 +. exp (~-. a));;
val sigm : float -> float = <fun>

となります.ここで,+./.~-. はそれぞれ float 用の加算,除算,符号反転を行う演算子です.OCaml では,浮動小数点数と整数を混同しないように,整数用の演算(+, -, *, /, etc.)と浮動小数点数用の演算(+., -., *., /., etc.)を分けています.

関数を呼び出すときは 関数 引数1 引数2 ... 引数n とします.例えば,

# sigm 0.6;;
- : float = 0.645656306225795396

という感じです(式の計算結果を渡したい時は sigm (0.2 -. 0.8) のように括弧で囲みます).

さらに,シグモイド関数をベクトルの各要素に適用する $\mathrm{sigmv}$ 関数を作ってみます.

# let sigmv v = Vec.map sigm v;;
val sigmv : ('a, 'b) vec -> ('a, 'c) vec = <fun>

OCaml では,関数自体を値として扱うことができ,関数を関数に渡したり,関数を返り値にしたりできます.Vec.map は受け取った関数(sigm)をベクトル(v)の各要素に適用する関数です.

ベクトルと行列の計算

試しに,ベクトルを作ってみましょう.

# let x = Vec.make1 input_dim;;
val x : (z s s, 'a) vec = R1 R2
                           1  1

Vec.make1 という関数に three という引数を渡して,結果を変数 x に束縛しています.関数 Vec.make1 は受け取った「次元」(ここでは input_dim = two)の長さを持つベクトルを作り,1で初期化して返します.two はその名の通り,「2次元」を表しているので,2次元ベクトルが帰ってきます.

さらに,重み行列とバイアスも作ります(注: OCaml では変数名・関数名の先頭に大文字は使えません)

# let w1 = Mat.random hidden_dim input_dim;;
# let b1 = Vec.random hidden_dim;;
# let w2 = Mat.random output_dim hidden_dim;;
# let b2 = Vec.random output_dim;;

隠れ層の出力信号 $\boldsymbol{y} = \mathrm{sigmv}(\boldsymbol{W}_1\boldsymbol{x}+\boldsymbol{b}_1)$ を計算するために,まずは $\boldsymbol{W}_1\boldsymbol{x}+\boldsymbol{b}_1$ を計算する必要があります.この計算は gemv という関数を使って

# let u = gemv ~trans:normal w1 x ~beta:1.0 ~y:(Vec.copy b1);;
val u : (z s s s s s, 'a) vec =
      R1        R2       R3        R4      R5
0.182888 -0.358936 -1.68752 -0.948303 1.05315

と書けます.gemv ~trans:normal ~alpha a x ~beta ~y で $\boldsymbol{y} \gets \alpha\boldsymbol{A}\boldsymbol{x} + \beta\boldsymbol{y}$ を計算しますが,~y:b1 と書いてしまうと b1 が上書きされてしまうので,ここではコピーしてから渡しています(~alpha, ~beta, ~y は省略可能で,省略時は ~alpha:1.0, ~beta:0.0~y は新しく確保されたベクトルになる).ちなみに,gemv ~trans:trans ~alpha a x ~beta ~y だと $\boldsymbol{y} \gets \alpha\boldsymbol{A}^\top\boldsymbol{x} + \beta\boldsymbol{y}$ を計算します.あとは,この usigmv に渡せば,$\boldsymbol{y}$ が計算できます.

ここまで紹介してきたパーツを組み合わせると,2層ニューラルネットの出力を計算する関数 exec を作れます.

# let exec w1 w2 b1 b2 x =
    let y = sigmv (gemv ~trans:normal w1 x ~beta:1.0 ~y:(Vec.copy b1)) in
    let z = sigmv (gemv ~trans:normal w2 y ~beta:1.0 ~y:(Vec.copy b2)) in
    z;;
val exec :
  ('a, 'b, 'c) mat ->
  ('d, 'a, 'e) mat ->
  ('a, 'f) vec -> ('d, 'g) vec -> ('b, 'h) vec -> ('d, 'i) vec = <fun>

ちなみに,let 変数名 = 式 in ... はローカル変数を作る構文です.exec w1 w2 b1 b2 x とすると,ニューラルネットの出力が返ってきます.

2層ニューラルネットの教師あり学習

さっきは,重み行列とバイアスを乱数で与えていましたが,これでは意味がありません.あくまで,ベクトルの分類に使いたいわけですから,ランダムではなく「いい感じの」重み行列とバイアスが欲しいですよね?ここからが,機械学習の本領発揮です.特徴ベクトル $\boldsymbol{x}$ と(モデル出力 $z$ の)目標値 $t$ = 1 or 0(あるクラスに属するか,属さないか)の組を訓練事例として大量に与えて,そこから訓練事例をうまく分類するような $\boldsymbol{W}_1$, $\boldsymbol{W}_2$, $\boldsymbol{b}_1$, $\boldsymbol{b}_2$ を見つけます.この計算は誤差逆伝搬法で行うことができます.アルゴリズムの導出については 村上・泉田研究室 ニューラルネットワーク 第6章 誤差逆伝搬法 が分かりやすいと思うので,興味のある人は読んでみるといいかもしれません.2層の場合について,簡単にアルゴリズムをまとめると,

  1. $\boldsymbol{W}_1$, $\boldsymbol{W}_2$, $\boldsymbol{b}_1$, $\boldsymbol{b}_2$ を乱数で初期化する
  2. 収束するまで,以下のパラメータの更新を繰り返す
    1. 訓練事例の集合から組 $(\boldsymbol{x}, \boldsymbol{t})$ を取ってくる(復元抽出)
    2. $\boldsymbol{y} \gets \mathrm{sigmv}(\boldsymbol{W}_1\boldsymbol{x}+\boldsymbol{b}_1)$
    3. $\boldsymbol{z} \gets \mathrm{sigmv}(\boldsymbol{W}_2\boldsymbol{y}+\boldsymbol{b}_2)$
    4. $\boldsymbol\delta_2 \gets (\boldsymbol{z}-\boldsymbol{t}) \otimes \boldsymbol{z} \otimes (\boldsymbol{1}-\boldsymbol{z})$
    5. $\boldsymbol\delta_1 \gets \boldsymbol{W}_2^\top \boldsymbol\delta_2 \otimes \boldsymbol{y} \otimes (\boldsymbol{1}-\boldsymbol{y})$
    6. $\boldsymbol{W}_1 \gets \boldsymbol{W}_1 - \eta \boldsymbol\delta_1 \boldsymbol{x}^\top$
    7. $\boldsymbol{W}_2 \gets \boldsymbol{W}_2 - \eta \boldsymbol\delta_2 \boldsymbol{y}^\top$
    8. $\boldsymbol{b}_1 \gets \boldsymbol{b}_1 - \eta \boldsymbol\delta_1$
    9. $\boldsymbol{b}_2 \gets \boldsymbol{b}_2 - \eta \boldsymbol\delta_2$

ただし,$\eta > 0$ は学習率,$\otimes$ はベクトルの要素ごとの乗算を表します.

誤差逆伝搬法でパラメータを1ステップ更新する関数を作る

まず,4. と 5. の式に出てくる $\boldsymbol{z} \otimes (\boldsymbol{1}-\boldsymbol{z})$ と $\boldsymbol{y} \otimes (\boldsymbol{1}-\boldsymbol{y})$ を計算する関数を作ります.この式はシグモイド関数の微分の計算に対応しています.

# let sigmdv v =
    let ones = Vec.make1 (Vec.dim v) in
    Vec.mul v (Vec.sub ones v);;

これを使って,パラメータ更新を行う関数 train を作ります.

# let train eta w1 w2 b1 b2 x t =
    let y = sigmv (gemv ~trans:normal w1 x ~beta:1.0 ~y:(Vec.copy b1)) in
    let z = sigmv (gemv ~trans:normal w2 y ~beta:1.0 ~y:(Vec.copy b2)) in
    let delta2 = Vec.mul (Vec.sub z t) (sigmdv z) in
    let delta1 = Vec.mul (gemv ~trans:trans w2 delta2) (sigmdv y) in
    ignore (ger ~alpha:(~-. eta) delta2 y w2);
    ignore (ger ~alpha:(~-. eta) delta1 x w1);
    axpy ~alpha:(~-. eta) b2 ~x:delta2;
    axpy ~alpha:(~-. eta) b1 ~x:delta1;;

基本的に,さっきのアルゴリズム(箇条書き)と行単位で対応しているので,何をしているかは解ると思いますが,少しクセのある関数だけ紹介します.ger は $\boldsymbol{x}\boldsymbol{y}^\top$ のような outer product を計算する関数で,ger ~alpha x y a で $\boldsymbol{A} \gets \boldsymbol{x}\boldsymbol{y}^\top + \alpha\boldsymbol{A}$ を計算します.axpy は積和演算の関数で,axpy ~alpha ~x y で $\boldsymbol{y} \gets \alpha\boldsymbol{x} + \boldsymbol{y}$ を計算します.
あと,ignore は返り値を捨てる関数で,返り値が不要なときは明示的にこの関数を書くことが多いです.

訓練させてみる

ここでは,訓練集合として XOR pattern を作ってみます.

# let samples = [|
      (Vec.of_array_dyn input_dim [|1.0; 1.0|], Vec.make output_dim 0.0);
      (Vec.of_array_dyn input_dim [|1.0; 0.0|], Vec.make output_dim 1.0);
      (Vec.of_array_dyn input_dim [|0.0; 1.0|], Vec.make output_dim 1.0);
      (Vec.of_array_dyn input_dim [|0.0; 0.0|], Vec.make output_dim 0.0);
   |];;

これは入力ベクトルと出力ベクトルの組を要素とする配列で,要素を取り出すときは,samples.(i) とします(i は添字).

# let (x, t) = samples.(0);;

これを train に渡せば,パラメータを1ステップだけ更新してくれます.ただし,たった1ステップの更新で全ての訓練事例を分類できるようになるわけではないので,とりあえず,5000回くらい繰り返してみます.

# for i = 1 to 5000 do
    let (x, t) = samples.(i mod (Array.length samples)) in
    train 0.5 w1 w2 b1 b2 x t
  done;;

では,いい感じのパラメータが求まったか確認してみます.

# let (x, t) = samples.(0);;
val x : (z s s, '_a) vec = R1 R2
                            1 1
val t : (z s, '_a) vec = R1
                          0
# exec w1 w2 b1 b2 x;;
- : (z s, 'a) vec =        R1
                    0.0742082

まだ少し誤差が残ってはいますが,だいたい目標値 0 に近い値(0.074)が出ています.クラス分類の用途では,0.5 を閾値として判定するので,このくらいの値でも十分です.ぜひ,他のデータに対してもうまく分類できるか調べてみて下さい.もし,うまくいかないようであれば,再度重みの更新式を反復します.

ここまでに作ったプログラムに gradient checking を加えたサンプルプログラムが https://github.com/akabe/slap/blob/master/examples/neural-network/two-layer/two_layer_neural_network.ml に置いてあります.gradient checking は誤差逆伝搬法で求めた勾配が正しいかどうか,愚直な数値微分を使って確認する方法です(非常に遅いのでデバックの時だけ使います).Gradient checking については Gradient checkingは超大事 が解りやすかも.

静的サイズ検査

では,SLAP の静的サイズ検査により,次元の合わない計算がコンパイル時に検査されるか確認してみましょう.対話環境だとコンパイルと実行の違いが解りにくいので,先ほど紹介したサンプルプログラム two_layer_neural_network.ml を使ってみましょう.試しに,コンパイルして実行してみると,

$ git clone https://github.com/akabe/slap.git
$ cd slap/examples/naural-network/two-layer/
$ ocamlfind ocamlopt -linkpkg -package slap -short-paths two_layer_neural_network.ml
$ ./a.out
w1 = [ 0.0326249 0.594539 
         4.67399  4.62221 
         5.39787 -3.44113 
        -0.77698  3.17245 
         2.27682 -4.23508 ]
w2 = [ -0.314653 6.01105 -6.23296 -2.79963 4.70526 ]
b1 = [ -0.693542 -1.07445 1.23425 -0.234505 -0.811487 ]
b2 = [ 0.0273244 ]
x = [ 1 1 ]; t = [ 0 ]; z = [ 0.08827 ]
x = [ 1 0 ]; t = [ 1 ]; z = [ 0.932276 ]
x = [ 0 1 ]; t = [ 1 ]; z = [ 0.922299 ]
x = [ 0 0 ]; t = [ 0 ]; z = [ 0.0405274 ]

というように求まった重みと分類結果を表示します.うまくいってますね!
このtwo_layer_neural_network.ml の120行目をわざと変に書き換えてみます.元のコードは

120:    train 0.5 w1 w2 b1 b2 x t

ですが,b2b1 に書き換えてみます.

120:    train 0.5 w1 w2 b1 b1 x t

b1 は5次元ベクトルで b2 は1次元ベクトルなので,このように書き換えると,次元が合わなくなるはずです.これを保存して,もう一度同じようにコンパイルしてみると...

$ ocamlfind ocamlopt -linkpkg -package slap -short-paths two_layer_neural_network.ml
File "two_layer_neural_network.ml", line 120, characters 17-19:
Error: This expression has type
         (z s s s s s, 'a) vec = (z s s s s s, float, rprec, 'a) Slap.Vec.t
       but an expression was expected of type
         (z s, 'b) vec = (z s, float, rprec, 'b) Slap.Vec.t
       Type z s s s s is not compatible with type z

素晴らしい!次元が合わなくなったので,コンパイルエラー(型エラー) が起きました!エラーメッセージの z s s などは次元を表していて,s の個数が次元に対応します.上のエラーメッセージは「1次元ベクトル((z s, 'b) vec)を渡すはずが,5次元ベクトル((z s s s s s, 'a) vec)が渡されています.次元が合いません.」という意味です.こんな感じで,SLAP を使うと次元のあわない計算を実行せずに検出できます.エラーメッセージを読めるようになるまで時間がかかるかもしれませんが,それでも,静的サイズ検査でバグを事前に検知できるのは大きなメリットだと思います.

多層ニューラルネットワークと Dropout の実装

今回は簡単のために2層で固定にしましたが,もう少し工夫すると,層の数も変更可能な多層ニューラルネットも実装できます(少し難しいので,この記事では扱いません).この多層ニューラルネットに深層学習のアルゴリズムの1つである Dropout [Hinton, 2012] を実装したサンプルプログラムが https://github.com/akabe/slap/tree/master/examples/neural-network/multilayer にアップしてあります.Dropout は過学習を防ぐ方法の1つで,ユニットをランダムに無視しながら学習することで,バギングのような効果を生み出すらしいです.コンパイルは

$ ocamlfind ocamlopt -linkpkg -package slap -short-paths neuralNetwork.mli neuralNetwork.ml multilayer_neural_network.ml

で出来ます.

今回はSLAPを使って,とりあえず何か作ってみようという事でニューラルネットを作ってみました.この記事ではSLAPの機能の一部しか紹介していませんが,そのうち他の機能についても記事を書くかもしれません.質問・ツッコミは歓迎します.

この投稿は ML Advent Calendar 201418日目の記事です。