Haskellでいってみよう

日曜プログラマにも満たないレベルでもHaskellで何かソフトウェアを作りたい!

レイトレーシング(11): 画質向上策2点

今回は寄り道

前回から半年も経ってしまった orz

前回で完全拡散面の反射まで実装した。当然次は鏡面反射、と思われたら 申し訳ないところだが、反射(と屈折)は手間がかかるので後回し。 今回は寄り道として生成画像の画質改善に取り組む。

1. 直接光は古典レイトレーシング

フォトンマッピングでは、放射するフォトン数に限界があり、多数のフォトン から輝度推定しても"まだら模様"のようになってしまう。また、点光源では本来 影の縁はくっきりするものだが、残念ながらぼやけてしまう。

例の本(フォトンマッピング、Herik Wann Jensen、オーム社)では、 実践的なアルゴリズムとして「直接光はフォトンマップから推定するのではなく 古典レイトレーシングで」と書いてある(ようだ)。

確かにランダムなフォトンから推定するより数式できっちり計算できる 古典レイトレーシングを用いたほうが画質的には有利だろう。

前回までで、フォトンマッピング法の画像と比較するため古典レイトレーシングも 実装してきた。これをそのまま流用すればよい。基本的な考え方は次の通り。

実装

ではまずフォトン追跡の実装から。

tracePhoton :: [Object] -> Int -> Photon -> IO [PhotonCache]
tracePhoton os l (wl, r) = do
  let is = calcIntersection r os
  if is == Nothing
    then return []
    else do
      let (p, n, m) = fromJust is
      i <- russianRoulette wl [reflectance m]
      pcs <- if i > 0
        then reflect p n os l wl
        else return []
      if (useClassicForDirect == False || l > 0) && diffuseness m > 0.0
        then return $ ((wl, initRay p (getDir r)) : pcs)
        else return pcs

tracePhotonに$l$という第二引数を設けている。これは追跡の深さを 表すもので、物体に当たって反射するたびに1ずつ増やして反射方向の追跡を するようにした。つまり$l=0$ならまだ反射していないから「光源からの 直接光」である。

  if (useClassicForDirect == False || l > 0) && diffuseness m > 0.0

これがその判定をしているところ。

  • useClassicForDirectは、光源からの直接光を古典レイトレーシングで計算 するかどうかのフラグ。Falseはこれまでどおりフォトンマップから推定する という意味。
  • l > 0は前述の通り一回以上反射したフォトンを意味する。
  • diffuseness m > 0.0フォトンが衝突した物体表面が、少しでも 拡散反射の要素を持っている=フォトンを記録する必要があるという意味。

これらのいずれにも合致しない場合は、そのフォトンの衝突点ではフォトンマップに 記録しない。

次にレイトレーシング部分を見てみよう。

traceRay :: Int -> Double -> KT.KdTree Double PhotonInfo -> [Object]
         -> [Light] -> Ray -> IO Radiance
traceRay 10 _ _ _ _ _ = return radiance0
traceRay l pw pmap objs lgts r
  | is == Nothing = return radiance0
  | otherwise     = return (em + di + ii)
  where
    is = calcIntersection r objs
    (p, n, m) = fromJust is
    em = sr_half *> emittance m
    di = if useClassicForDirect
      then brdf m $ foldl (+) radiance0 $ map (getRadianceFromLight objs p n) lg ts
      else radiance0
    ii = estimateRadiance pw pmap (p, n, m)

where節の後半でuseClassicForDirectというフラグをチェックしている。 ここの仕組みはこうだ。

  • useClassicForDirectが真の場合は、物体上の点に注ぐ直接光を古典 レイトレーシングで計算する。getRadianceFromLight関数がそれ。光源の 数だけ繰り返している。偽の場合は、古典では輝度を計算しないので値ゼロ (=radiance0)としている。
  • その点には間接光も入ってくるので、その輝度は従来通りフォトンマップから 推定する。(ifの下のestimateRadiance)

なおフォトンマップから求められる輝度は、古典レイトレーシングを使う場合には 「直接光のフォトンを記録しない」ようにしてあるので、ちゃんと間接的に注ぐ 光だけで輝度推定される。古典を使わない場合は直接光のフォトンも記録されて いるので従来と同じ結果となる。

画像比較

では生成された画像を比べてみる。まずは点光源の場合。

f:id:eijian:20180108095118p:plain

効果は一目瞭然だろう。直接光を古典で計算した場合(図2)は奥の壁、床、球の 上面が滲みなく滑らかに描画されている。直接光の当たらない天井や影の部分は 間接光が支配的なので致し方ないが。注目は球の影。従来はどうしても縁がぼやけて しまったがそこがくっきりしている。点光源はこうでないと。次に平面光源の場合。

f:id:eijian:20180108095144p:plain

こちらも同様にスムーズだ。最後に平行光源の場合。

f:id:eijian:20180108095153p:plain

図5は壁面の直接光の当たる境界がぼやけてしまっているが、古典を使うと くっきりしたエッジになっている。本来こうあるべきだ。 (こちらの画像の場合天窓からの光以外は全て間接光なので、ほとんどマダラが 残ったままだが、ここは我慢)

処理時間も比べてみよう。直接光を古典レイトレーシングに任せることで 直接光のフォトンを記録しないと書いた。つまりフォトンマップが小さく なるということだ。結果は次の表のとおり。

画像 直接光 フォトン 処理時間(MM:SS)
1 点光源 フォトンマップ 318,164 6:50
2 点光源 古典レイトレ 117,852 1:22
3 面光源 フォトンマップ 288,976 6:14
4 面光源 古典レイトレ 89,420 1:09
5 平行光源 フォトンマップ 321,142 23:07
6 平行光源 古典レイトレ 120,841 2:01

どうやらフォトンマップが小さいと輝度推定時にフォトン収集の時間が短くて 済むようで、圧倒的に処理時間が短くなることがわかった! 特に平行光源の例では1/10に短縮できている。なぜこれほど短縮できるのか (もしくはフォトンマップだけだと時間がかかるのか)はよくわからないが。

2. えせアンチエリアシング

もう一つの画質改善ネタであるアンチエリアシングに移ろう。 アンチエリアシングの詳しい説明はWebを漁れば山ほど出てくるのでそちらを 見ていただくとして、今回実装した仕組みを説明しよう。

仕組み

レイトレーシングは画面を格子状に分割し、その分割された小さい四角 (ピクセル)の中心に向かって光を逆向きに追跡する(追跡レイとする)。 各ピクセルの輝度値がある閾値を超えたら、そのピクセルの内側のいくつかの 点を追加で追跡し、その輝度を平均する。

節タイトルで「えせ」としたのは、ちゃんとアンチエリアシングしたければ 追加する追跡レイはピクセル内の「ランダム」な場所に対して行うべきだから である。規則的に追跡レイを追加するとどうしてもエリアシングの問題が 残るらしい。しかし今回は、小難しいことを考えるのがだるかったので実装の 簡便さを優先した。。。時間があったらちゃんと実装しよう。

まず、すべてのピクセルで追跡レイを追加していたら時間ばかりかかって非効率 なので追加すべきかどうかをまず判断する。そのため、とあるピクセル(p)とその 周りの8つのピクセルについてそれぞれ差を求め、閾値を超えているか判定する。 1点でも閾値を超えていたら追加レイを飛ばす。(下図)

f:id:eijian:20180108095205p:plain

なお、このプログラムでは差を求める際に輝度を比較するのではなくRGB値、 すなわち「画面に出力する色」に変換してから比較している。なぜなら、 結局目に見えるところでの色の違いがある場合に追加レイを飛ばしたいから。 最初は輝度で試したが、閾値の設定がシーンの明るさ・光の量で変える必要が あり設定が難しく、実用的でないと判断した。

次に追加レイを飛ばす場所だが、上記の通り規則的にした。ピクセル内の 中心以外の4点だ。単純にピクセルを4分割してそれぞれの中心点に向かって 追跡している。

f:id:eijian:20180108095212p:plain

これら4点と最初に計算した中心点の全5点から得られた輝度を平均するのだが、 その方法には単純平均するか重み付けするかがある。結果に微妙な差が出るが 大きな問題はないと判断し、今回は単純平均する。ランダムに追加レイを飛ばす 場合は中心からの距離に差が出ると思うのでちゃんと重み付けした方が よいだろう。

実装

アンチエリアシングを実施するかどうかはメインループの最後のところ。

  :
  forM_ [0..(V.length pixels - 1)] $ \i -> do
    rgb <- smooth antiAliasing tracer pixels i
    putStrLn $ rgbToString rgb

smooth関数で処理している。引数のantiAliasingアンチエリアシングを 実施するかどうかのフラグ。smoothは次の通り。

smooth :: Bool -> (Ray -> IO Radiance) -> V.Vector Rgb -> Int -> IO Rgb
smooth False _ ims i = return (ims V.! i)
smooth True tracer ims i
  | isDifferent i ims imgOffset == False = return (ims V.! i)
  | otherwise = do
    l <- retrace tracer i
    return $ avg ((ims V.! i):l)

まずフラグが偽なら何もせず最初に求めた色を返す。真の場合、まずは追跡レイを 追加するかどうか判断する(isDifferent)。追加が必要なければやはりそのまま 色を返す。追加が必要となったら再トレースだ(retrace)。そこで得られた色を avgで平均して返せばOK。

追加が必要かどうかの判断は次の通り。

isDifferent :: Int -> V.Vector Rgb -> [Int] ->  Bool
isDifferent _ _ [] = False
isDifferent p rs (i:is)
  | p' <  0         = isDifferent p rs is
  | p' >= length rs = isDifferent p rs is
  | df              = True
  | otherwise       = isDifferent p rs is
  where
    p' = p + i
    df = diffRgb (rs V.! p) (rs V.! p')

少し説明が必要かもしれない。2番目の引数は、すでに 計算済みの画面全部の色値(RGB)の配列(Vector)である。着目しているピクセルは 第一引数pに配列内の位置として与えられる。p'が比較対象となる 周りのピクセル(の位置)は、第3引数で与えられている。 これは上の方で次のように定義されている。

imgOffset :: [Int]
imgOffset = [-xres-1, -xres, -xres+1, -1, 1, xres-1, xres, xres+1]

画面の横方向画素数xresを使い、着目点pに対するオフセットを 示しているわけだ。このようにして、周りの8個の点と比べ、閾値を 超えたら直ちにTrueを返すようにしてある。周りの点が画面をはみ出す 場合はもちろん無視だ(条件 p' < 0p' >= length rs)。

次に再トレース部分を説明しよう。

retrace :: (Ray -> IO Radiance) -> Int -> IO [Rgb]
retrace tracer p = do
  let p' = (fromIntegral (p `div` xres), fromIntegral (p `mod` xres))
  ls <- mapM tracer $ map generateRay' (map (badd p') blur)
  return $ map radianceToRgb ls

といっても大したことはない。現在の着目点pに対し、1/4ずつずれた4点を 求め、それぞれ普通に追跡して結果をリストにして返しているだけだ。 ずれた点を計算しやすくするため、blurをあらかじめ定義している。

blur :: [(Double, Double)]
blur = [(-0.25, -0.25), (-0.25,  0.25), ( 0.25, -0.25), ( 0.25,  0.25)]

生成画像の比較

では実際にアンチアリアスを効かせた画像を生成してみよう。従来の画像との 比較はつぎのとおり。

f:id:eijian:20180108095220p:plain

f:id:eijian:20180108095229p:plain

ぱっと見はよくわからないかもしれないが、特に球の上部や平面光源の縁が気持ち なだらかになったのがわかるだろうか?これがアンチエリアシングの効果だ。 ただ、いくつかのピクセルについて追加レイを飛ばすので、当然計算時間が 余計にかかる。

アンチエリアシング フォトン 処理時間(MM:SS)
89420 1:09
89420 1:52

画面の複雑さや色合いにもよるので一般にどのぐらい時間が増えるかは なんとも言えないが、この例ではほぼ2倍の計算時間がかかったことになる。 この辺は画質とのトレードオフだが、世間のCGソフトは普通アンチエリアシングが 効いているからなぁ。

まとめ

今回は多少なりとも画質を向上する策として、以下の2点の改善を実施した。

両方とも、簡便な方法を採った割にはそれなりにうまくいったと思う。 また、古典レイトレーシングを使うことで大幅な計算時間の短縮ができた ことは収穫だった。

次回はとうとう(?)鏡面反射、屈折に対応しよう。うまく行けばいわゆる 「集光模様」が得られるはずだ。


ここまでのソースはこちら

レイトレーシング(10): 拡散反射と面光源

DeepLearning回からだいぶ間があき、さらにレイトレーシングは1年半以上も 時間が空いてしまった。それは拡散反射でとある処理の解法で壁にぶつかって いたからだ。それが非常に簡単に解決できるとわかったのでちょっと進める ことにした。

0. 前回まで

前回までで、フォトンマッピング法を用いたレイトレーシングプログラムを 作成した。といってもフォトン追跡は何か物体にぶつかったらそれでおしまい、 反射は未対応という中途半端な状態だ。

久しぶりなので気分を一新して新しいシーンに変えてみよう。 実はこのページのパクリだ。 ちなみにこの連載は知りたいことが多く書いてありとても参考にさせて いただいている。売り物の3DCGソフトを使っての説明だが理論的なところが 詳しくよくわかる。

さて、前回までで作ったフォトンマッピングプログラムと古典レイトレーシングで 描画したものがこちら。

f:id:eijian:20170624204748p:plain

一目瞭然、フォトンマッピングのほうは影がぼやけてしまう。 点光源なので影は古典のようにくっきりとした縁、かつ真っ黒でなければ ならない。前の記事でも書いたがフォトンマッピングはある点の周囲のフォトンを 集計して輝度を計算するので、影のなかにフォトンがないと収集する範囲を 広げてしまうのだ。

ということでさらの前回のおさらいだが、フォトンを収集する範囲を 限定して描画してみる。設定は0.2 [m]だ。それより外のフォトンは無視して 輝度を判定する。

f:id:eijian:20170624204758p:plain

残念ながら影の縁はぼやけてしまうが、上の画像よりははるかに本物っぽい。 今後はひとまずフォトンの収集範囲を0.2 [m]として、色々試すことにしよう。 ちなみにフォトン数は20万個、これでも残念ながら"まだら模様"になるのは 今は我慢しよう。

なおフォトンマッピングと古典レイトレーシングでは輝度の計算方法が全く 異なるわけだが、両者の明るさや色の加減はほぼ同等であることに注目したい。 レイトレーシングといっても要するに物理(光学)シミュレーションなので 同じ設定なら(ほぼ)同じ結果が出てこないといけないのだ。その点で、 ひとまず両者のプログラムに大きな間違いはなさそうだ。

1. 拡散反射の実装

さて本題の拡散反射である。現時点ではシーン中の全物体は完全拡散反射面 としているので、フォトンを追跡して物体にぶつかったらフォトンマップに 記録し、別方向へ反射させればよい。ところが長い間詰まっていたのは その反射方向をどうするか、であった。

(1) 拡散反射ベクトルの生成法

完全拡散反射面ではフォトンが物体にぶつかるとその面の半球状のいずれかに 等確率で反射する。当初は、まず上方向の平面を考え、半球のどこかへ向かう ベクトルを求めてからそれを法線方向によって回転させようとした。

f:id:eijian:20170624204806p:plain

反射ベクトルは点光源でフォトンをランダムに生成した経験から、それを 上半分だけに限定すれば生成できる。あとは回転だが、回転行列を作るのか四元数を 使うのか、よく分からないまま止まっていた。一般的にどうしているのか検索も してみたがどうもそれらしい情報を見つけられない。

で、先日極めて容易に生成できる方法に気がついた。球状のランダムなベクトルを 生成したあと、求めたい面の法線とのなす角が90度以下ならそのまま採用、90度 以上なら反転させればいいのだ。

f:id:eijian:20170624204816p:plain

なんと簡単なことに長時間悩んだことか、恥ずかしくなるぐらいだ。。。 この処理をコードにすると次の通り。generateRandomDir?は以前作って いたもの(生成方法により四種類あり、?はその番号、一例として四番目を示す)。

generateRandomDir4 :: IO Direction3
generateRandomDir4 = do
  y' <- MT.randomIO :: IO Double
  p' <- MT.randomIO :: IO Double
  let y = y' * 2.0 - 1.0
      p = p' * 2.0 * pi
      r = sqrt (1 - y * y)
      x = r * (cos p)
      z = r * (sin p)
      v = initPos x y z
      len = norm v
  return $ fromJust $ normalize v

diffuseReflectionが反射ベクトルの生成だ。

diffuseReflection :: Direction3 -> IO Direction3
diffuseReflection n = do
  d <- generateRandomDir4
  let cos = n <.> d
  return $ if cos > 0.0 then d else negate d

こんなたった数行でできることに悩まされていたとは・・・。

(2) フォトン追跡をいじる

反射ベクトルを作ることができたので、フォトン追跡にも手を 入れよう。従来の追跡処理は次のようであった。

  • フォトンの進行方向にある物体についてその距離を取得する
  • フォトンより前にあるものの中から最も近いものを選びPhotonCacheとして返す

コードはこちら。

tracePhoton :: [Object] -> Photon -> IO [PhotonCache]
tracePhoton os (wl, r) = do
  let iss = filter (\x -> fst x > nearly0) (concat $ map (calcDistance r) os)
      (t, s) = head $ sortBy (comparing fst) iss
  return [(wl, initRay (target t r) (getDir r))]

ここは大幅に変更している。やりたいことは次の通り。

  • フォトンがぶつかる最初の物体との交点を得る(calcIntersection)
  • ぶつからなければ空リストを返す。
  • ロシアンルーレット法を使い物体の反射率以下なら反射(reflect)し さらにフォトン追跡をしてPhotonCache情報を得る。
  • 反射率以上ならフォトンが吸収されたとして追跡を止める。
  • 交点が拡散反射面(*1)ならその点のPhotonCahce情報を作成する。
  • 交点および反射方向から全PhotonCache情報を集めて返す。

(*1)について、今は拡散反射面という前提だが、将来の拡張のため このような判定を入れている。上記処理のコードは次の通り。

tracePhoton :: [Object] -> Int -> Photon -> IO [PhotonCache]
tracePhoton os l (wl, r) = do
  let is = calcIntersection r os
  if is == Nothing
    then return []
    else do
      let (p, n, m) = fromJust is
      i <- russianRoulette wl [reflectance m]
      pcs <- if i > 0
        then reflect p n os l wl
        else return []
      if diffuseness m > 0.0
        then return $ ((wl, initRay p (getDir r)) : pcs)
        else return pcs

reflect :: Position3 -> Direction3 -> [Object] -> Int -> Wavelength
        -> IO [PhotonCache]
reflect p n os l wl = do
  dr <- diffuseReflection n
  let r' = initRay p dr
  pcs <- tracePhoton os (l+1) (wl, r')
  return pcs

なお、後々の拡張のため、反射回数を引数に追加している(l)。 ちなみにrussianRouleteはこのようになっている。

russianRoulette :: Wavelength -> [Color] -> IO Int
russianRoulette wl cs = do
  r <- MT.randomIO :: IO Double
  return $ rr wl cs 0.0 r (length cs)

rr :: Wavelength -> [Color] -> Double -> Double -> Int -> Int
rr _ [] _ _ len = len
rr wl (c:cs) c0 r len
  | r < c'    = len
  | otherwise = rr wl cs c' r (len - 1)
  where
    c' = c0 + selWl wl c
    selWl :: Wavelength -> Color -> Double
    selWl Red   (Color r _ _) = r
    selWl Green (Color _ g _) = g
    selWl Blue  (Color _ _ b) = b

まず[0,1]の一様乱数rを作り、Colorのリストから指定の波長の値を取り出して そのrを上回ったところでリスト要素の順番(逆順)を返す。 分かりにくいので例を示す。物質の反射率、透過率が次のようで あったとする。ただし各色の反射率+透過率は1.0を超えない。

反射率: 赤=0.5、緑=0.1、青=0.6

透過率: 赤=0.2、緑=0.6、青=0.4

1.0-(反射率+透過率)は吸収率だ。ここで赤を例に判定しよう。その場合、 リスト[Color]は[[0.2,0.6,0.4],[0.5,0.1,0.6]]である。赤について取り出すと [0.2,0.5]だ。rが0.6であれば、

  • 1番目の透過率=0.2と比べるとrが大きいので次へ
  • 2番目の反射率=0.5だが先の0.2とあわせて0.7としてrと比較すると、 今度はrが小さいので、この時の要素数1を返す

ということで1になる。

f:id:eijian:20170624204858p:plain

もしrが0.1なら1番目の透過率0.2より小さいので結果は2なわけだ。この番号は、 先のフォトン追跡の際に吸収、反射、透過のどれかを表すのに使っており、 このプログラムでは0=吸収、1=反射、2=透過と決めた。なので、 russianRouletteの返値が0より大きければフォトンを反射してさらに追跡 するようにしている。

      i <- russianRoulette wl [reflectance m]
      pcs <- if i > 0
        :

さて、それでは次に画像を生成してみよう。

(3) 間接光の効果はすごい

拡散反射を実装したので、あらためてフォトンマップを生成する。

$ cabal build
  :
$ dist/build/pm/pm > scene1.map
$ dist/build/rt/rt < scene1.map > scene1.ppm

以下に、機能拡張したプログラムで生成した画像と比較のため直接光のみの 画像および古典版を示す。古典版はいわゆる「環境光」に2 [mW/m2/sr]の輝度を与えた。

f:id:eijian:20170624204916p:plain

結果は明らかだ。まだら模様と影の縁がぼやけているのは仕方ないとして、

  • 球の下の方は床からの照り返しでほんのり明るくなっている
  • 球の左右はそれぞれ赤・青の壁からの反射で少し赤み・青みがかっている
  • 天井も壁からの照り返しがある

など、だいぶリアルな画像が得られたと思う。まだら模様もまあいわゆる 「味わい」と言えなくもない(と自分を慰めてみる)。

2. 面光源も

さて、拡散反射による間接光が思いの外いい効果を出したので、もう少し 拡張してみる。先の例では点光源なのに影の縁がぼやけていた。でも大きさの ある光源ならぼやけていて当たり前、違和感がないはず?

まずは実装が簡単そうな平面光源を追加してみる。

(1) 平面光源

点光源は一点から球状のあらゆる方向へフォトンが放射されるが、面光源は 面のそれぞれの点から半球状に放射される。そう、拡散反射の実装で拡散反射 ベクトルを求めたがそれがそのまま使える。あとは平面をどう定義するかだが、 あまり手をかけたくないので2つのベクトル($\boldsymbol{d_1}, \boldsymbol{d_2}$)で表される 「平行四辺形」にしたいと思う(下図)。

f:id:eijian:20170624204933p:plain

フォトンの放出方法だが、まず面光源のどこから放出されるかを2つの一様 乱数($0 \leq u \leq 1, 0 \leq v \leq 1$)で決める。放出点を $p$ とすると、

$ \boldsymbol{p} = \boldsymbol{p_0} + u \cdot \boldsymbol{d_1} + v \cdot \boldsymbol{d_2} $

となる。あとは拡散反射と同様に半球状のランダムな方向ベクトルを生成して フォトンを飛ばしてやれば良い。コードにするとこんな感じ。

data Light =
  PointLight
    { lcolor :: Color
    , lflux  :: Flux
    , pos    :: Position3
    }
  | ParallelogramLight
    { lcolor :: Color
    , lflux  :: Flux
    , pos    :: Position3
    , nvec   ::   Direction3
    , dir1   :: Direction3
    , dir2   :: Direction3
    }

generatePhoton (ParallelogramLight c _ p n d1 d2) = do
  wl <- MT.randomIO :: IO Double
  t1 <- MT.randomIO :: IO Double
  t2 <- MT.randomIO :: IO Double
  d  <- diffuseReflection n
  let r = initRay (p + t1 *> d1 + t2 *> d2) d
      w = decideWavelength c wl
  return (w, r)

生成してみたフォトンマップがこちら。天井に黒い四角が見える。これが面光源。

f:id:eijian:20170624204946p:plain

さて、これを使ってレイトレーシングしたいのだが、問題は大きさのある光源の 場合「目に見える」ということだ。点光源は大きさがないので無視していたが、 今度はそうはいかない。なのでレイトレーシング時に他の物体と同じように 交差判定とかしてやらないといけない。というか本来は逆だ。物体がたまたま 光を放出しているから光源なのだ。というわけで、物体の表面属性に「発光」を 加えてやろう。他にも今後対応する(かもしれない)反射・透過関係も入れて、 次のように定義した。

data Material = Material
  { reflectance   :: Color
  , transmittance :: Color
  , specularRefl  :: Color      -- specular reflectance
  , emittance     :: Radiance
  , ior           :: Color      -- index of refraction
  , diffuseness   :: Double     -- 拡散性
  , smoothness    :: Double
  } deriving Eq

このemittanceがそれだ。ここまでくると、物体定義を使って 光源としての処理をしてやるのが筋ってもんだ、と思う。つまり Light型を物体と別に定義するのではなく物体を使ってフォトン放出を するということだ。が、ちょっと手間なので今回は割愛、今後のがんばりに期待。

で、代わりに「平行四辺形」という物体を定義してやろう。

data Shape = Point !Position3
           | Plain !Direction3 !Double
           | Sphere !Position3 !Double
           | Parallelogram !Position3 !Direction3 !Direction3 !Direction3
           deriving Eq

Parallelogramがそれ。あとは他の物体と同じく各関数を定義してやる。

getNormal p (Parallelogram _ n _ _) = Just n 

distance r@(pos, dir) (Parallelogram p n d1 d2)
  | res == Nothing = []
  | otherwise      = [t]
  where
    res = methodMoller 2.0 p d1 d2 pos dir
    (u, v, t) = fromJust res

ポリゴンの当たり判定はこのページを 参考にさせていただいた。「Tomas Mollerの交差判定」だそうだ。これを次のように 定義した。

methodMoller :: Double -> Position3 -> Direction3 -> Direction3
             -> Position3 -> Direction3
             -> Maybe (Double, Double, Double)
methodMoller l p0 d1 d2 p d
  | detA == 0.0        = Nothing
  | u < 0.0 || u > 1.0 = Nothing
  | v < 0.0 || v > 1.0 = Nothing
  | u + v > l          = Nothing
  | otherwise          = Just (u, v, t)
  where
    re2 = d <*> d2
    detA = re2 <.> d1
    p'   = p - p0
    te1  = p' <*> d1
    u    = (re2 <.> p') / detA
    v    = (te1 <.> d)  / detA
    t    = (te1 <.> d2) / detA

なお、もともとは三角ポリゴンの判定に使うもので、計算される $u$, $v$ は

 { \displaystyle
0 \leq u \leq 1, 0 \leq v \leq 1
}

 { \displaystyle
u + v \leq 1
}

を満たす必要がある。しかし今回は平行四辺形の当たり判定をしたいので、

$ 0 \leq u \leq 1, 0 \leq v \leq 1 $

$ u + v \leq 2 $

とした。そのため以下のように謎の"2.0"が入っているのだ。

    res = methodMoller 2.0 p d1 d2 pos dir

さあ、準備は整った。画像を生成してみよう。

$ dist/build/pm/pm > scene2.map
$ dist/build/rt/rt < scene2.map > scene2.ppm

生成された画像と同じシーンを古典で生成したもの、および点光源の 画像を並べてみる。

f:id:eijian:20170624204959p:plain

まあ正直なところこの程度の平面光源では点光源と大差ないのだが。

(2) 平行光源=太陽光

このまま勢いに乗って平行光源に手を出そう! 典型的な平行光源は「太陽」だが、その特徴は一方向にしかフォトンが放出 されないことだ。太陽は地球から見れば点光源のようなものだが、あまりに 遠いので地球上ではフォトンの放射方向が(ほぼ)一方向になってしまう。 古典レイトレーシングでは大きさをもたずあらゆる場所で光が届いているか 判定するのだが、フォトンマッピングでは光源の範囲を制限してやらないと 無限にフォトンを放出しなければならなくなってしまう。 ということで今回はちょっとせこい手だが「面光源から一方向にのみフォトンが 放出される」ものを平行光源とすることにした。 イメージは「窓から入り込む太陽光」だ。

f:id:eijian:20170624205011p:plain

とすれば実装は平行光源のそれを参考に比較的簡単にできあがる。

data Light =
  (中略)
  SunLight
    { lcolor :: Color
    , lflux  :: Flux
    , pos    :: Position3
    , nvec   :: Direction3
    , dir1   :: Direction3
    , dir2   :: Direction3
    , ldir   :: Direction3
    }

SunLightがそれ。以下、各関数の定義を示す。

generatePhoton (SunLight c _ p n d1 d2 d0) = do
  wl <- MT.randomIO :: IO Double
  t1 <- MT.randomIO :: IO Double
  t2 <- MT.randomIO :: IO Double
  let r = initRay (p + t1 *> d1 + t2 *> d2) d0
      w = decideWavelength c wl
  return (w, r)

getDirection (SunLight _ _ lp n d1 d2 dt) p
  | cos > 0.0      = []
  | res == Nothing = []
  | otherwise      = [t *> dt']
  where
    d = lp - p
    cos = n <.> d
    dt' = negate dt
    res = methodMoller 2.0 lp d1 d2 p dt'
    (u, v, t) = fromJust res

getRadiance l@(SunLight (Color r g b) f _ _ _ _ _) (d:ds) =
  (Radiance (r * f) (g * f) (b * f)) : getRadiance l ds

generatePhotonは、面のどこからフォトンが放出されるかをランダムに 決めたらあとは指定した一方向へ放出している。 getDirectionはややこしそうだが、実のところ平行四辺形からその地点へ 光が放出されるかどうか、放出方向へ向かって交差判定しているのだ。 なおこの光源も「目に見える」が、窓という想定で空色にしておこう。

これで平行光源の定義もできた。画像を生成してみよう。また、フォトン数を 20万個から30万個に増やした時の画像も示す。まだら模様がいくぶん緩和されて いるのがわかる。影も滑らかで本物っぽい!

f:id:eijian:20170624205023p:plain

古典と比べるとその効果は歴然である。古典では間接光が表現できず ボールの形すら認識できない。一方フォトンマッピングでは間接光により 壁やボールがやんわりと照らされ形になっている!これこそフォトンマッピングの 醍醐味と言えるだろう。ただ、直接壁に当たる光については輪郭がぼやける。 これは点光源での影と同じである一定の範囲のフォトンを収集してしまうため 「直射日光」の当たる少し外側も引っ張られて非常に明るくなってしまうのだ。 フォトンマップ(下図)では輪郭がくっきりしているので、やはり収集方法の 問題だ。

f:id:eijian:20170624205033p:plain

とはいえ、平行光源については概ね満足できる結果といえよう。

3. まとめ

やっと拡散反射を実装できた。手間取ったが、拡散反射により得られる間接光には それだけの効果がある!とくに閉じた部屋へ太陽光(平行光源)が入射した時の表現は 思った以上に綺麗で(個人的には)感動ものだった。

次回は画質の向上をねたにしよう。


ここまでのソースはこちら

DeepLearning(5): スペースリーク解消!

前回で形だけは写経が完成したわけだが、最後に書いたとおりひどいスペースリークが起こって使い物にならないことがわかった。今回はスペースリークの解消に取り組もう。(前回記事ではメモリリークと書いたが、Haskell界隈では(?)スペースリークと言うようだ。どう違うのか、同じなのかよくわからない)

(今回のソース

1. 再帰の改善

いくつかのWebサイトを読んでみたが、Haskellは遅延評価のためにスペースリークを起こしやすいそうだ。素人には細かいことはわからないが、遅延評価のため必要になるまで計算を保留し、その情報をメモリ(ヒープ領域)に置くようだ。だから保留される計算が増えるほどメモリを消費するらしい。再帰といえばスタックオーバーフローだが、それは右再帰での話のようで、、、

さて計算が保留されるような処理がどこにあるのか、ということだがじつはさっぱりわからない。が、どうもfoldlのような左再帰で起こるらしいので、まずはメインループを疑ってみた。もとの実装は次のような感じ。

(Main-train.hs)

main = do

  (中略)

  let
    is = [(epoch_ed st), (epoch_ed st - 1) .. (epoch_st st)]
    getTeachers = getImages (poolT st) (batch st)
    putF = putStatus tm0 (epoch_st st) (epoch_ed st) (opstep st)

  putStrLn "Training the model..."
  putF 0 (layers st) sampleE
  layers' <- trainLoop getTeachers sampleE putF (learnR st) (layers st) is

  (中略)

trainLoop :: (Int -> IO [Trainer]) -> [Trainer]
  -> (Int -> [Layer] -> [Trainer] -> IO ()) -> Double -> [Layer] -> [Int]
  -> IO [Layer]
trainLoop _ _ _ _ ls []             = return ls
trainLoop getT se putF lr ls (i:is) = do
  ls' <- trainLoop getT se putF lr ls is
  teachers <- getT i
  let
    rls = tail $ map reverseLayer $ reverse ls'    -- fist element isn't used
    dls = map (train ls' rls) teachers
    ls'' = update lr (transpose dls) ls'           -- dls = diff of layers
  putF i ls'' se
  return ls''

このtrainLoopがメインループというかメイン再帰であるが、ご覧の通り珍妙なコードになっている。trainLoopは処理の最初でさらにtrainLoopを呼んでいる。これはこのループが学習によってレイヤをどんどん更新していくためで、書いた当時はその更新を表現する方法がこれしか思い浮かばなかったのだ。これが左再帰なのか深く考えていないが、わかりにくく汚いコードには違いないので書き直してみよう。foldlはアキュームレータという"変数"を使って計算結果で更新していくそうだが、これは使えそうだ。

trainLoopはIOがあるのでfoldlを使うことができないがfoldMが使えそうだ。一方で、Haskellfoldlはスペースリークを起こしやすいとの記事もありfoldl'を使えとある。foldl'版のfoldMはないのかと思ったら下記のページに簡単な実装が書いてあった。

StackOverflow: Does Haskell have foldlM'?

これを使って書き直したのが下記のコード。

(Main-train.hs)

main = do

  (中略)

  let
    is = [(epoch_st st) .. (epoch_ed st)]
    getTeachers = getImages (poolT st) (batch st)
    putF = putStatus tm0 (epoch_st st) (epoch_ed st) (opstep st)
    loopFunc = trainLoop' getTeachers sampleE putF (learnR st)

  putStrLn "Training the model..."
  putF 0 (layers st) sampleE
  layers' <- foldM' loopFunc (layers st) is

  (中略)

foldM' :: (Monad m) => ([a] -> b -> m [a]) -> [a] -> [b] -> m [a]
foldM' _ z [] = return z
foldM' f z (x:xs) = do
  z' <- f z x
  z' `seq` foldM' f z' xs

trainLoop' :: (Int -> IO [Trainer]) -> [Trainer]
  -> (Int -> [Layer] -> [Trainer] -> IO ()) -> Double -> [Layer] -> Int
  -> IO [Layer]
trainLoop' getT se putF lr ls i = do
  teachers <- getT i
  let ls' = updateLayers lr teachers ls
  putF i ls' se
  return ls'

updateLayers :: Double -> [Trainer] -> [Layer] -> [Layer]
updateLayers lr ts ls = update lr ls (transpose dls)
  where
    rls = tail $ map reverseLayer $ reverse ls    -- fist element isn't used
    dls = map (train ls rls) ts             -- dls = diff of layers

ついでにtrainLoop'内も小ぎれいにしてコアとなる更新処理をupdateLayersに出しておいた。

trainLoop'再帰がなくなり非常にすっきりした。また、ループのために追番リストを渡しているが、前は逆順だったのが正順に改まった。

で、大事なのはこれでスペースリークが改善するかだが、、、「ダメ」だった。やはりこの問題をきちんと理解せず表面的に再帰処理をどうこうしてもダメなんだろうなと思う。ただ、メインループの気色悪い実装が多少は改善したので、それでよしとしよう。

2. BangPattern

スペースリークの件でググるとBangPatternなるものに言及されている記事がいろいろ見つかる。先述の「計算が保留されてメモリに一時置かれたもの」(これをthunkというらしい)を"とっとと計算しつくす"ように指示するもののようだ(間違ってたらご指摘を)。

ということで、上記のメインループに使ってみた。BangPatternを使うにはソースファイルの先頭にその旨を書くようだ。

(Main-train.hs)
  :
{-# LANGUAGE BangPatterns #-}

  (中略)

foldM' :: (Monad m) => ([a] -> b -> m [a]) -> [a] -> [b] -> m [a]
foldM' _ <font color="red">!z</font> [] = return z
foldM' f <font color="red">!z</font> (x:xs) = do
  z' <- f z x
  z' `seq` foldM' f z' xs

  (中略)

trainLoop' :: (Int -> IO [Trainer]) -> [Trainer]
  -> (Int -> [Layer] -> [Trainer] -> IO ()) -> Double -> [Layer] -> Int
  -> IO [Layer]
trainLoop' getT se putF lr <font color="red">!ls</font> i = do
  :

  (中略)

updateLayers :: Double -> [Trainer] -> [Layer] -> [Layer]
updateLayers lr ts <font color="red">!ls</font> = update lr ls (transpose dls)
  :

foldM'trainLoop'updateLayersの引数についている"!“がその指定。とにかく手当たり次第入れてみた。が、「全滅」した。全く効果がない。。。

thunkをつぶすために deepseq を使うという記事もいくつかあったが、やってみようとすると対象のデータ型をNFData型クラスのインスタンスにする(?)とか、結局素人にはなすすべがなく断念した次第である。

3. hmatrixの採用

効果がありそうなことを試してもちゃんと原因を特定していないので無関係な場所に適用してしまっているのだろう。やっぱりもう少しちゃんと原因を探してみたい。

-1. 問題発生箇所はどこか

スペースリークに気づいたのは畳み込み層を含めた逆伝播処理(Aとしよう)が仕上がった時であり、その前の全結合層だけ逆伝播(Bとする)していた時は気にならなかった。なので両方の実行時のメモリ割り当てを比べてみた。

AもBも処理が進むたびにメモリをどんどん食っていくので、残念ながらスペースリークは前から起こっていたらしい。しかし、Aでは繰り返し1回で約30MB増えていくのに比べ、Bでは約1MBしか増えない。だからBでは500回程度の繰り返しではメモリが足りて問題に気づかなかったようだ。一方Aでは500回も繰り返すと15GB必要になる計算なのでとてもではないが私のPCでは無理だ。

Aの方が極端にメモリを食うので、Aで追加したコード(畳み込み・プーリングの逆伝播)を重点的に調べることにした。特に重そうな畳み込み層の逆伝播処理は次の通り。

(ConvLayer.hs)

deconvolve :: Int -> [FilterC] -> Image -> Delta -> (Delta, Layer)
deconvolve s fs im d = (delta, ConvLayer s (zip dw db))
  where
    delta = convolve s fs $ addBorder s d
    db = map (sum . map sum) d  -- delta B
    sim = slideImage s im
    dw = map (hadamard sim) d

delta =...はほぼ順伝播の処理そのままなので無視、その他を「ダミー(全部0のデータを代わりに与えたり空リストを渡したり)」に変更して実行してみた。

すると、使用メモリの増加が30MBから10MB程度になった!やはりここが悪いらしい。これら処理にはmapsumなど再帰関数が多用されている。データは全て「リスト」で表現しており、処理内容からしてmapsumを多用することは避けられない。どうしようもない?

-2. 思い切ってリストからhmatrixに書き換え

このプログラムでは型(ImageとかFilterCとか)は全てリストで定義してきた。これはまず標準の機能・ライブラリでやれるところまでやってみたかったから。リストでは処理が遅いのはわかっているのでゆくゆくはhmatrixとかRepaとかに移そうと思ってはいたが、速度改善のネタとしてとっておくつもりだったんだがなぁ。。。

だが私のレベルではリストを使いつつ問題を解決することができそうにない。とにかくリストをできるだけ使わないようにするしかないが、代わりに行列ライブラリを使うとして何を使うか。本命はRepaだった。高速化の一環でGPUを使った並列化も考えているため、Accelerateに書き方が似ているからだ(逆か)。

しかしこの記事ではRepaはどうもhmatrixに比べだいぶ遅いようだ。またRepaの解説記事も読んだが書き方が私にはとても難解で奇妙で理解しがたい。一方こちらのhmatrixの解説を読むと素直な書き方ができるようだし、ライブラリのページをみるといろいろ使えそうな関数が揃ってそうだ。ということで、hmatixに決定!

導入と設定は次の通り。

$ cabal install hmatrix

cabalファイルにも追加が必要だ。

(deeplearning.cabal)

  :

executable train
  main-is:             Main-train.hs
  -- other-modules:       
  -- other-extensions:    
  build-depends:       base >=4.8 && <5.0
                     , mersenne-random >= 1.0
                     , containers
                     , time
                     , deepseq
                     , hmatrix

  :

hmatrixは一次元配列相当のVectorと二次元のMatrixからなり三次元以上はなさそう。これらを使って既存の型を定義し直す。

(Image.hs,LayerType.hs)

(before)

type Plain = [[Double]]    -- 2D: X x Y pixels
type Class = [Double]    -- teacher vector
type ActFunc = [Double] -> [Double]
type Kernel = [[Double]]
type FilterF = [Double]

(after)

type Plain = Matrix R    -- 2D: X x Y pixels
type Class = Vector R    -- teacher vector
type ActFunc = Matrix R -> Matrix R
type Kernel = Matrix R
type FilterF = Matrix R

できるだけMatrixVectorに変更した。RDoubleと同じらしい。注意が必要なのがFilterFで、前の実装ではこれは[FilterF]とリストにしていた。今回これをFilterFに改めた。要は[[Double]]Matrix Rに変えたのだ。

この変更により、既存の関数は引数の型合せが大変だったが、それよりも従来リストで書いていたものを全部新しい型を使って全面的に書き直すのに苦労した。特に畳み込み層、プーリング層。これらを見ていこう。

deep learning では行列を使うのが有効だと思うが、前回までのプログラムはあえてリストにしたため、例えば畳み込みで画像の一部を切り取るためにややこしくわかりにくい関数を作ってしまった。

今回型は行列なので、hmatrixの便利な関数subMatrixを使うことでそれが超簡単にできてしまった。逆に、これまで「一部を切り出して畳み込み、次を切り出して畳み込み」と直列処理?をしていたのを、切り取りを最初に一括で行い、それを順に畳み込むよう変更することになった。今回コードの畳み込み処理は次のよう。

(ConvLayer.hs: new)

convolve :: Int -> [FilterC] -> Image -> Image
convolve s fs im = map (convolveImage x iss) fs
  where
    pl = head im
    x  = cols pl - (s - 1)
    y  = rows pl - (s - 1)
    ps = [(i, j) | i <- [0..(x-1)], j <- [0..(y-1)]]
    iss = map subImage im
    subImage :: Plain -> Matrix R
    subImage i = fromColumns $ map (\p -> flatten $ subMatrix p (s, s) i) ps

convolveImage :: Int -> [Matrix R] -> FilterC -> Plain
convolveImage x iss (k, b) = reshape x $ cmap (+ b) $ vsum vs
  where
    vs = zipWith (<#) (toRows k) iss

面倒な行列の切り出しがsubMatrixで済んでしまうことは大きい。前のコードに比べかなり短くなった(28行から14行)。

(参考:前のコードの同じ部分)

(ConvLayer.hs: old)

convolve :: Int -> [FilterC] -> Image -> Image
convolve s fs im = map (convolveImage s im) fs

convolveImage :: Int -> Image -> FilterC -> Plain
convolveImage s im (ks, b) = sumPlains b (zipWith (convolvePlain s) ks im)

sumPlains :: Double -> [Plain] -> Plain
sumPlains b [] = repeat $ repeat b  -- 2 Dimensions (X*Y)
sumPlains b (p:ps) = zipWith f p (sumPlains b ps)
  where
    f :: [Double] -> [Double] -> [Double]
    f [] _          = []
    f (a:[]) (b:_)  = [a + b]
    f (a:as) (b:bs) = (a + b) : f as bs

convolvePlain :: Int -> [Double] -> Plain -> Plain
convolvePlain s k ps
  | len < s   = []
  | otherwise = convolveLine s k p : convolvePlain s k ps'
  where
    len = length ps
    p   = take s ps
    ps' = tail ps

convolveLine :: Int -> [Double] -> [[Double]] -> [Double]
convolveLine s k ps
  | len < s   = []
  | otherwise = sum (zipWith (*) vs k) : convolveLine s k ps'
  where
    len = minimum $ map length ps
    vs  = concatMap (take s) ps
    ps' = map tail ps

プーリング層も同様で、行列を小領域に区分けして処理していくため、subMatrixが活躍する。

こちらも始めに小領域を全部作って順に処理する形に変えたので、行数は31行から19行にへった。

(PoolLayer.hs: new)

poolMax :: Int -> Image -> [Image]
poolMax s im = [os, is] 
  where
    pl = head im
    x  = cols pl `div` s
    y  = rows pl `div` s
    ps = [(i*s, j*s) | i <- [0..(x-1)], j <- [0..(y-1)]]
    (os, is) = unzip $ map (toPlain x y . maxPix s ps) im

toPlain :: Int -> Int -> [Pix] -> (Plain, Plain)
toPlain x y pls = ((x><y) op, (x><y) ip)
  where
    (op, ip) = unzip pls

maxPix :: Int -> [(Int, Int)] -> Plain -> [Pix]
maxPix s ps is = map (max' . toPix) ps
  where
    toPix :: (Int, Int) -> [Pix]
    toPix p = zip (concat $ toLists $ subMatrix p (s, s) is) [0.0..]

(参考:前の実装はこう)

(PoolLayer.hs: old)

poolMax :: Int -> Image -> [Image]
poolMax s im = [fst pixs, snd pixs] 
  where
    pixs = unzip $ map unzipPix (poolMax' s im)

unzipPix :: [[Pix]] -> ([[Double]], [[Double]])
unzipPix pixs = unzip $ map unzip pixs

poolMax' :: Int -> Image -> [[[Pix]]]
poolMax' s = map (poolMaxPlain s)

poolMaxPlain :: Int -> Plain -> [[Pix]]
poolMaxPlain s p = map (poolMaxLine s) ls
  where
    ls = splitPlain s p

splitPlain :: Int -> Plain -> [[[Double]]]
splitPlain s [] = []
splitPlain s p  = p' : splitPlain s ps
  where
    (p', ps)  = splitAt s p

poolMaxLine :: Int -> [[Double]] -> [Pix]
poolMaxLine _ [] = []
poolMaxLine s ls
  | len == 0  = []
  | otherwise = pixs : poolMaxLine s ls'
  where
    len  = length $ head ls
    pixs = max' $ zip (concatMap (take s) ls) [0.0 ..]
    ls'  = map (drop s) ls

hmatrixの使用前後を比べてみるとわかるが、かなり作り方というか処理の流れが変わってしまった。関数は、ほとんど一から作りなおしたようなものだ。

ここで、テストコードに大いに助けられた。私はズボラなのでテストコードを書くのは嫌いで苦手だが、deep learningでは思った通りの変換ができているか確認しておかないと各層の実装間違いが後に響くため、いつもより多めにテストコードを入れていた。

前述の通り一部の関数はほぼ作り直しだったが、テストで得られる結果は前の実装と同じはずだ。データ型を変更したので入力・出力の形は異なるが結果の値は同じだ(見難いので一部改行を入れてある)。

(old code)

>>> let filter = [
  ([[1.0,2.0,3.0,4.0],[5.0,6.0,7.0,8.0]],0.25),
  ([[3.0,4.0,5.0,6.0],[7.0,8.0,1.0,2.0]],0.5),
  ([[5.0,6.0,7.0,8.0],[1.0,2.0,3.0,4.0]],0.75)] :: [FilterC]
>>> let img1 = [
  [[1.0,2.0,3.0],[4.0,5.0,6.0],[7.0,8.0,9.0]],
  [[4.0,5.0,6.0],[7.0,8.0,9.0],[1.0,2.0,3.0]]] :: Image
>>> convolve 2 filter img1
[[[200.25,236.25],
  [173.25,209.25]],
 [[152.5,188.5],
  [233.5,269.5]],
 [[152.75,188.75],
  [197.75,233.75]]]

(new code)

>>> let filter = [
  (fromLists [[1.0,2.0,3.0,4.0],[5.0,6.0,7.0,8.0]],0.25),
  (fromLists [[3.0,4.0,5.0,6.0],[7.0,8.0,1.0,2.0]],0.5),
  (fromLists [[5.0,6.0,7.0,8.0],[1.0,2.0,3.0,4.0]],0.75)] :: [FilterC]
>>> let img1 = [
  fromLists [[1.0,2.0,3.0],[4.0,5.0,6.0],[7.0,8.0,9.0]],
  fromLists [[4.0,5.0,6.0],[7.0,8.0,9.0],[1.0,2.0,3.0]]] :: Image
>>> convolve 2 filter img1
[(2><2)
 [ 200.25, 236.25
 , 173.25, 209.25 ],(2><2)
 [ 152.5, 188.5
 , 233.5, 269.5 ],(2><2)
 [ 152.75, 188.75
 , 197.75, 233.75 ]]

このように、テストで前と同じ結果であることを保証できないと内部実装をガラッと変えるのは怖い。今後テストコードをこまめに書くモチベーションアップにはなったと思う。

-3. 評価

(1) 学習結果(認識精度)

改善したプログラム(新CNNと呼ぼう)を例によって10回実行し、それぞれ認識精度と実行時間を計測した。前回は10回の平均を取ったが、これだと立ち上がりのタイミングが大きくぶれるため平均したら実際よりなだらかな微妙な曲線になってしまう。そこで今回は平均ではなく中間値でグラフを書いてみることにした。前回の、全結合層のみの値と畳み込み層までの値も同じく中間値でプロットしてみた。

f:id:eijian:20170205114950p:plain

横軸が処理時間なので、新CNNの方が全結合層のみの場合より立ち上がるタイミングが遅いものの、学習効率が高く早期に95%を超えるのがわかる。今回の新CNNでは学習の進み方が圧倒的に速いためとても性能が高いと思う。

前回同様、処理にかかる時間を比較表にした。今回は500 epoch まで計測できたのでそれも追加した。

新CNN CNN 全結合層のみ学習
処理時間(s) 177 517 157
処理性能(epoch/s) 1.130 0.387 1.273
認識精度(%) 200epochs 97.8 97.5 74.2
認識精度(%) 500epochs 99.6 - 95.1

単位時間あたりの処理量が、かなり全結合層のみ学習の時のそれに近づいた。一方で認識精度はさすが、かなり良い。

ということで、hmatrixを使った今回のバージョンでやっと本当の意味でのyusugomori実装の「写経」が完成したと言えるだろう。

(3) スペースリーク

そうそう、今回はスペースリークの解消が目的であった。新CNN実行時のメモリ使用量をtopコマンドで眺めていたが、繰り返し回数によらずだいたい60MBから70MBの間で推移し、多少増減するものの70MBを超えて増加することなく安定していた。スペースリークは解消したようだ!

-4. hmatrixを使った感想

このdeep learningのプログラムを実装するにあたり初めてhmatrixを知ったのだが、以下の点でこのライブラリをとてもよい。

  • わかり易い: 行列を作ったり、いろいろな操作をするのに込み入った記述はいらない。例えば 2x3行列を作りたければ (2><2) [要素のリスト]と書けばよい。

  • Numクラスのインスタンス(?): 足し算(+)、引き算(-)の関数はNumクラスのインスタンスでないと実装できないが hmatrixのMatrixVectorはそのインスタンスなので普通に足したり引いたりできる。以前別のプログラムでNumericPreludeを使って四則演算を自作したことがある。hmatrixはどうやっているのだろうか。今度ソースを眺めてみよう。

  • 便利関数多し: 高校や大学で習うような行列・ベクトルの演算はいわずもがな、CNNを作るのに「こんな関数が欲しいなぁ」と思うとちゃんと用意されている。subMatrixはもちろん、要素を全部足すsumElementsや2つのベクトルの要素をそれぞれ掛けた値を行列にするouterなど、deep learningを実装するために作られたのではと思えるほど。なお、hmatrixを使った実装が終わり評価もほぼ終わったところでcorr2という関数を見つけた。これはどうやら「畳み込みをする関数」のようだ・・・。せっかく実装したのだが・・・。しかしcorr2を使うにはデータ型を変更しなくてはいけないみたいなので今回はパス。次の改善で考えよう。

4. おまけ

これまで何回か処理性能を評価してきたが、試行毎に学習の立ち上がりや認識精度の改善度合いに大きなばらつきがあることがわかった。試行毎の違いは全結合層および畳み込み層のフィルタの初期状態による。初期状態は乱数で設定するが、その値の範囲は次のように設定している。($x$をフィルタの設定値、$in$と$out$はその層の入力数・出力数とする)

(全結合層)

 {\displaystyle
\left|x\right| \leq \frac{1.0}{in}
}

(畳み込み層)

 {\displaystyle
\left|x\right| \leq \sqrt{\frac{6}{in+out}}
}

一方でこの記事によると、deep learningの論文では学習効率を上げるためのフィルタの初期値として次の式が示されているという(私は原論文を読んでいないが)。

(活性化関数が$\tanh$のとき)

 {\displaystyle
\left|x\right| \leq \sqrt{\frac{6}{in + out}}
}

(活性化関数が$\mathrm{sigmoid}$のとき)

 {\displaystyle
\left|x\right| \leq 4 \cdot \sqrt{\frac{6}{in + out}}
}

特に全結合層では大きく違いがある。これを論文で示された値にしたらどうなるか全結合層の実装で下記のように試してみた。

initFilterF :: Int -> Int -> IO FilterF
initFilterF k c = do
  rs <- forM [1..k] $ \i -> do
    w <- forM [1..c] $ \j -> do
      r <- MT.randomIO :: IO Double
      return ((r * 2.0 - 1.0) * a)
    return (0.0:w)
  return $ fromLists rs
  where
    a = 1.0 / fromIntegral c

    ↓

    a = 4.0 * sqrt (6.0 / fromIntegral (c + k))

実行結果を表にしたものがこちら。

f:id:eijian:20170205114957p:plain

効果は歴然、学習効率が格段に向上していることがわかる。フィルタの初期値の違いがこれほど大きな差を生むとは驚きだ。逆に言えば、初期値を適切に与えないとなかなか学習が進まないということでもあり、おそらく個々の問題毎に最適な値が違うと思うので、これはもう試行錯誤しかないのだろう。ただ上記の論文ではその中でも比較的良い結果を得る指針かと。

5. まとめ

今回の改善で、やっとまともに動作するyusugomori実装の「写経」が完成し、初回に書いたステップ1が完了した。さすがに深層学習を一から勉強しないといけなかったのでしんどかった。。。

それでは次回、MNISTデータを使った学習(ステップ2)に取り組もう。

DeepLearning(4): CNNの逆伝播完成?

また時間が経ってしまった。畳み込み層の逆伝播のアルゴリズムを 理解するのに手間取ってしまった。。。気を取り直して、逆伝播の後半戦といこう。

1. 各層のおさらい

今回作っているプログラムは次のように11層のレイヤで画像を学習・評価している。

No. 実装
0 入力 -
1 畳み込み層1 ×
2 活性化層1
3 プーリング層1 ×
4 畳み込み層2 ×
5 活性化層2
6 プーリング層2 ×
7 平坦化層 ×
8 全結合層1
9 活性化層3
10 全結合層2
11 活性化層4(出力)

前回までで全結合層の逆伝播処理は説明したので、それらの処理は実装済み、 残るは×のついている層だ(平準化層、プーリング層、畳み込み層)。

その前に、前回のプログラムでは逆伝播において伝播させる誤差($\delta$)を 次のようにDoubleの一次元配列(というかリスト)とした。

type Delta = [Double]

しかしプーリング層や畳み込み層では、誤差は単純な一次元配列にならず、画像と 同じくXY平面がNチャネルという三次元のデータになる。ということで、次のように 改めた。すでに実装済みとした部分もこれに合わせて変更してある。

type Delta = Image  -- (= [[[Double]]])

2. 各層の逆伝播処理

2-1. 平坦化(Flatten)層

前回は全結合層までの実装だったので、上記の表で言えば No.11 → No.8 までしか できていなかった。その次の平坦化層(No.7)はCNNの醍醐味である畳み込み層へ誤差を伝えるための 重要な転換点だ。が、実装は極めて単純だった。

unflatten :: Int -> Int -> Image -> Image
unflatten x y [[ds]] = split y $ split x ds
  where
    split :: Int -> [a] -> [[a]]
    split _ [] = []
    split n ds = d : split n ds'
      where
        (d, ds') = splitAt n ds

この層の逆伝播処理(unflatten)は、下層から一次元配列で伝わってきた誤差を 画像と同じ三次元配列に組み替えて上層に伝えること。 画像一枚の画素数(x,y)は引数で与える。 第三引数が微妙な形[[ds]]になっているのはDelta型を変更したから。 dsの型は[Double]だ。

ちなみに、入力されるデータ数がx,y画素できちんと割り切れることを想定して いて何のエラー処理もしていない。。。

2-2. プーリング層

プーリング層にはフィルタ情報がないので、下層から伝播してきた誤差($\delta$)を 更新して上層へ伝えれば良い。ただし、プーリング層の更新処理は特殊で、 他の層では誤差に何らかの計算処理をすればよいが、プーリング層では次のような ことをしないといけない。

  • 順伝播の際に小領域(2x2とか3x3とか)から最大値を選択し、その画素の位置を 覚えておく。
  • 逆伝播時は、誤差を順伝播で選択した画素の位置に分配し、小領域のそれ以外の画素には 0を設定する。

これを図にしたらこんな感じ。

f:id:eijian:20170108203040p:plain

前回のプログラムでは、プーリング層への入力画像は記録して いるものの、どの画素を選択したかという情報は記録していなかった。 yusugomori実装では、 選択画素を記録しておく代わりにプーリング層の「出力値=選択された画素」の値を使って、 小領域の値と順次比較することでどの画素を選択したかを特定するようになっている。 つまり誤差を求めるために「プーリング層の入力値と出力値の両方を使っている」 のだ。

しかし私の実装では、逆伝播を再帰で綺麗に(?)処理しようとしたために 誤差の計算は「各層の入力値のみ使う」形にしている。このシンプルな構造は 維持したい・・・。

結局「プーリング層の入力値を改竄する」ことにした(笑)。そのために 順伝播での処理をこう変えた。

  • 画像を小領域に分割する(poolMaxLine)
  • 小領域の各画素に番号を付ける(zipで画素値と追番を組(Pix)にする)
  • 小領域中の画素の最大値を選択する(max')
  • 選択した情報を集めると、(最大値,小領域中の追番)の組の集合になる
  • 組の第一要素は下層へ伝える出力値、第二要素は「最大値画素の小領域内の 追番(の集合)」になる(poolMax)
  • プーリング層は引数でもらった入力値をリストから削り、代わりに 「追番の集合」と出力値を返す

プログラムはこうなった。

type Pix = (Double, Double)

   :

poolMax :: Int -> Image -> [Image]
poolMax s im = [fst pixs, snd pixs] 
  where
    pixs = unzip $ map unzipPix (poolMax' s im)

   :

poolMaxLine :: Int -> [[Double]] -> [Pix]
poolMaxLine _ [] = []
poolMaxLine s ls
  | len == 0  = []
  | otherwise = pixs : poolMaxLine s ls'
  where
    len  = length $ head ls
    pixs = max' $ zip (concatMap (take s) ls) [0.0 ..]
    ls'  = map (drop s) ls

max' :: [Pix] -> Pix
max' [] = error "empty list!"
max' [x] = x
max' (x:xs) = maximum' x (max' xs)

maximum' :: Pix -> Pix -> Pix
maximum' a@(v1, i1) b@(v2, i2) = if v1 < v2 then b else a

Pixの組の定義で追番がDoubleなのは、追番の集合を 「画像」データとして保持するため型を合わせたから。poolMaxで 組のリストから「出力値」と「追番集合」を分離して返している。 poolMaxLine中の

    pixs = max' $ zip (concatMap (take s) ls) [0.0 ..]

が各画素に追番を付けて最大値を選択しているところだ。

forwardLayer :: Layer -> Image -> [Image]
forwardLayer (NopLayer)         i = [i]
forwardLayer (ActLayer f)       i = [activate f    i, i]
forwardLayer (MaxPoolLayer s)   i =  poolMax  s    i
forwardLayer (ConvLayer s fs)   i = [convolve s fs i, i]
   :

forwardLayerにおいてプーリング層だけ他と異なっているのがわかるだろう。 他は全て引数のImage(=i)をその層で計算された出力値の後ろにくっつけて 返している。プーリング層では"i"をくっつけずPoolMaxの戻り値を使う。 これが「入力値を削る」という意味だ。

ここまでが順伝播での「準備段階」。逆伝播の処理は次の通り。

  • 誤差と追番を組にする(depoolMax内のconcatWith)
  • その組から小領域を生成する(追番のところは誤差値、それ以外は0)(expand)

順伝播で頑張った分、逆伝播は単純になった。プログラムではこう。

depoolMax :: Int -> Image -> Delta -> (Delta, Layer)
depoolMax s im d = (zipWith (concatWith (expand s 0.0)) im d, MaxPoolLayer s)
  where
    concatWith :: ([Int] -> [Double] -> [[Double]])
               ->  [[Double]] -> [[Double]] -> [[Double]]
    concatWith f i d = concat (zipWith f (map (map truncate) i) d)

expand :: Int -> Double -> [Int] -> [Double] -> [[Double]]
expand s r ps ds = map concat (transpose $ map split' $ zipWith ex ps ds)
  where
    split' = split s
    ex :: Int -> Double -> [Double]
    ex p d = take (s*s) (replicate p r ++ [d] ++ repeat r)

split :: Int -> [a] -> [[a]]
split s [] = []
split s xs
  | length xs < s = [xs]
  | otherwise     = l : split s ls
    where
      (l, ls) = splitAt s xs

concatWithは少々ややこしいが追番集合を実数から整数に変換して誤差と 組にしたうえでexpandに渡している。

これで誤差を上層へ伝播できるわけだ。

2-3. 畳み込み層

(1) 誤差の計算

畳み込み層の逆伝播処理についてわかりやすく説明している書籍やWebサイトを 見つけられなかったので理解するのにかなり時間がかかってしまった。 唯一yusugomoriさんのblogが あったが、数式が苦手な私には理解は厳しく、結局yusugomori実装のソースを 追いかけることになった。。。

そしてやっと、逆伝播における畳み込み層の誤差計算の仕組みがわかった(つもり)。 順伝播と比較するとわかりやすいので図にしてみた。 (専門の方には初歩すぎるだろうが、本文章は自身の忘備録を兼ねて いるため軽く流していただきたい)

f:id:eijian:20170108203047p:plain

この例はフィルタ( $w$ )が2x2の場合である。まず順伝播では、入力画像の4つの画素から 出力の画素値を決めている。どの入力画素がどのフィルタ値と掛けられたか矢印を 色分けした。たとえば $x_{11}$ は $w$ の1番目(赤、 $w_1$ としよう)を掛けて $y_{11}$ の一部となっているので矢印は赤。 $y_{11}$ を式で書くと、

{ \displaystyle
y_{11} = x_{11} \cdot w_1 + x_{12} \cdot w_2 + x_{21} \cdot w_3 + x_{22} \cdot w_4
}

である。一方右の逆伝播の図では、今度は $y_{11}$ が下層から伝播された誤差である。 その誤差から上層へ伝える誤差の一つ $x_{11}$ への矢印はやはり赤であるから $w_1$ を掛けるのだ。つまり順伝播と逆をしているだけだ。

しかしプログラム上は $w$ の逆変換(?)である $w'$ (全結合層における $W$ の転置に相当)を 求めて使っている。これは「順伝播と同じ畳み込み処理をそのまま使って誤差を求める ことができる」からだ(と理解した)。

畳み込み層の順伝播処理では、出力画像の解像度は 入力画像のそれより小さくなる。例えばフィルタが2x2で入力が5x5ドットとすると 出力は(5-(2-1))x(5-(2-1))=4x4だ。逆伝播では下層から伝播した誤差データの 解像度(例えば4x4)が、求めたい誤差データの解像度(たとえば5x5)より小さいので そのままでは順伝播と同じ処理が使えない。そこで、誤差データに「外枠」を付けて 求めたい解像度より大きなサイズにする。右図では上と左に1ドットずつ"0"が並んでいる。 これがその「外枠」で、上下左右に同じ幅で追加する必要がある。幅はフィルタサイズ で決まる。2x2ならそれぞれ-1して1ドットずつ、3x3ならやはり-1して2ドットずつだ。

ここで、右図の左上の青点線で囲った中を見てみる。各画素から $x_{11}$ が 得られているが、 $y_{11}$ 以外は0なので $x_{11} = y_{11} \cdot w'_4 = y_{11} \cdot w_1$ だ。 左の順伝播での青点線の中のちょうど逆だ。一方、右図の赤点線の中は $x_{23}$ を 周りの $y$ からフィルタ $w'$ を掛けて求めている。これも順伝播での $x_{23}$ の 逆なだけ、ただ赤点線の中だけみるとフィルタが順伝播の逆順になっている。 これが $w$ の逆変換 $w'$ を使う理由。この $w'$ さえ用意してしまえば、あとは 順伝播の処理をそのまま呼び出せばよい。プログラムを見てみよう。

deconvolve :: Int -> [FilterC] -> Image -> Delta -> (Delta, Layer)
deconvolve s fs im d = (delta, ConvLayer s (zip dw db))
  where
    delta = convolve s fs $ addBorder s d
      :

この最後のdeltaが誤差を計算しているところ。fsには予め計算しておいた $w'$ が 入っている。下層からの誤差daddBorderで「外枠」を付けたら、 順伝播のconvolveを呼び出すだけ。ネタがわかってしまえば プログラム自身はとても簡単であった。

(2) 勾配の計算

次に勾配( $\Delta W$ と $\Delta b$ )の計算について考える。 まず簡単な方から。 $\Delta b$ は、下層からの誤差を画像の各チャネルで 集計するだけ。各チャネルはX,Yの二次元データなのでそれを全て足せば良い。 式で書くとしたら下記のような感じか($c$はチャネル番号とする)。

{ \displaystyle
\Delta b^{(c)} = \sum_{x=1}^N \sum_{y=1}^M \delta_{xy}^{(c)}
}

$\Delta W$ の方はちょっとややこしい。下図を見て欲しい。

f:id:eijian:20170108203054p:plain

これは誤差、入力画像とも、とあるチャネルを取り出して $\Delta w$ (一つのフィルタということで小文字)を求める場合だ。それぞれ二次元データに なっている。要は $\delta$ と $a$ の該当する画素を掛け合わせて全て足すのだ。 ベクトルの内積のような感じ(行列も内積と言うのかな?)。 たとえば $\Delta w_1$ (赤)を求めるには、

{ \displaystyle
\Delta w_1 = \sum_{x=1}^4 \sum_{y=1}^4 \delta_{xy} \cdot a_{xy}
}

を計算する。それ以外の要素もほぼ同じだが、 $a$ をフィルタの各要素の 場所ぶんだけずらした値を使う。各色の点線で囲った部分を使ってフィルタの 同じ色の要素を計算するわけだ。

これを誤差のチャネル数( $k$ とする)と入力画像のチャネル数( $c$ とする)を 掛けた $k \times c$ 枚計算すれば $\Delta W$ が出来上がる。

プログラムではこうなっている。

deconvolve :: Int -> [FilterC] -> Image -> Delta -> (Delta, Layer)
deconvolve s fs im d = (delta, ConvLayer s (zip dw db))
  where
    delta = convolve s fs $ addBorder s d
    db = map (sum . map sum) d  -- delta B
    sim = slideImage s im
    dw = map (hadamard sim) d   -- delta W

slideImage :: Int -> Image -> [[Plain]]
slideImage s im = map (concat . slidePlains s . slidePlain s) im

(slide部分は省略)

hadamard :: [[Plain]] -> Plain -> [[Double]]
hadamard sim ds = map (map (mdot ds)) sim

まずslideImageで上記説明した各色の枠内の $a$ を取り出している。 ただし、右や下の使わない部分は、誤差dと突き合わせる際に無視されるため わざわざ切り取っていない(例えば赤点線内は $a_{11}$ から $a_{14}$ までだがslideImageで得られるデータには $a_{15}$ も含まれている)。 そして、hadamardで誤差dとスライドした入力値を掛け合わせて $\Delta W$ を求めている。

(3) フィルタ更新

(2)までで、ある教師データについて逆伝播の処理内容を説明した。 あとはバッチ1回で複数の教師データから得られた各 $\Delta W_i, \Delta b_i$ から フィルタの更新をすればよい。これは全結合の時と同じ式で計算できる。

{ \displaystyle
\boldsymbol{W} \leftarrow \boldsymbol{W} - \eta \frac{1}{N} \sum_{i=1}^{N} \Delta \boldsymbol{W}_i
}

{ \displaystyle
\boldsymbol{b} \leftarrow \boldsymbol{b} - \eta \frac{1}{N} \sum_{i=1}^{N} \Delta \boldsymbol{b}_i
}

この部分はプログラムではこうしている。

updateConvFilter :: Int -> [FilterC] -> Double -> [Layer] -> Layer
updateConvFilter s fs lr dl = ConvLayer s $ zip ks' bs'
  where
    sc = lr / fromIntegral (length dl)
    (ks , bs ) = unzip fs
    (kss, bss) = unzip $ map unzip $ strip dl
    dbs = vscale sc $ vsum bss
    dks = map (mscale sc . msum) $ transpose kss
    bs' = vsub bs dbs
    ks' = zipWith msub ks dks

strip :: [Layer] -> [[FilterC]]
strip [] = []
strip (ConvLayer s fs:ds) = fs:strip ds
strip (_:ds) = strip ds  -- if not ConvLayer

引数fsが更新前のフィルタ($W$や$b$)、dlに各教師データから 学習で得られた $\Delta W$ と $\Delta b$ が入っている。これを上式に 当てはめて計算するのだ。

さあ、プログラムは完成した!

3. 評価

それでは早速出来上がったプログラムを実行して学習をさせてみよう。 (ソースはここ)

$ cabal clean
$ cabal build
  :
$ ./dist/build/train/train
Initializing...
Training the model...
iter =     0/200 accuracy = 0.3333333333 time = 0.13968s
iter =     5/200 accuracy = 0.3333526481 time = 11.021087s
iter =    10/200 accuracy = 0.3333741834 time = 21.513648s
iter =    15/200 accuracy = 0.3334006436 time = 33.707092s
iter =    20/200 accuracy = 0.3334352183 time = 45.034608s
iter =    25/200 accuracy = 0.3334818706 time = 55.7984s
  :

3-1. 前回分(全結合層のみ)との比較

(1) 学習結果(認識精度)

上記のプログラムを10回実行し、それぞれ認識精度と実行時間を計測した。 今回epoch数は200(前回は500)としたが理由は後述する。結果は次のグラフの とおり。

f:id:eijian:20170108203100p:plain

点線が各試行時の実測値、太いオレンジの実線が10回の平均だ。各試行では 認識精度が向上する時期(epoch)に大変ばらつきがある。これは全結合層の 学習でも同じだったが、最初にランダムに設定するフィルタ初期値が良くない からだと思われる。このように向上する時期がバラバラなので、 平均を取るとなんだか緩やかに学習が進むように見えてしまう。。。

次にyusugomori実装と比較してみる。上記の平均がオレンジ、yusugomori実装は 赤だ。参考までに前回の全結合のみの場合の平均(青)も比較のため入れた。

f:id:eijian:20170108203108p:plain

オレンジは上述の通り平均値のためなんだかふらついて変だが、最終的には ほぼ同程度の認識精度が得られている。前のグラフを見れば、認識精度の 向上具合はyusugomori実装と比べても遜色ないと思われる。 (実は逆伝播処理の「値の正しさ」はきちんと検証していないのだが、 この結果をみるとそんなに外していないのかな?ということでよしとする)

(2) 処理時間

ただ、いいことばかりではない。 畳み込み層まで逆伝播させたため、学習にかかる時間は逆に増えてしまった。 200 epoch実行時の結果を次表に示す。

CNN 全結合層のみ学習
処理時間(s) 517 157
処理性能(epoch/s) 0.387 1.273
認識精度(%) 200epochs 97.5 74.2
認識精度(%) 500epochs - 95.1

このように、同じ時間で処理できるepoch数が1/3以下に落ちてしまった。 さすがに複雑な畳み込み層/プーリング層の逆伝播処理は負荷が高い・・・。 同じepoch数ではさすがにCNNの方が認識精度が格段に高い。しかし全結合層のみの 学習でも500 epochs続ければ95%を超えるわけで、結局時間と得られる精度は どうなるかと思ってプロットしてみたのが次のグラフ。

f:id:eijian:20170108203114p:plain

このように、処理時間と認識精度を秤にかけると、畳み込み層まで逆伝播させたときは、 得られる精度に対しあまりアドバンテージがないように思える。もちろん、 認識精度を99%以上まで非常に高くするためには畳み込み層の学習が 必要なのだろうと思うが。

一方で、当然ながら私の実装のまずさもあるだろう。やたら再帰を使うとか 自分でも頭が混乱するような処理になっているので、もっと簡潔な実装に 改められたら改善するかもしれない。今後の課題の一つだ。

3-2. メモリリーク!?

さて、表題に「逆伝播完成?」とクエスチョンを付けた理由がこれだ。 epoch数を200に下げた理由も。

やっとできたと思ってepoch数500で実行してみるとなかなか終わらない。 開始時は1 epochの処理に約2秒かかったので単純計算では1000秒になるが 30分経っても350程度。おかしいと思ってtopコマンドでプロセスを見て みたらメモリを10GBも食っていたのですぐ停止した!

気を取り直してもう一度実行したときのtopコマンドの出力がこれ。train プロセスが大量のメモリを消費しているのがわかる。これは繰り返し回数が だいたい180ぐらいのときの状態。

Processes: 162 total, 3 running, 159 sleeping, 639 threads             00:22:32
Load Avg: 1.99, 1.98, 1.80  CPU usage: 8.48% user, 32.11% sys, 59.40% idle
SharedLibs: 92M resident, 33M data, 6516K linkedit.
MemRegions: 23895 total, 477M resident, 29M private, 103M shared.
PhysMem: 4080M used (1131M wired), 14M unused.
VM: 471G vsize, 623M framework vsize, 71855224(48041) swapins, 73935000(0) swapo
Networks: packets: 631497/460M in, 456172/62M out.
Disks: 2706786/300G read, 1541857/302G written.

PID    COMMAND      %CPU TIME     #TH   #WQ  #PORT MEM    PURG   CMPRS  PGRP
25694  train        81.4 08:49.94 1/1   0    10    1064M- 0B     4044M+ 25694
0      kernel_task  76.5 02:28:47 114/7 0    2     443M-  0B     0B     0
  :

Macのメモリ管理についてはよく知らないが、trainプロセスの使用メモリが 1064M、CMPRSが4044Mで、合わせて5100M=5GBほども消費している。 これはいただけない。

本実装ではやたら(無駄に?)再帰を使ってしまっているが、その効率が 悪いのかもと思いつつ、消費メモリが増える一方なのはヒープ領域を やたら使っているのだろう。ググると、foldl'を使えとかサンクがどうとか いろいろあるみたいだが、簡単には解決しそうになかったので今回はここまで。 対策は次回に頑張ろうと思う。

4. まとめ

今回でやっと一通りCNNが 実装 できた。学習能力も手本としたyusugomori実装と同等のものができた。 しかしながら、上述のとおりメモリを大量に消費するという問題を抱えたままで、 このままでは全く実用に耐えない。

次回はこの問題を解決させたいと思う(できないかも・・・)。

DeepLearning(3): そして逆伝播(でも全結合層まで)

第1、2回で順伝播の処理は実現した。今回やっと「学習」処理に入る。 (それにしても時間が取れず、またしてもだいぶあいだが空いてしまった。記憶が飛んでいる)

前回作成したプログラムでは、多段のニューラルネットワークを構成したが、 その終盤には2つの全結合層(隠れ層と出力層)が含まれている。本稿では、 逆伝播による学習で全結合層のフィルタ(ウェイト)を更新する部分について 解説し、学習によりフィルタの分類能力が向上することを確認したい。

今回のソースはこちら

1. 全結合層の逆伝播処理

1-1. 処理の流れ

逆伝播の理論や数式の詳細は各種書籍やWebサイトを見ていただくとして、 コード実装の説明に必要な事項をざっと説明しよう。 畳み込み層の逆伝播はちょっと複雑っぽいので次回にし、今回は全結合層だけに しておく。

損失関数

「学習」とは、例えば画像認識では入力画像から「犬」や「猫」を正しく判断 できるように各層のフィルタを修正してくことだ。 (だから本プログラムでは全結合層と畳み込み層がその対象なのだが、 今回は全結合層だけ。)

修正するにはまず「どれぐらい正解とずれているか」を確認しないといけない。 これを損失関数 $E$ で測る。今回は出力層の活性化関数にsoftmaxを使っており、 出力値と正解値の誤差を簡単にするため $E$ には交差エントロピー関数を 採用した。

{ \displaystyle
E = -\sum_k t_k \ln y_k \qquad (1)
}

ここで、$t$ は正解値(ベクトル)、 $y$ は出力値(ベクトル)だ。例えば 本プログラムでは画像を3種類に分類するのでベクトルの要素数は3、 仮に教師データを分類1とするなら、

$t = [1.0, 0.0, 0.0], y = [0.7, 0.1, 0.2]$

といった感じだ。

フィルタの更新

さて損失関数で表される「ずれ」を小さくするため、フィルタ $\boldsymbol{W}$ を更新する。大雑把には次の通り。

{ \displaystyle
 \boldsymbol{W} \leftarrow \boldsymbol{W} - \eta \Delta \boldsymbol{W} \qquad (2)
}

ここで、 \eta は学習率(learning rate)と呼ばれるパラメータだ。 さて問題は $\Delta \boldsymbol{W}$ をどう求めるかだが、損失関数 $E$ の勾配で表す ことができるようだ。

全結合層のフィルタ $\boldsymbol{W}$ は入力数 $\times$ 出力数の行列と見ることができる。 その個々の要素を $w_{ij}$ とすると、

 { \displaystyle
w_{ij} \leftarrow w_{ij} - \eta \frac{\partial E}{\partial w_{ij}} \qquad (3)
}

と書ける(らしい)。だから、最後の偏微分を解けばめでたしめでたしとなる(が、 そこが大変)。直接この偏微分を各層で解くと本当に大変だそうだが、 実は一つ前の層(出力に近い層)の計算結果をうまく使うとこれが簡単にできる、 というのが「誤差逆伝播法」のミソらしい。つまり出力層側から 誤差を計算して次の層に渡す。それを元に次の誤差を計算して さらに次の層に渡す。これを繰り返す。

出力層(Lとする)では、上記の偏微分は以下のようになる。

{ \displaystyle
\frac{\partial E}{\partial w_{ij}^{(L)}} = \frac{\partial E}{\partial y_i^{(L)}}
 \frac{\partial y_i^{(L)}}{\partial a_i^{(L)}}
 \frac{\partial a_i^{(L)}}{\partial w_{ij}^{(L)}}
= \delta_i y_j^{(L-1)} \qquad (4)
}

{ \displaystyle
\left( \delta_i = \frac{\partial E}{\partial y_i^{(L)}}
 \frac{\partial y_i^{(L)}}{\partial a_i^{(L)}} \right)
}

この $\delta_i$ が次の層へ伝える「誤差」、これに一つ手前の層の結果 $y_j$ を掛ける。

順伝播/逆伝播の流れと数式

これらを踏まえ、各層の順伝播・逆伝播を表にしてみた (逆伝播は、今回は全結合層まで)。 右上のカッコ付き数字は層の番号。 太字はベクトル/行列で、 $\boldsymbol{W}_{[3,2]}$ は 3x2行列、のように右下に要素数を記載した。 なお、 $\boldsymbol{c}$ は教師データである。 (層番号と要素数をつけたら数式が非常に醜くなってしまった。。。)

No. 順伝播 逆伝播 勾配($\Delta \boldsymbol{W}$)
0 入力 $\boldsymbol{y}_{[12,12,1]}^{(0)}$
1 畳み込み層1 $\boldsymbol{a}_{[10,10,10]}^{(1)}, a_{cij}^{(1)} = \sum_c \sum_s \sum_t w_{st}^{(1)} y_{(s+i)(t+j)}^{(0)}$
2 活性化層1 $\boldsymbol{y}_{[10,10,10]}^{(2)} = \max (\boldsymbol{a}^{(1)}, 0)$
3 プーリング層1 $\boldsymbol{y}_{[5,5,10]}^{(3)} = \max (y_{(li+s)(lj+t)}^{(2)}), \; {\scriptsize \mathrm{where} : s,t\in : \mid 0,l\mid}$
4 畳み込み層2 $\boldsymbol{a}_{[4,4,20]}^{(4)}, a_{cij}^{(4)} = \sum_c \sum_s \sum_t w_{st}^{(4)} y_{(s+i)(t+j)}^{(3)}$
5 活性化層2 $\boldsymbol{y}_{[4,4,20]}^{(5)} = \max (\boldsymbol{a}^{(4)}, 0)$
6 プーリング層2 $\boldsymbol{y}_{[2,2,20]}^{(6)} = \max (y_{(li+s)(lj+t)}^{(5)}), \; {\scriptsize \mathrm{where} : s,t\in : \mid 0,l\mid}$
7 平坦化層 $\boldsymbol{y}_{[80]}^{(7)} \leftarrow \boldsymbol{y}^{(6)}$ $\boldsymbol{\delta}_{[80]}^{(7)} = {}^t\boldsymbol{W}_{[20,80]}^{(8)} \cdot \boldsymbol{\delta}^{(8)}$ -
8 全結合層1 $\boldsymbol{a}_{[20]}^{(8)} = \boldsymbol{W}_{[80,20]}^{(8)} \cdot \boldsymbol{y}^{(7)} + \boldsymbol{b}_{[20]}^{(8)}$ $\boldsymbol{\delta}_{[20]}^{(8)} = f'(\boldsymbol{a}^{(8)}) \odot \boldsymbol{\delta}^{(9)}, \; {\scriptsize where : f = \max}$ $\Delta w_{kl}^{(8)} = \delta_k^{(8)} \cdot y_l^{(7)},\; b_k^{(8)} = \delta_k^{(8)}$
9 活性化層3 $\boldsymbol{y}_{[20]}^{(9)} = \max (\boldsymbol{a}^{(8)}, 0)$ $\boldsymbol{\delta}_{[20]}^{(9)} = {}^t\boldsymbol{W}_{[3,20]}^{(10)} \cdot \boldsymbol{\delta}^{(10)}$ -
10 全結合層2 $\boldsymbol{a}_{[3]}^{(10)} = \boldsymbol{W}_{[20,3]}^{(10)} \cdot \boldsymbol{y}^{(9)} + \boldsymbol{b}_{[3]}^{(10)}$ $\boldsymbol{\delta}_{[3]}^{(10)} = \boldsymbol{y}^{(11)} - \boldsymbol{c}_{[3]}$ $\Delta w_{ij}^{(10)} = \delta_i^{(10)} \cdot y_j^{(9)},\; b_i^{(10)} = \delta_i^{(10)}$
11 活性化層4(出力) $\boldsymbol{y}_{[3]}^{(11)} = \mathrm{softmax} (\boldsymbol{a}^{(10)})$ - -

※ なぜNo.10:逆伝播で $\boldsymbol{\delta} = \boldsymbol{y} - \boldsymbol{c}$ となるのかはそっち系のサイト(こことか)や書籍を読んでほしい。

誤差 $\delta$ が各層で計算され、次の層へ伝わっていくイメージが伝わるだろうか? これをすべての層に対して遡っていけば、めでたく一個の教師データに対する 学習が完了するわけだ。

実際には複数個(N)の教師データについて $\Delta \boldsymbol{W}$ を求め、それを平均して $\boldsymbol{W}$ の更新に使う。数式(2)は実際は次のような感じになる。 $\boldsymbol{W}$ の右下の $i$ はi番目の教師データから得られたものという意味。

{ \displaystyle
\boldsymbol{W} \leftarrow \boldsymbol{W} - \eta \frac{1}{N} \sum_{i=1}^{N} \Delta \boldsymbol{W}_i \qquad (5)
}

1-2. 誤差逆伝播の実装

面倒な理論・数式の説明はこれぐらいにして、どう実装したか説明しよう。

処理の全体像

下記のmain関数内で呼ばれているtrainLoopが中心部だ。

(Main-train.hs)

main = do
     :
  layers' <- trainLoop getTeachers sampleE putF (learnR st) (layers st) is
     :

学習を繰り返すためにepoch番号を入れた配列isを渡していて、学習済みの状態 (層リスト)を次に引き渡すためちょっと面倒だが再帰になっている。

(Main-train.hs)

trainLoop :: (Int -> IO [Trainer]) -> [Trainer]
  -> (Int -> [Layer] -> [Trainer] -> IO ()) -> Double -> [Layer] -> [Int]
  -> IO [Layer]
trainLoop _ _ _ _ ls []             = return ls
trainLoop getT se putF lr ls (i:is) = do
  ls' <- trainLoop getT se putF lr ls is
  teachers <- getT i
  let
    rls = tail $ map reverseLayer $ reverse ls'    -- fist element isn't used
    dls = map (train ls' rls) teachers
    ls'' = update lr (transpose dls) ls'           -- dls = diff of layers
  putF i ls'' se
  return ls''

trainLoop関数の最初の行でisを一つ減らしてtrainLoopを呼んでいるところ。 このため、isの要素がなくなるところまでtrainLoopの呼び出しが続き、 最後に「学習前の層リスト」が返され、以後少しずつ層リストを学習で修正する。 最後に、学習が完了した層リストがmain中のlayers'に入るという寸法だ。

trainLoopの中では、上記の通り一つ前のepoch番号の学習結果を得てから 今回のepochで使う教師データを生成している。yusugomori実装では教師データ数が 150個に対しバッチ数が150のため毎回全件で学習しているのだ。

次のrlsはちょっとややこしい。学習では誤差逆伝播法を使うが、 そのため層リストを逆順に辿りたい。ただ単純な順序の反転ではダメなので、 まず層リストを反転した上で各層を「逆変換」している(reverseLayer)。 全結合層ではフィルタ行列の転置が逆変換に相当する。

次の行に出てくるtrainがメイン処理だ。これは後述する。 この行で教師データ全部に対し勾配($\Delta \boldsymbol{W}_i$)を計算し、 次の行で層リストを更新して終わりだ。なおputFは学習結果を使って テストデータの認識精度を計算し、表示する処理だ。5 epoch毎に結果を出力する。

次にtrainについて。

(Trainer.hs)

train :: [Layer] -> [Layer] -> Trainer -> [Layer]
train [] _ (i, c) = []
train ls rls (i, c) = dls
  where
    (y, op') = splitAt 1 $ forwardProp ls [i]
    d = (head $ head $ head y) `vsub` c
    (_, dls) = backwardProp (zip (tail op') rls) (d, [])

where句の次の行は前回までに作成した順伝播処理だ。得られた 結果リスト(各層の処理結果)を最後の出力=yとそれ以外に分離している。 次の行、yがImage型のため余計な次元を除去して一次元ベクトルにした後で $\boldsymbol{y} - \boldsymbol{c}$ を計算している (にしても3個もheadがあるのは格好悪い)。 これが最初の誤差 $\boldsymbol{\delta}^{(10)}$ だ。

3行目が逆伝播の本体。一つ目の引数zip (tail op') rlsは 各層の出力値と逆変換用の各層を層毎に組みにしたもののリストだ。 逆伝播はこのリストに対して繰り返す。 二つ目の引数で前層で計算した誤差と各層の勾配、それを更新して次へ渡す。

(Trainer.hs)

backwardProp :: [(Image, Layer)] -> (Delta, [Layer]) -> (Delta, [Layer])
backwardProp [] (_, ls) = ([], ls)
backwardProp ((im,l):ols) (d, ls) = backwardProp ols (d', l':ls)
  where
    (d', l') = backwardLayer l im d

ここまでが、先述の表と図で説明したものを処理している部分である。 うーん、今回も話が長くなりそうだ。

誤差と勾配の計算

逆伝播処理は層の種類で異なるため、backwardLayer関数のパターンマッチで 仕分けている。今のところ活性化層(ActLayer)と全結合層(FullConnLayer)のみ 実装している。他は次回になんとかしようと思う。

(Layer.hs)

backwardLayer :: Layer -> Image -> Delta -> (Delta, Layer)
backwardLayer (NopLayer)           _  d = (d, NopLayer)
backwardLayer (ActLayer f)         im d = deactivate f    im d
backwardLayer (MaxPoolLayer s)     im d = depoolMax  s    im d
backwardLayer (ConvLayer s fs)     im d = deconvolve s fs im d
backwardLayer (FullConnLayer fs)   im d = deconnect    fs im d
backwardLayer l@(FlattenLayer x y) im d = (head $ head $ unflatten x y [[d]], l)

中身を示そう。まず簡単な活性化層の方から。

(ActLayer.hs)

deactivate :: ActFunc -> Image -> Delta -> (Delta, Layer)
deactivate f im delta
  | c == [0.0]  = (zipWith (*) delta f', ActLayer relu')
  | c == [1.0]  = ([], ActLayer f)
  | otherwise = ([], ActLayer f)
  where
    c = f [0.0]
    f' = relu' (head $ head im)

出力に近い方、全結合層の後の活性化層のうち最後の層=出力層では softmax関数を使っているが、誤差を計算する必要がないのでほったらかしている。 一方その二つ前の活性化層3(No.9)はReLUを採用しているので、逆関数 $f'$ は ステップ関数(relu'として別に定義)になり、前層からきた誤差と 掛け合わせている(zipWith (*) delta f'の部分)。

ところで、活性化関数の種類によって処理を変えたいのだが、普通の引数の 型(DoubleとかStringとか)や値でのパターンマッチが使えないみたいだ。 引数fの値(reluとかsoftmaxとか)で分けられたら楽なんだが。。。 それができないので、関数fにダミー値([0.0])を与えた結果で場合分け している。reluなら0.0に、softmaxなら1.0になるからだ(たぶん)。

次は全結合層の説明。

(FullConnLayer.hs)

deconnect :: [FilterF] -> Image -> Delta -> (Delta, Layer)
deconnect fs im delta = (mmul delta fs', FullConnLayer $ calcDiff delta im')
  where
    fs' = tail fs
    im' = head $ head im

calcDiff :: Delta -> [Double] -> [FilterF]
calcDiff delta im = map (mulImage im') delta
  where
    im' = 1.0:im
    mulImage :: [Double] -> Double -> [Double]
    mulImage im d = map (*d) im'

誤差は、フィルタfsと前層からの誤差deltaを掛けている(mmul delta fs')。 fsはすでに転置行列にしてあるが、その最初の要素はバイアスbに相当するため 削る(tail fs)。全結合層は学習によりフィルタを更新するから、誤差から 学習に必要な勾配 $\Delta W$を計算している。calcDiffがそれ。 誤差値deltaと順伝播での入力値imのそれぞれの要素を掛け合わせて 作り出している。本実装ではバイアスbはフィルタに統合されているので、 入力値に1.0を加えているが(im' = 1.0:im)。

全結合層のフィルタ更新

さあ最後に、各教師データから得られた勾配を平均してフィルタを更新する 部分を見てみよう。

(Trainer.hs)

update :: Double -> [[Layer]] -> [Layer] -> [Layer]
update lr _ [] = []
update lr [] ls = ls
update lr (dl:dls) (l:ls) = (updateLayer lr l dl):(update lr dls ls)

dlは処理対象の層lに対する各教師データからの勾配が入っている。 これをupdateLayerに渡して計算している。lrは学習率だ。

次が、全結合層での勾配を集計する処理である。

(FullConnLayer.hs)

updateFullConnFilter :: [FilterF] -> Double -> [Layer] -> Layer
updateFullConnFilter fs lr dl = FullConnLayer fs'
  where
    ms = strip dl
    delta = mscale (lr / (fromIntegral $ length ms))  $ msum ms
    fs' = msub fs delta

ほぼ数式(5)そのままなので、ほぼ説明不要だろう。ms = strip dlに ついてだけ。各教師データから得た勾配のリストから、この層での対象外の データがあれば除外するようにしている。普通は問題ないはずだが。

2. 「学習」の評価

プログラムも出来上がったことだし、どれだけ学習するか評価してみよう。 コンパイルと実行は次の通り。

% cabal build
   :
% ./dist/build/train/train
Initializing...
Training the model...
iter =     0/500 accuracy = 0.3333333333 time = 0.149965s
iter =     5/500 accuracy = 0.3333413554 time = 5.132365s
   :

途中経過を 5 epoch 毎に画面に出している。前回から少しフォーマットを 変更し、正答率を「認識精度=accuracy」に改めている。

2-1. 学習結果(ノーマル版)

上記のプログラムを10回実行し、それぞれ認識精度と実行時間を計測した。 結果は次のグラフの通り。

f:id:eijian:20161108015850p:plain

点線が各試行時の実測値、太い実線が10回の平均値だ。これを見ると、認識精度は 各試行でかなりばらついているのがわかる。また、認識精度の向上が 急峻なものと緩慢なものがあり、かなり差がある。

プログラムでは、学習用データ、テスト用データ、および各フィルタに乱数を 使っているが、それらの初期状態によりこれだけの差が生じてしまっていると 思われる。(学習効果が薄い=うまく特徴が出ない画像群が生成されたとか、 フィルタの初期状態がイケてないとか)

ただ、総じて 500 epoch ぐらい繰り返せば認識精度がほぼ9割を超えることも 分かった。 まだ畳み込み層の学習をしていない段階でここまで精度が上がるとは思っておらず、 嬉しい誤算だった。(こんなものなのかな?)

2-2. 学習結果(学習率の差)

先のグラフで、認識精度の向上は最初からではなく50 epochを超えたあたりから 急に上がりだしている。私は素人なのでよくわからないが、いろいろな情報を見ると 「学習率」はどう設定するかとか、学習進行に応じて変化させるとか書かれて いて、重要なパラメータであるらしい。

このプログラムでは学習率を0.1に固定しているが、初期段階ではこれが小さすぎる のかもしれない。ということで、学習率を0.2と0.05に変えて試行した結果も出して みた。(0.2と0.05の試行は5回、その平均をとった)

f:id:eijian:20161108015855p:plain

予想通りかというとちょっと微妙だが、学習率を大きくした方がより速く認識精度が 向上していることがわかる。ただ、初っ端から急上昇しているわけではないので もっと大きくしないといけないのだろう。ただ、大きいままで学習を続けると 収束しないという話もある(?)ので、やはり動的に変化させるのがいいのだろう。 あとでそれにも取り組みたい。

今回のプログラムはyusugomori実装を「写経」することが第一課題なので 学習率は0.1に戻しておこう。

2-3. yusugomori実装との比較

ここまで、学習が進んだ!と喜んできたわけだが、「写経」元のyusugomori実装と 同じことができているのだろうか?

ということでyusugomori実装との比較をしてみたのが次のグラフだ。

f:id:eijian:20161108015900p:plain

avg-0.1が本プログラム(学習率0.1、試行10回の平均値、青線)、py-cnnが yusugomori実装(学習率0.1、試行1回、赤線)、py-cnnfcがyusugomori実装で 畳み込み層の学習を抑制した改造版(学習率0.1、試行1回、緑線)だ。 yusugomori実装では乱数のseedが固定されているため毎回同じ結果になるのと、 後で述べるようにめちゃくちゃ処理時間が掛かるのとで試行は1回だけにした。

先述の通り乱数の生成値により試行結果はかなりぶれるのだが、青と緑の線は 350 epoch近辺からほとんど同じような値になっていて、(楽観的に見れば) 本プログラムはyusugomori実装と同等な処理ができているのではないかと思う。 一方本来の学習(=畳み込み層も学習)であれば認識精度が急上昇し、かつ 200 epoch ぐらいでも認識精度が98%を超えてくる。それに比べ本プログラムでは 畳み込み層のフィルタがいつまでたっても初期状態(ランダム)のままなので、 この差は致し方ないところだろう。まあ、ここまで性能差があると早く 畳み込み層の学習部分も完成させたいと思うからモチベーションの向上になる。

次に処理時間について見てみる。200 epoch 実行時の結果は次表のとおりである。 (iMac early 2009, 3.0GHz Core 2 Duo, 4GB memory)

avg-0.1 py-cnn py-cnnfc
処理時間(s) 157 24838 6152
処理性能(epoch/s) 1.273 0.008 0.033
認識精度(%) 74.2 98.4 62.9

yusugomori実装の処理時間が半端ではないことがわかる。本プログラムも コンパイルされたネイティブコードであることを考えると決して十分な処理速度 ではないが、それでもまったく比べものになっていない。

NumPyライブラリのインストールに何か抜けている手順があるのだろうか? インストール時に何か特別なことをした記憶もないので、以下の記事にあるように 「ちゃんとした」インストールができていないのかもしれない。

(参考1) Kesin's diary: NumPy, SciPyをちゃんとインストールする

もしくは次の記事のように多重ループと配列要素への直接読み書きが 影響しているのかもしれない。

(参考2) numpyで明示的にループを書くと極端に遅くなる

いずれにせよ、当方は今の所Pythonを使う気はないので深入りはしない...

3. まとめ

今回は全結合層部分の逆伝播処理を実装した。これにより学習ができるように なったので、精度はともかく画像認識ができるようになった。ただ、これはまだ 半分。次回は畳み込み層の逆伝播を実装し、プログラムを完成させようと思う。

DeepLearning(2): まずは順伝播(下)

前回に続き、CNNの順伝播処理について。前回は下準備という感じだったが、 今回は各層の処理を考よう。

なお今回の記事に限らないが、私の記事は自分自身の忘備録であることが 主目的なため、各種技術・理論については正しくないかもしれないし プログラミングに必要なことに絞って書いている。よってCNNやDeep Learningの 正しいところは前回紹介したような、各種書籍・Webサイトにあたってほしい。

CNNの組み立て

まずはyusugomoriさんの Python実装 (yusugomori実装)を写経するので、そこで定義されているCNNを作りたいと思う。

CNNの構造

前回出した図をもう一度使うが、作成するCNNの構造は次のようになっている。

f:id:eijian:20160830205123p:plain

  • 第1畳み込み層
    • 入力画像はモノクロ12x12ドット(1チャネル)
    • 畳み込みフィルタは3x3ドット、10枚
    • 畳み込みにより10x10ドット、10チャネルの「画像」に変換される
  • 第1活性化層
    • 活性化関数はランプ関数(ReLU)、全ドットにReLUを適用
    • 画像サイズ、チャネル数に変化なし
  • 第1プーリング層
    • プーリングサイズは2x2ドット
    • 最大プーリングを適用
    • 出力画像の大きさは10÷2で5x5ドット、10チャネルとなる
  • 第2畳み込み層
    • 入力画像はモノクロ10x10ドット(10チャネル)
    • 畳み込みフィルタは2x2ドット、20枚
    • 畳み込みにより4x4ドット、20チャネルの「画像」に変換される
  • 第2活性化層
    • 活性化関数はReLU、全ドットにReLUを適用
    • 画像サイズ、チャネル数に変化なし
  • 第2プーリング層
    • プーリングサイズは2x2ドット
    • 最大プーリングを適用
    • 出力画像の大きさは4÷2で2x2ドット、20チャネルとなる
  • 一次元化
    • 2x2ドット、20チャネルを一次元化して2x2x20=80個のリストに変換
  • 隠れ層(全結合層)
    • 入力要素は80個、出力要素は20個
  • 第3活性化層
    • 全要素にReLUを適用
  • 出力層(全結合層)
    • 入力要素は20個、出力要素は3個
  • 第4活性化層
    • 全要素にSoftmaxを適用

各層は前回説明したようにmain関数内にLayerのリストとして 定義している。

トレーニング

各教師データ(画像)に対して上記のCNNを適用するのだが、一枚毎に 次のtrain関数で変換している。前回示したloop関数内で呼び出されている。 中身は各層の順伝播処理を呼んでいるだけなので簡単だ。

(Trainer.hs)

train :: [Layer] -> Trainer -> IO (Image, [Layer])
train [] (i, c) = return (i, [])
train ls (i, c) = do
  let op  = forwardProp ls [i]
      ls' = backwordProp ls op
  return (head op, ls')

forwardPropが順伝播、backwordPropが逆伝播だが、まだ順伝播のみの実装だ。 forwardProp再帰により層を順に適用し、それぞれの順伝播処理を呼び出す ようにした。

(Layer.hs)

forwardProp :: [Layer] -> [Image] -> [Image]
forwardProp [] is = is
forwardProp (l:ls) is = forwardProp ls (forward l is)

forward :: Layer -> [Image] -> [Image]
forward _ [] = []
forward NopLayer i = i
forward l@(ActLayer f) im@(i:is)     = (activate l i):im
forward l@(MaxPoolLayer s) im@(i:is) = (poolMax l i):im
forward l@(ConvLayer s fs) im@(i:is) = (convolve s fs i):im
forward l@(HiddenLayer fs) im@(i:is) = (connect fs i):im
forward l@(FlattenLayer f) im@(i:is) = (f i):im

各層の処理

何もしない層(NopLayer)

入力をそのまま返すダミーだ。(プログラムの動作を確認するため)

畳み込み層(ConvLayer)

この層はCNNのキモになるので自分の忘備録のために少し丁寧に記述しておく。 畳み込み層では、画像の各チャネルについてフィルタを適用する。 フィルタの大きさは3×3などの比較的小さいものだ。一般的でないかもしれないが 簡便のため正方形とする。

各チャネルの二次元データ上の各画素について、そこを起点にフィルタと 同じ大きさの正方形を取り出してフィルタの対応する画素と掛け合わせて 合計するのだ。下図では入力を$x$、出力を$u$、フィルタを$w$とする。

f:id:eijian:20160910202832p:plain

この処理、手続き型言語であれば多重ループで難なく書ける処理なのだが、 Haskellでどう書くか? 画像をリストでなくArrayで表現して、各画素を何度も順序ばらばらにアクセス すればすれば多重ループで書けるが、今回は再帰で処理するために次のようにした。

f:id:eijian:20160910202840p:plain

フィルタの大きさをsとすると、画像データの最初のs行をtakeで取り出して 各行の処理(convolveLine)を行う。そのあと、1行目を取り除き(tail)、 また同じことを繰り返すのだ。残りの行数がフィルタの大きさ(s)より小さく なったら終了する。

(ConvLayer.hs)

-- s: filter size

convolvePlain :: Int -> [Double] -> Plain -> Plain
convolvePlain s k ps
  | len < s   = []
  | otherwise = (convolveLine s k p:convolvePlain s k ps')
  where
    len = length ps
    p   = take s ps
    ps' = tail ps

各行についてtakeで行の最初からs個取り出すと$s \times s$個の 画素データが得られる。これにフィルタを適用する。 そのあと、各行の最初の要素を取り除き、それを 使って同じ処理を行の要素がs個未満になるまで繰り返す。 sum $ zipWith (*)の部分がフィルタを適用しているところである。

(ConvLayer.hs)

convolveLine :: Int -> [Double] -> [[Double]] -> [Double]
convolveLine s k ps
  | len < s   = []
  | otherwise = ((sum $ zipWith (*) vs k):convolveLine s k ps')
  where
    len = minimum $ map (length) ps
    vs = concat $ map (take s) ps
    ps' = map (tail) ps

この場合、入力画像のサイズに対し「フィルタのサイズ-1」だけ 小さな画像が出力される。今回の例では、元画像が12x12、フィルタが3x3のため、 12-(3-1)=10となり10x10の画像が得られる。

もう一つ、フィルタサイズより小さくなった時は黒や白のドットを 補うことで入力画像のサイズを維持する手法もあるようだ。 が、一般的な画像認識で使われるであろう数百×数百程度の画像では誤差の ようなものなので前者の手法を使う。

こうすれば、画像の各画素に対して再帰的にフィルタを適用できる。 入力画像の1チャネル毎にこの処理を繰り返すと、全体としては下図のように 複数チャネルの画像に変換されて出てくるわけだ。

f:id:eijian:20160910202845p:plain

その後、入力画像の各チャネルの処理結果を足し合わせる必要がある。 これは各チャネルの同じ位置にある画素データを足すだけなので簡単な ことのはずだが、頭がこんがらがって以外と手間取った。

(ConvLayer.hs)

sumPlains :: Double -> [Plain] -> Plain
sumPlains b [] = repeat $ repeat b  -- 2 Dimensions (X*Y)
sumPlains b (p:ps) = zipWith f p (sumPlains b ps)
  where
    f :: [Double] -> [Double] -> [Double]
    f [] _          = []
    f (a:[]) (b:_) = [a + b]
    f (a:as) (b:bs) = (a + b):(f as bs)

特に全チャネル分を足し合わせ、さらにバイアス値を加える必要があるので 多少複雑になってしまった。全画素に同じバイアス値を加えるために、 repeat $ repeat バイアス値により画像サイズがわからなくても 処理できた。普通ならちゃんと画像サイズ分の二次元データを用意する ところだろう。細かいところだが、遅延評価の恩恵か。

活性化層/活性化(ActLayer)

先述の通り、ニューラルネットについての書籍やサイトでは活性化層という 独立した層として説明されておらず、畳み込み層や隠れ層の後処理的な扱いだ。 が、今回は独立した層としたほうが「私にとって」わかり易そうだったので 分けることにした。 ただ、後でチューニングのために他の層に統合するかもしれない・・・。

活性化関数として、ランプ関数(Rectified Linear Unit, ReLU)とsoftmax関数を 用意した。それぞれ下記で表される。 これらはActLayer層の作成時にどちらか選ぶのだ。

(ReLU) $f(x_i) = \max(x_i, 0)$

(Softmax) $f(\boldsymbol{x}) = \frac{e^{x_i}}{\sum_{k=1}^K e^{x_k}}$

ここで、$x_i$は$i$番目の出力値、$K$は出力値の個数(次元)である。 実装はこう。

(ActLayer.hs)

relu :: ActFunc
relu as = map f as
  where
    f :: Double -> Double
    f a = if a > 0.0 then a else 0.0

softmax :: ActFunc
softmax as = map (\x -> x / sume) es
  where
    amax = maximum as
    es   = map (\x -> exp (x - amax)) as
    sume = sum es

記事を書いていて、ReLUのほうは素直にmaxを使えばいいと思った。 後で直そう。Softmaxの方も、ほぼ定義通りだ。ただし入力xに対し、 xの最大値を引いているがこれは発散しないための工夫だそうだ。 activate本体はすべての画素に対して活性化関数を適用する。 まあ、難しいものではない。

(ActLayer.hs)

activate :: Layer -> Image -> Image
activate (ActLayer f) im = map (map f) im

プーリング層(MaxPoolLayer)

本プログラムでは最大プーリングのみ対応する。他も必要になったらその時考えよう。 プーリングの処理では画像を小ブロックに分け、その中の各値から出力値を 計算する。最大プーリングでは、最大値を選択して出力するわけだ。 幾つかの画素から一つだけ値を出力するのだから、当然出力値の数は少なくなる。 小ブロックが2x2ドットだとすると、元の画像が縦横半分のサイズになるのだ。

f:id:eijian:20160910202851p:plain

これをプログラムにしてみよう。まず、1チャネルの画像の処理について。

(PoolLayer.hs)

poolMaxPlain :: Int -> Plain -> [[Pix]]
poolMaxPlain s p = map (poolMaxLine s) ls
  where
    ls = splitPlain s p

splitPlain :: Int -> Plain -> [[[Double]]]
splitPlain s [] = []
splitPlain s p  = (p':splitPlain s ps)
  where
    (p', ps)  = splitAt s p

poolMaxPlainが1チャネル分の処理の主体だ。入力データ(Plain)を splitPlainによってブロックサイズで縦方向に分割している。 縦sドットの横に細長いデータをそれぞれpoolMaxLineに引き渡して 処理させているのだ。

(PoolLayer.hs)

poolMaxLine :: Int -> [[Double]] -> [Pix]
poolMaxLine _ [] = []
poolMaxLine s ls
  | len == 0  = []
  | otherwise = (pixs:poolMaxLine s ls')
  where
    len  = length $ head ls
    pixs = max' $ zip (concat $ map (take s) ls) [1..]
    ls'  = map (drop s) ls

max' :: [Pix] -> Pix
max' [] = error "empty list!"
max' [x] = x
max' (x:xs) = maximum' x (max' xs)

maximum' :: Pix -> Pix -> Pix
maximum' a@(v1, i1) b@(v2, i2) = if v1 < v2 then b else a

ここは(今は)余計な処理が入っているためちょっと複雑だ。説明しよう。 大まかには、与えられた細長いデータからsドット分を取り出し、max'関数で 最大値を選択している。残りのデータは再帰により同じことを繰り返す。 ここで処理対象がDoubleではなくPixという型に変わっているが、 その定義はこう。

(PoolLayer.hs)

type Pix = (Double, Int)

これは小ブロックの各要素に番号をつけるために用意した型だ。 poolMaxLineでは、$s \times s$ドット分取り出したら一次元化して それに番号を振っている。以下の部分。

    zip (concat $ map (take s) ls) [1..]

これは最大プーリングでどの画素が最大だったかを覚えておいて、 後々逆伝播処理で伝播先(?)とするために必要だとか (とどこかで読んだはず、だがまだ勉強中で本当かどうかわかっていない。 いつのまにかしれっと消えているかもしれない。。。)。

二次元データについてpoolMaxPlainで処理できれば、それを画像データ全部に 適用すればよいだけだ(poolMaxpoolMax')。

一次元化(FlattenLayer)

これは層というよりデータの整形だ。畳み込み層やプーリング層の 出力は二次元平面、複数チャネルという、いわゆる「画像」の形をしているが、 隠れ層では「一次元配列」を入力に取る。ということでその間に変換が必要なのだ。

f:id:eijian:20160910202856p:plain

(Image.hs)

flatten :: Image -> Image
flatten im = [[concat $ concat im]]

単に全チャネルの画素値を一列に並べているだけである。 ただし、Layer型なので結果は「画像=3次元配列」にする必要がある。 なのでわざわざ前後に2つの角カッコをつけて3元配列にみせかけている。

隠れ層(HiddenLayer)

これはいわゆる全結合層なのだが、他の記事やプログラムを見ると隠れ層と いう名で作られているのでそれに倣った。隠れ層としながら、最終の出力層も このHiddenLayerを使っているので名前は変えるべきかもしれない。。。

この層では、入力(channel)c個を出力(kernel)k個に変換する。 c個の要素全てがk個の各出力に関連付けられているため全結合というらしい。 その関連付けをフィルタとして与える必要があるのでそのサイズは$c \times k$に なる。初期はランダム値を設定し、学習によってその値を少しずつ補正することで 画像認識などができるらしい。フィルタは$c \times k$の二次元リストで 与えることにしよう。

(HiddenLayer.hs)

initFilterH :: Int -> Int -> IO [FilterH]
initFilterH k c = do
  f <- forM [1..k] $ \i -> do
    f' <- initKernel c
    return f'
  return f
  where
    a = 1.0 / fromIntegral c
    initKernel :: Int -> IO FilterH
    initKernel c = do
      w <- forM [1..c] $ \i -> do
        r <- MT.randomIO :: IO Double
        return ((r * 2.0 - 1.0) * a)
      return (0.0:w)

各要素は入力要素数cに対し$1/c$ ~ $-1/c$ の間の乱数としている。 ただし、最初の要素はバイアス値で、初期値は0.0とするらしい。

この層では入力値にフィルタ値を掛けてすべて足しこむ(k個分)。 式で書くとこうなる。なおフィルタの要素は $w_{ij}$ としている。

[u_i = \sum{j=0}^c w{ij} x_j + b_i]

ここで$x_j$は$j$番目の入力値、$u_i$は$i$番目の出力値、$b_i$は$i$番目の バイアス値である。

実装はこう。

(HiddenLayer.hs)

connect :: [FilterH] -> Image -> Image
connect [] _ = error "invalid FilterH"
connect _ [] = error "invalid Image"
connect fs [[im]] = [[map (dot im) fs]]
  where
    dot :: [Double] -> [Double] -> Double
    dot im f = sum $ zipWith (*) (1.0:im) f

dot関数のzipWith第一引数に1.0が追加されているのは先の$b_i$を 足し込むためである(FilterHには$b_i$が含まれている)。

実行してみた

以上でひとまず順伝播については実装できた(と思う)。まずはコンパイル

$ cd deeplearning
$ cabal configure
  :
$ cabal build
  :

では実行してみよう。

$ dist/build/deeplearning/deeplearning
Building the model...
Training the model...
iter = 0/200 time = 0s ratio = 0.0
iter = 5/200 time = 0.0030001s ratio = 0.33333713817332733
iter = 10/200 time = 0.333019s ratio = 0.33333713817332733
iter = 15/200 time = 0.682039s ratio = 0.33333713817332733
     :
iter = 200/200 time = 13.0807482s ratio = 0.33333713817332733
Finished!

反復5回毎に状態を出力しており、左から

  • epoch(反復回数)
  • 経過時間(秒)
  • 正解率
    • 正解率は、全テストデータを順伝播にかけ得られた出力値から正解にあたる 要素の値(その要素が正解である確率)を取り出して平均したもの。
    • 3種類に分類されるデータを使ったので何も学習していない状況だと 出力は均等割合となり約1/3になる。

である。 正解率は期待通り約1/3になっているものの、学習してないのでこの結果だけで はなんとも言えない。ただ、何となくそれらしい出力が出てくると楽しい。

まとめ

前回から取り組んでいたCNN画像認識プログラムの「順伝播」のみ実装完了した。 出力結果にはほとんど意味はないが、各層の処理を経てなんとなく結果が出力 されたので雰囲気はつかめたかなと。

次回以降で本当のメインである学習処理を実装し「意味のある結果」が得られる ことを期待しよう。

(ここまでのソースはこちら

DeepLearning(1): まずは順伝播(上)

さて、今回からは新たなネタとしてDeepLearningに取り組もう。 DeepLearningといえばすでにブームは去って、今は現実の課題に どんどん活用されている状況だ。使いやすく性能のいいツールや ライブラリがいろいろあり、ちょっと学べば簡単にその恩恵に 与ることができる(かもしれない)。

しかし、仕組みがわからないのをただ使うのは面白くないし 個人的に写真の自動分類をしたいというニーズもあるし、遅まきながら 自分で作ってみよう、というわけだ。

目標と参考資料

最終的には「一般物体認識」まで行きたい(というかそうでないと 写真の分類にならない)が、目標が高すぎると頓挫するのでここでは 以下の2つのステップまでとする。

  • ステップ1: yusugomoriさんのpython実装を写経する
  • ステップ2: 上を使ってお決まりのMNISTデータを使った性能検証をする

「自分で作ってみよう」と大口を叩いていながら他の方の実装を写経する とは詐欺的だが、そこは当方素人なので。。。以後、yusugomoriさんのpython実装を 「yusugomori実装」と表記する。

作るのは、畳み込みニューラルネットワーク(Convolutional Neural Network, CNN)による画像認識プログラムとする。以下、参考にした書籍や Webサイトを列挙。

ひとつの本やWebサイトだけでは理解が進まないので、多数に目を通して同じ概念を いろいろな説明で読むほうがよいと思う。

まずは型を考える

(ちゃんと理解できているとは言い難いが)CNNによる画像認識は 次の図のように「画像」を各層で変換していって最終的に分類結果 (どれに該当するかを確率で示す)を得るようだ。

f:id:eijian:20160830205123p:plain

そこでまずは「画像(Image)」型と「層(Layer)」型を定義しよう。 カラー画像だと、$X \times Y$二次元のドットで構成された平面(Plain)が 赤緑青の3色(チャネルとする)集まってできている。一般の画像ファイル では各ドットを0-255の整数で表すことが多いが、CNNでは強さを実数で表すようだ。 畳み込み層の変換後のデータも多チャンネルの「画像」とみなせそうなので、 Imageは次のように定義した。実体はDoubleの3次元リストだ。

(Image.hs)

type Plain = [[Double]]    -- 2D: X x Y pixels
type Image = [Plain]       -- n channel

わざわざImagePlainのリストにしたのは、Plainごとに処理することが 多そうだからだ。 あとついでに、「学習」の際に必要となる教師データも定義しておく。

(Image.hs)

type Class = [Double]    -- trainer vector
type Trainer = (Image, Class)

実数のリストだが、教師データは要素のどれか一つが1.0、他が0.0となる。

次にレイヤを考える。CNNでは次のようなレイヤが使用される。

  • 畳み込み層 (convolution layer)
  • プーリング層 (pooling layer)
  • 活性化層/活性化関数 (activation layer)
  • 全結合層 (fully connected layer)

このうち、活性化層というのは一般的な表現かどうかよくわからない。 最初は型クラスとしてLayerを定義し各層を型としてやればいいと考えた。 が、各層をつなげて一連の変換処理を表現するのに「リスト」を使いたかったので (私の知識範囲では)うまい方法が思いつかなかった。 ということで、Layer型を次のように代数データ型で定義することにした。

(LayerType.hs)

data Layer = NopLayer
           | ConvLayer Int [FilterC]
           | MaxPoolLayer Int
           | ActLayer ActFunc
           | HiddenLayer [FilterH]
           | FlattenLayer (Image -> Image)

NopLayerは何も変換しないダミーの層(テスト用)、 プーリング層は実用上Max Poolingだけでよいと考えてMaxPoolLayer、 全結合層は隠れ層ということでHiddenLayerとしている。 各行に出てくる引数の型はそれぞれ以下のように定義してある。 ActFuncは活性化関数、FilterCは畳み込み層のためのフィルタ、 FilterHは全結合層(隠れ層)のためのフィルタ、をそれぞれ表す型だ。

(LayerType.hs)

-- for Activation Layer
type ActFunc = [Double] -> [Double]

-- filter for Convolution Layer
type Kernel = [[Double]]
type Bias   = Double
type FilterC = (Kernel, Bias)

-- filter for Hidden Layer
type FilterH = [Double]

学習用プログラム

deep learningで画像認識させるためのステップは大雑把には以下の手順を踏む 必要があると思っている。

  1. 大量のサンプル画像を用意し、それらを分類する。
  2. サンプル画像から「教師データ」を作成する。
  3. 教師データを使って「学習」を行う。(学習とは、各層で使われるフィルタを 教師データを元に調整していくこと、と理解)
  4. 画像を調整済みフィルタを使って分類する(各分類の確率を計算する)。

1と2はプログラム以前の準備段階だ。が、yusugomori実装では教師データと 学習状況の確認に使うテストデータをプログラム内で生成している。 最終的には学習データとテストデータは外から与える形にしたいが、 まずは同じようにプログラム内で生成するようにしよう。

教師データの生成

ゆくゆく教師データを外から与えることも考え、「画像データの集合」 の型??Poolを作ろう。これらの集合から幾つかの画像を取り出して 学習に使いたい。そこで、データに追番をつけてMapに登録 することにした。とりあえずはオンメモリで処理できれば良いので 以下のようにプール(オンメモリ版、MemPool)を定義した。

(Pool.hs)

class Pool p where
  getImages :: p -> Int -> Int -> IO [Trainer]
  nSample :: p -> Int

newtype MemPool = MemPool { m :: Map.Map Int Trainer }

yusugomori実装に倣って教師データを生成する関数も用意した。

(Pool.hs)

initSamplePool :: Int -> (Int, Int) -> Int -> Double -> Int -> IO MemPool
initSamplePool c (sx, sy) o p n = do
  s0 <- forM [0..(n-1)] $ \i -> do
    let cl = i `mod` o  -- class of this image

    -- Image data
    s1 <- forM [1..c] $ \j -> do
      s2 <- forM [0..(sy-1)] $ \y -> do
        let p' = if y `div` st == cl then p else (1-p)
        s3 <- forM [1..sx] $ \x -> do
          a <- pixel p'
          return a
        return s3
      return s2

    -- Trainer data
    e1 <- forM [0..(o-1)] $ \j -> do
      return $ if j == cl then 1.0 else 0.0
    return (s1, e1)

  return (MemPool (Map.fromList $ zip [0..] s0))
  where
    st = sy `div` o
    pixel :: Double -> IO Double
    pixel p = do
      v <- MT.randomIO :: IO Double
      let v' = if v < p then 0.5 else 0.0
      return v'

結構ごちゃごちゃしているが、以下のような3種類の画像データを生成してプール として返しているのだ。

f:id:eijian:20160830011222p:plain

実際は上記のようなキッチリした画像ではなく、灰色部分に一定割合で白が 混ざっていて、白色部分は逆に灰色が混ざっている。教師データではその 割合を5%、テストデータは10%としている。

つぎに、プールから指定枚数の画像を取り出す関数getImagesとプール内の 画像総数を返す関数nSampleを示そう。

(Pool.hs)

{-
getImages
  IN : pool
       epoch number
       batch size
-}

instance Pool MemPool where
  getImages p@(MemPool m) e b = do
    let s = nSample p
        o = (e-1) * b `mod` s
        mx0 = o + b - 1
        mx2 = mx0 - s
        mx1 = if mx2 < 0 then mx0 else s - 1
        im0 = mapMaybe (\x -> Map.lookup x m) [o..mx1]
        im1 = mapMaybe (\x -> Map.lookup x m) [0..mx2]
    return (im0 ++ im1)

  nSample (MemPool m) = Map.size m

getImagesは若干ややこしいが、要は教師データを満遍なく学習させたいので epoch毎に違う画像を取り出すようにしている(最初のepochで1番目から10個 取り出したら次のepochでは11番目から10個取り出す)。 ゆくゆくは順番に取り出すのではなくランダムに取り出すオプションも実装したい。

main関数(全体の流れ)

ここで、全体の流れを説明しておきたいのでmain関数の話をしよう。 大雑把には次のような流れだ。

  • パラメータ定義(main以前):フィルタの大きさとかバッチサイズとか (将来的にパラメータを設定ファイルから読み込むようにしたい)
  • 教師データ、テストデータの準備
  • 各レイヤの初期フィルタの生成
  • レイヤの組み立て
  • 学習の繰り返し

まずパラメータの説明から。画像は3種類に分類する(k)。教師データが 50x3=150個、テストデータが10x3=30個だ。画像サイズ、フィルタは 先に示した各レイヤの仕様どおり。なおバッチサイズは教師データが少ないので 毎回全教師データを学習するために150個としている。

(Main-cnn.hs)

-- PARAMETERS

-- training/test dataset
k = 3   -- number of class
n = 50  -- number of teacher data for each class
m = 10  -- number of test data for each class

train_N    = n * k
test_N     = m * k

-- image size
image_size = [12, 12]
channel    = 1

-- filter specs
n_kernels    = [10, 20]
kernel_sizes = [3, 2]
pool_sizes   = [2, 2]
n_hidden     = 20
n_out = k

-- loop parameters
epochs = 200
opStep = 5
epoch0 = 1
batch  = 150

-- others
learning_rate = 0.1

次にメインルーチンだ。教師データとテストデータは次のように生成 している。説明はいらないだろう。

(Main-cnn.hs)

main :: IO ()
main = do
  putStrLn "Building the model..."
  poolT <- initSamplePool 1 (12, 12) 3 0.95 train_N   -- trainer data pool
  poolE <- initSamplePool 1 (12, 12) 3 0.90 test_N    -- test data pool
  sampleE <- getImages poolE 1 test_N

次はレイヤの定義だ。まずフィルタを初期化してからレイヤを複数並べる。

(Main-cnn.hs)

  fc1 <- initFilterC 10 1 12 12 3 2
  fc2 <- initFilterC 20 10 5 5 2 2
  fh1 <- initFilterH n_hidden (2*2*20)
  fh2 <- initFilterH n_out n_hidden

  let layers = [ConvLayer 3 fc1, ActLayer relu, MaxPoolLayer 2,
                ConvLayer 2 fc2, ActLayer relu, MaxPoolLayer 2,
                FlattenLayer flatten,
                HiddenLayer fh1, ActLayer relu,
                HiddenLayer fh2, ActLayer softmax]

ここでは、畳み込み層1(活性化関数ReLU)→プーリング層1→ 畳み込み層2(活性化関数ReLU)→プーリング層2→隠れ層(活性化関数ReLU)→ 出力層(活性化関数softmax) という多層構造にした。 なお、FlattenLayerは畳み込み+プーリングの処理から全結合へデータ形式を 変換する処理である。

最後が学習の繰り返しだ。

(Main-cnn.hs)

  tm0 <- getCurrentTime
  putStatus 0 tm0 [([0.0], 0.0)] layers
  loop is tm0 layers batch poolT sampleE

時間の記録と学習前の状態をダミーで表示したあと、繰り返し学習する loopを呼んでいる。loopの詳細は次の通り。

(Main-cnn.hs)

loop :: Pool p => [Int] -> UTCTime -> [Layer] -> Int -> p -> [Trainer]
     -> IO ()
loop [] _ _ _ _ _ = putStr ""
loop (i:is) tm0 ls b pt se = do
  teachers <- getImages pt i b
  ops <- mapM (train ls) teachers
  let (_, dls) = unzip ops
      ls' = updateLayer ls dls   -- dls = diff of layers
  if i `mod` opStep == 0 then putStatus i tm0 (evaluate ls' se) ls'
                         else putStr ""
  loop is tm0 ls' b pt se

大したことはやってなくて、epoch毎に繰り返しloopを呼び、 epoch番号がなくなったら(空リストになったら)、繰り返しを終了する。 中身は、学習する教師データをプールから取り出して、それぞれ学習をし、 各層を学習結果を使って更新する。とはいえ、まだ「学習」は実装しておらず、 この部分はダミーである。

ここまでがmain関数だ。

長くなってしまったので、各層(レイヤ)の処理については次に 回すことにする。

まとめ

今回はCNNの画像認識プログラムの初回として以下に取り組んだ。

  • CNNで使われるいろいろな型を定義した。
  • 入力値を変換するいろいろな「層」を定義した。
  • これらを使って順伝播の処理を実装した。

次回は各層についての説明。

(ここまでのソースはこちら