Haskellでいってみよう

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

レイトレーシング(12): 集光模様!

0. 以前の記事の訂正

第8回の記事で完全拡散反射の場合のBRDF($f_r$)を下記のように示した。

 \displaystyle
  f_r = \frac{1}{2\pi}

が、間違いだった。shocker-0x15 さんからご指摘をいただきました。 ありがとうございます。正しい式は以下。

 \displaystyle
  f_r = \frac{1}{\pi}

これにより、得られる輝度が全体的にほぼ倍になるが、輝度を画面の 色に変換するところで調整してやればいいので大きな問題はなかった。 やれやれ。

1. 今回の改善

今回は例トレーシング回の最終回ということで、やっと目標だった集光 模様の実現に取り組むことにする。また、ついでに球以外の形状の サポート、設定ファイルの対応もやってみた。

-1. Material型

まず最初にMaterial型の拡張について説明したい。Material型は 物体の物性や模様を表すために定義した型だ。前回までは単純な 完全拡散面を想定していたためパラメータは少ししかなかった(コード 自体には書いていても使っていないものもあった)。鏡面反射・屈折を サポートするにあたりいくつかのパラメータを追加したし、実際に 値を使うようにコードを拡張した。下表に従来と今回の対応を示す。

parameter 説明 前回版 今回
emittance Radiance 自己発光強度
reflectance Color 拡散反射率(≒物体色?)
specularRefl Color 鏡面反射率(金属反射) -
transmittance Color 透過率 - -
ior Color 屈折率 -
diffuseness Double 拡散度※1
metalness Double 金属性※2 -
smoothness Double 平滑度 - -

※1: 拡散反射と鏡面反射のあり合いを表す。1.0だと完全拡散反射、0.0 だと完全鏡面反射となる。鏡面反射方向からの光をどれだけ反射するか。

※2: 金属性と書いたが金属反射と透過の割合を表す。1.0だと完全に 不透明で金属反射のみ、0.0だとガラスのような透明物体となる。

なお、まだ対応できていないパラメータもある。transmittanceは透明物質内での 光の通過割合で、色ガラスのような物体を表すため、smoothnessは表面のざらつき 具合でぼやけたハイライトなどを表すためのパラメータだ。ぜひともこれらも実装 したいが、今後の課題とする。

iorの型がColorになっているのは、色(光の波長)毎に屈折率が違う 場合があるため。三角プリズムによる光の分散を再現するために必要だ。 (とはいえ、現状は赤、緑、青の3波長しか扱っていないため、いわゆる虹色に ならないが)

各パラメータがどのように使われるかは後の節で説明する。

-2. 鏡面反射・屈折ベクトル

反射・屈折ベクトルの求め方については、いろいろなところに情報があるので ここでは割愛し、数式とコードのみ示すことにしよう。各記号・文字がわかり 易いように、全体を図にするとこんな感じ。

f:id:eijian:20180707022753p:plain

(反射)

入射ベクトル  \boldsymbol{e} 、法線ベクトル $\boldsymbol{n}$ とする時、求める反射ベクトル $\boldsymbol{v_r}$ は、

 \displaystyle
\boldsymbol{v_r} = \boldsymbol{e} - 2 \langle \boldsymbol{e}, \boldsymbol{n} \rangle \boldsymbol{n}

である。ただし、$\boldsymbol{e}, \boldsymbol{n}$ は正規化(=長さが1)されているものとする。 コードは次の通り。

specularReflection :: Direction3 -> Direction3 -> (Direction3, Double)
specularReflection n e
  | v == Nothing = (n, 0.0)
  | c < 0.0      = (fromJust v, -c)
  | otherwise    = (fromJust v,  c)
  where
    c = e <.> n
    v = normalize (e - (2.0 * c) *> n)

まあ、コードは上式をそのまま表現しただけである。戻り値は反射ベクトルの他に 内積 $\langle \boldsymbol{e}, \boldsymbol{n} \rangle$ すなわち $\cos \theta$ ( $\theta$ は $\boldsymbol{e}$ と $\boldsymbol{n}$ のなす角) も含めている。$\cos \theta$ は あとで使われることが多いのだ。

(屈折)

絶対屈折率 \eta_1 の物質から絶対屈折率  \eta_2 の物質に光が入射した場合、 入射角  \theta に対し、屈折ベクトルの角度が法線に対して  \phi とすると それらの間には以下の関係が成り立つという(理屈はよくわかっていない)。

 \displaystyle
\eta_1 \sin \theta = \eta_2 \sin \phi

この関係式を使うと、(途中の式変形などはよくわからないが)屈折ベクトル $\boldsymbol{v_t}$ は次の式で求められる。

 \displaystyle
\boldsymbol{v_t} = \frac{\eta_1}{\eta_2} \left( \boldsymbol{e} + \left( \cos \theta - \sqrt{{\left(\frac{\eta_2}{\eta_1} \right)}^2 + {\cos}^2 \theta - 1} \right) \boldsymbol{n} \right)

ここで $\cos \theta$ は反射のときと同じ。コードは下記のとおり。

specularRefraction :: Double -> Double -> Double -> Direction3 -> Direction3
                   -> (Direction3, Double)
specularRefraction ior0 ior1 c0 ed n
  | r <  0.0     = (o3, 0.0)
  | t == Nothing = (o3, 0.0)
  | otherwise    = (fromJust t, ior')
  where
    ior' = ior0 / ior1
    r = 1.0 / (ior' * ior') + c0 * c0 - 1.0
    a = c0 - sqrt r
    n' = if ed <.> n > 0.0 then negate n else n
    t = normalize (ior' *> (ed + a *> n'))

引数のior0ior1c0は数字が紛らわしいが、それぞれ \eta_1\eta_2 、 $\cos \theta$ を意味する。 edは入射ベクトル $\boldsymbol{e}$ である。where節を 少し説明しよう。屈折ベクトルを求める式中に平方根が含まれている。この中身が マイナスになる場合、光は全反射する(=屈折ベクトルは存在しない)。 これをチェックするため、一旦平方根の中身をrとし、ガードでマイナスの場合は 零ベクトルを返している。さて、屈折ベクトルを求めるには入射光が通ってきたのと 同じ側に向いた法線ベクトルが必要だが、本プログラムでは交点での法線は常に物体の 「外側」に向くよう定義している。つまり球なら表面から外側に法線ベクトルが向く。 しかし屈折では物体内部を光が通って物体表面にて反射と屈折が起きる。 この場合、反射屈折ベクトルの計算に必要なのは物体の内側に向いた法線ベクトルである。 そのためわざわざ $ \cos \theta $ の符号をみて $\boldsymbol{n}$ を反転させている(n')。 このことは屈折だけでなく反射も同じだが、ガードによる条件分岐で済ましている。 あとは上述の数式通りなので特に問題ないだろう。屈折はおまけとして比屈折率 ( \frac{\eta_1}{\eta_2} )も戻り値に入れている。

(Schlickの近似式)

鏡面反射・屈折では、物体表面の反射率が入射角により変化する。正確には Fresnelの公式というので求めるのだが、これがなかなかややこしい。おそらく計算 負荷も高い。そこで便利な"Schlickの近似式"というのを代わりに使うそうだ。

 \displaystyle
F_r \approx F_0 + (1 - F_0) {(1 - \cos \theta)}^5

ここで、$F_0$ は物質の法線方向からの光の反射率である。物質によりその 値が異なる。上記近似式とともに、下記のWebページにいくつかサンプルがある。

https://ja.wikipedia.org/wiki/フレネルの式 http://d.hatena.ne.jp/hanecci/20130525/p3 https://yokotakenji.me/log/math/4501/

ここで、上述の鏡面反射ベクトル計算で返される $ \cos \theta $ が必要なのだ。 コードはこちら。$ F_0 $ (Color型)と $ \cos \theta $ が引数である。 あとは上記の近似式そのままだ。

reflectionIndex :: Color -> Double -> Color
reflectionIndex (Color r g b) c =
  Color (r + (1-r) * c') (g + (1-g) * c') (b + (1-b) * c')
  where
    c' = (1.0 - c) ** 5.0

-3. フォトン追跡

鏡面反射・屈折が追加されたので、フォトン追跡もそれに合わせて いろいろ拡張せねばなるまい。これまでは拡散面だけだったので、 その反射率によって拡散反射か吸収(つまり追跡終了)しかなかったが、 フォトン追跡の条件分岐のバリエーションが増えた。

  • 拡散反射か鏡面反射か吸収か
  • 反射か屈折か

この場合分けを前述のMaterialの各パラメータとロシアンルーレットを 使って表現していくのだ。解り易いようにフローチャートで表そう。 下図に前回(拡散反射面だけ)と今回(鏡面反射・屈折追加)のフローを示そう。 前回のフローは今回のフローの一部になっている(赤点線の部分)のだ。

f:id:eijian:20180707022738p:plain

次にコードで比べてみる。前回のコードはこちら。

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

中程のrussianRouletteにて拡散反射率による分岐をしている。関数reflectは 拡散反射してさらにフォトンを追跡するものだ。最後のifは前回取り入れた 画質向上策である。

一方、今回拡張したコードはこちら。

tracePhoton :: Bool -> Material -> [Object] -> Int -> Photon
            -> IO [PhotonCache]
tracePhoton _   _   _   10 _        = return []
tracePhoton !uc !m0 !os !l !(wl, r)
  | is == Nothing = return []
  | otherwise     = do
    let
      is' = fromJust is
      (p, _, m) = is'
      d = diffuseness m
    i <- russianRoulette [d]
    ref <- if i > 0
      then reflectDiff uc m0 os l wl is'
      else reflectSpec uc m0 os l (wl, r) is'
    if (uc == False || l > 0) && d > 0.0
      then return $ ((wl, initRay p $ getDir r) : ref)
      else return ref
  where
    is = calcIntersection r os

reflectDiff :: Bool -> Material -> [Object] -> Int -> Wavelength
            -> Intersection -> IO [PhotonCache]
reflectDiff uc m0 os l wl (p, n, m) = do
  i <- russianRoulette [selectWavelength wl $ reflectance m]
  if i > 0
    then do  -- diffuse reflection
      dr <- diffuseReflection n
      tracePhoton uc m0 os (l+1) $ (wl, initRay p dr)
    else return [] -- absorption

reflectSpec :: Bool -> Material -> [Object] -> Int -> Photon -> Intersection
            -> IO [PhotonCache]
reflectSpec uc m0 os l (wl, (_, ed)) (p, n, m) = do
  let
    f0 = selectWavelength wl $ specularRefl m
    (rdir, cos0) = specularReflection n ed
    f' = f0 + (1.0 - f0) * (1.0 - cos0) ** 5.0
  j <- russianRoulette [f']
  if j > 0
    then tracePhoton uc m0 os (l+1) (wl, initRay p rdir)
    else do
      if (selectWavelength wl $ ior m) == 0.0
        then return []   -- non transparency
        else reflectTrans uc m0 os l wl ed (p, n, m) cos0

reflectTrans :: Bool -> Material -> [Object] -> Int -> Wavelength -> Direction3
             -> Intersection -> Double -> IO [PhotonCache]
reflectTrans uc m0 os l wl ed (p, n, m) c0 = do
  let
    ior0 = selectWavelength wl $ ior m0
    ior1 = selectWavelength wl $ ior m
    (tdir, ior') = specularRefraction ior0 ior1 c0 ed n
    m0' = if tdir <.> n < 0.0 then m else m_air
  tracePhoton uc m0' os (l+1) (wl, initRay p tdir)

一見して前回と比べてかなり複雑になったのがわかるだろう。russianRouletteが いくつか使われており、それが各分岐になっている。reflectDiff, reflectSpec, reflectTransはそれぞれ拡散反射、鏡面反射、屈折の場合のフォトン追跡の処理だ。

屈折をサポートしたので、最初の目的だった「集光模様」が得られる はずだ!というわけで、ガラス球を配置してフォトンマップを生成してみよう。

f:id:eijian:20180707022748p:plain

左が不透明な拡散反射面の球の場合、右がガラス球だ。左は拡散反射するので 球の形にフォトンが記録されていて下に球の影が見える。一方右はガラス球 なので球面状にフォトンは記録されず、代わりに球の影の真ん中にフォトンが 集中している部分がある。これが集光模様になるのだ!

-4. レンダリング方程式(?)

さて、いよいよレンダリングの拡張だ。フォトンマップをみる限り、想定したとおり 集光模様が描けそうな気がする!

レンダリングは、注目する点(視線レイと物体の交点)にいろいろな方向から届く 光を集計して輝度として返せば良い。ただ、言葉では簡単だがこれがなかなか難しい。 私はいまだに良い(正確な/物理的に正しい/表現力の高い/…)集計方法がわからない。 とはいえわからないなりに作るしかないので、今回採用した集計方法を説明しよう。

まずは前回までの式。

 \displaystyle
L_o = L_e + d \frac{1}{\pi} (L_l + L_d)

各文字はそれぞれ以下を表す。

  • $L_o$ : 交点から視線方向への輝度
  • $L_e$ : 物体の発光輝度
  • $L_l$ : 光源からの直接光(フォトンマップから推定することも可)
  • $L_d$ : 間接光(フォトンマップから推定)
  • $d$ : 拡散度(Material中のdiffuseness

これだとそんなに複雑ではない。該当するコードは次の通り。

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) lgts
      else radiance0
    ii = estimateRadiance pw pmap (p, n, m)

emが $L_e$ 、diが $L_l$ 、iiが $L_d$に相当する。直接光と間接光は 別々にBRDFを適用しているのでちょっとわかりにくいが、基本的に上記の式を そのまま適用している。

鏡面反射と屈折が追加されたことで、輝度計算式は次のように拡張した。

 \displaystyle
L_o = \frac{1}{2 \pi} L_e + d \frac{1}{\pi} (L_l + L_d) + (1 - d) ( f L_s + (1 - m) ( (1 - f) L_t ) )

各文字はそれぞれ以下を表す。

  • $L_o$ : 交点から視線方向への輝度
  • $L_e$ : 物体の発光輝度
  • $L_l$ : 光源からの直接光(フォトンマップから推定することも可)
  • $L_d$ : 間接光(フォトンマップから推定)
  • $L_s$ : 鏡面反射方向の輝度
  • $L_t$ : 屈折方向の輝度
  • $d$ : 拡散度(Material中のdiffuseness
  • $f$ : 交点の反射率 (Schlickの近似式で計算)
  • $ m $ : 金属性(Material中のmetalness

なお、$ d, f, m $ はいずれも $ 0 \sim 1 $の間の値をとる。コードは次の通り。 (そういえば物体の発光輝度に $\frac{1}{2 \pi}$ を掛けているのはなぜだっけ?)

traceRay :: Screen -> Material -> Int -> PhotonMap -> [Object] -> [Light] -> Ray -> IO Radiance
traceRay _    _   10 _     _     _     _  = return radiance0
traceRay !scr !m0 !l !pmap !objs !lgts !r
  | is == Nothing = return radiance0
  | otherwise     = do
    si <- if df == 1.0 || f == black
      then return radiance0
      else traceRay scr m0 (l+1) pmap objs lgts (initRay p rdir)
    ti <- if f' == black || ior1 == 0.0
      then return radiance0
      else do
        let
          ior0 = averageIor m0
          (tdir, ior') = specularRefraction ior0 ior1 cos0 (getDir r) n
          m0' = if tdir <.> n < 0.0 then m else m_air
        traceRay scr m0' (l+1) pmap objs lgts (initRay p tdir)
    return (sr_half    *> emittance m +
            df         *> brdf m (di + ii) +
            (1.0 - df) *> (f <**> si + (1.0 - mt) *> f' <**> ti))
  where
    is = calcIntersection r objs
    (p, n, m) = fromJust is
    di = if useClassicForDirect scr
      then foldl (+) radiance0 $ map (getRadianceFromLight objs p n) lgts
      else radiance0
    ii = estimateRadiance scr pmap (p, n, m)
    (rdir, cos0) = specularReflection n (getDir r)
    df = diffuseness m
    mt = metalness m
    f = reflectionIndex (specularRefl m) cos0
    f' = negateColor f                         -- this means '1 - f'
    ior1 = averageIor m

反射屈折方向にトレースするかどうかの判定など多少ややこしい部分はあるが、 基本的には上記式をそのまま再現している。

さて、それでは拡張した今回のプログラムで数種類の物質の球をレンダリング してみよう。

f:id:eijian:20180707022721p:plain

鏡面反射が加わったことで、かなり表現力が高まった気がする。また、ガラス球が 表現できるようになり念願の集光模様が再現できた!

-5. おまけ:設定ファイルに対応

実はこれまではレンダリングする画像の情報(スクリーン情報と物体の情報)は ソースコードに直書きしていた。流石にシーンを変更するたびにソースを修正して リコンパイルするのは面倒臭く、またさまざまなシーンを描画するために画像情報を 残しておきたいことから、設定ファイルから情報を読み込むように改めた。

まずスクリーン情報(カメラの位置や向き、解像度、アンチエリアス要否など)を 表すファイル(Screenファイル)の例。フォーマットは(エセ)YAMLとした。

nphoton       : 100000
xresolution   : 256
yresolution   : 256
antialias     : yes                # yes or no
samplephoton  : 100
useclassic    : yes                # yes or no
estimateradius: 0.3
ambient       : [ 0.001, 0.001, 0.001 ]
maxradiance   : 0.01
eyeposition   : [ 0.0, 2.0, -4.5 ]
targetposition: [ 0.0, 2.0, 0.0 ]
upperdirection: [ 0.0, 1.0, 0.0 ]
focus         : 2.7
photonfilter  : none               # cone, gauss

こちらは単純な変数名と値の組だけなのでなんのことはない。一方シーン情報 (物体の形や配置、色などを定義)の例としてガラス球の場合を示そう。 ちょっと長いが勘弁願いたい。

light:
  - type     : parallelogram
    color    : [ 1.0, 1.0, 1.0 ]
    flux     : 5.0
    position : [ -0.5, 3.99, 2.5 ]
    dir1     : [ 1.0, 0.0, 0.0 ]
    dir2     : [ 0.0, 0.0, 1.0 ]

material:
  - type         : solid
    name         : mwall
    emittance    : [ 0.0, 0.0, 0.0 ]
    reflectance  : [ 0.5, 0.5, 0.5 ]
    transmittance: [ 0.0, 0.0, 0.0 ]
    specularrefl : [ 0.8, 0.8, 0.8 ]
    ior          : [ 0.0, 0.0, 0.0 ]
    diffuseness  : 1.0
    metalness    : 0.0
    smoothness   : 0.0
  - type         : solid
    name         : mwallr
    emittance:     [ 0.0, 0.0, 0.0 ]
    reflectance:   [ 0.4, 0.1, 0.1 ]
    transmittance: [ 0.0, 0.0, 0.0 ]
    specularrefl:  [ 0.0, 0.0, 0.0 ]
    ior:           [ 0.0, 0.0, 0.0 ]
    diffuseness:   1.0
    metalness:     0.0
    smoothness:    0.0
  - type         : solid
    name: mwallb
    emittance:     [ 0.0, 0.0, 0.0 ]
    reflectance:   [ 0.1, 0.1, 0.4 ]
    transmittance: [ 0.0, 0.0, 0.0 ]
    specularrefl:  [ 0.0, 0.0, 0.0 ]
    ior:           [ 0.0, 0.0, 0.0 ]
    diffuseness:   1.0
    metalness:     0.0
    smoothness:    0.0
  - type         : solid
    name: mparal
    emittance:     [ 0.7958, 0.7958, 0.7958 ]
    reflectance:   [ 0.0, 0.0, 0.0 ]
    transmittance: [ 0.0, 0.0, 0.0 ]
    specularrefl:  [ 0.0, 0.0, 0.0 ]
    ior:           [ 0.0, 0.0, 0.0 ]
    diffuseness:   0.0
    metalness:     0.0
    smoothness:    0.0
  - type         : solid
    name         : glass
    emittance:     [ 0.0, 0.0, 0.0 ]
    reflectance:   [ 0.0, 0.0, 0.0 ]
    transmittance: [ 0.0, 0.0, 0.0 ]
    specularrefl:  [ 0.08, 0.08, 0.08 ]
    ior:           [ 1.5, 1.5, 1.5 ]
    diffuseness:   0.0
    metalness:     0.0
    smoothness:    0.0

vertex:
  - cl01 : [ -0.5, 3.99, 2.5 ]
  - cl02 : [ 0.5, 3.99, 2.5 ]
  - cl03 : [ -0.5, 3.99, 3.5 ]

object:
  - type    : plain
    name    : flooring
    normal  : [ 0.0, 1.0, 0.0 ]
    position: [ 0.0, 0.0, 0.0 ]
    material: mwall
  - type    : plain
    name    : ceiling
    normal  : [ 0.0, -1.0, 0.0 ]
    position: [ 0.0, 4.0, 0.0 ]
    material: mwall
  - type    : plain
    name    : rsidewall
    normal  : [ -1.0, 0.0, 0.0 ]
    position: [ 2.0, 0.0, 0.0 ]
    material: mwallb
  - type    : plain
    name    : lsidewall
    normal  : [ 1.0, 0.0, 0.0 ]
    position: [ -2.0, 0.0, 0.0 ]
    material: mwallr
  - type    : plain
    name    : backwall
    normal  : [ 0.0, 0.0, 1.0 ]
    position: [ 0.0, 0.0, -6.0 ]
    material: mwall
  - type    : plain
    name    : frontwall
    normal  : [ 0.0, 0.0, -1.0 ]
    position: [ 0.0, 0.0, 5.0 ]
    material: mwall
  - type    : sphere
    name    : ball_glass
    center  : [ 0.0, 0.8, 3.0 ]
    radius  : 0.8
    material: glass
  - type    : parallelogram
    name    : ceiling_light
    pos1    : cl01
    pos2    : cl02
    pos3    : cl03
    material: mparal

以前、CPU回で使ったParsecライブラリをここでも使ってみた。 CPU回でなんとなくパーサの感触をつかんでいたのであまりハマることなく 実装できた。Parsecを使ったパーサ作りは、乗ってくるとどんどん パーサを高度化していけるので楽ちんだ。先人に感謝。

実行方法は次の通り(ガラス球の例、ex-11.6)。

$ cabal build
$ dist/build/pm/pm example/screen0.scr example/ex-11.6.scene | dist/build/rt/rt example/screen0.scr example/ex-11.6.scene | convert - ~/tmp/ex-11.6.png

これで ~/tmp/ex-11.6.pngという画像ファイルが出来上がる。 ただ、同じ情報ファイルをフォトンマップ作成(pm)とレイトレーシング(rt)で 繰り返すのは面倒だし画像フォーマット変換(ImageMagickconvertを使用)も 毎回書くのは邪魔臭いのでシェルスクリプトにすることにした。

$ util/drv.sh example/screen0.scr example/ex-11.6.scene ~/tmp/ex-11.6.png

2. 作例

上記の他にも少し作例を示そう。

f:id:eijian:20180707022728p:plain

左上は天井からの日光を鏡面三角柱で壁に反射している図、右上は例の フォトンマッピング本でなんども出てくる図。これがやりたかった! 下の二つは球以外の形状をガラスで作ってみたもの。最後に集光模様の 面白い例として光の波長により屈折率が異なるガラス球の画像。屈折率の違いに よりプリズムのように虹色(?)に分かれた集光模様になる。ただし、本プログラムは RGB三種類のフォトンしか使っていないため、青より短波長の紫が表現できない。 そこまでやるにはRGBに紫を追加してフォトンの種類を増やして追跡しないとだめ だろう。ちょっと面倒なので宿題にしておこう。

f:id:eijian:20180707023003p:plain

3. まとめ

今回でやっと集光模様までたどり着いた。最初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画像認識プログラムの「順伝播」のみ実装完了した。 出力結果にはほとんど意味はないが、各層の処理を経てなんとなく結果が出力 されたので雰囲気はつかめたかなと。

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

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