読者です 読者をやめる 読者になる 読者になる

Haskellでいってみよう

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

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の逆伝播完成?

deep learning

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

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): そして逆伝播(でも全結合層まで)

deep learning

第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): まずは順伝播(下)

deep learning

前回に続き、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): まずは順伝播(上)

deep learning

さて、今回からは新たなネタとして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で使われるいろいろな型を定義した。
  • 入力値を変換するいろいろな「層」を定義した。
  • これらを使って順伝播の処理を実装した。

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

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

CPUの創りかた(10): おまけ、アセンブラ

logic circuit

CPU自体は前回までで完成したので、次のネタに行ってもよかったのだが、

  • CPU(の実行プログラム)に与えるのがマシンコードだと、いちいち ハンドアセンブルするのが面倒くさい。
  • 別件でパーサを書く必要がありParsecライブラリに興味を持っていたので そのうち勉強したいと思っていた。

という理由により、今回はおまけとして簡易アセンブラを作ってみよう。

まずは文法の定義から

このCPU(TD4)は命令数が非常に少ないので、何となれば、単純な文字列 変換でもよかったのだが、Parsecの勉強も兼ねているからそこはちゃんと 文法の定義から入らなければなるまい。

文法定義とくればBNFだろう。大学時代に、わからないなりに適当な文法を 定義して遊んでいたのを思い出す。その頃の記憶を頼りにBNFで書き始めたのだが、 どうやら世間には EBNF というのがあるようなのでどうせならこちらにしよう、BNFの問題を 解決しているということだし。

ところでアセンブリ言語(ニーモニック?)の文法はどのように定義したら いいのだろう? と調べたら、 こういうの が見つかった。作りたいのはこれよりだいぶ簡素なものなので、ありがたく 参考にさせていただいた。EBNFはこれが初めてなので書き方が間違っているかも しれないが、とりあえず次のように定義してみた。

# EBNF for TD4 assembler

  program = instcode , { instcode } ;

  instcode = ( code2 | code1 ) , linefeed ;

  code2 = inst2, space, operand2 ;
  code1 = inst1, space, operand1 ;

  inst2 = 'add' | 'mov' ;
  inst1 = 'in' | 'out' | 'jnc' | 'jmp' ;

  operand2 = register , "," , operand1 ;
  operand1 = register | imdata ;

  register = 'a' | 'b' ;

  imdata = 4 * digit2 ;
  digit2 = "0" | "1" ;

  space = white space , { white space } ;
  white space = '\x20' | '\x09' ;

  linefeed = '\x0a' | ( '\x0d' , '\x0a' ) ;

2オペランド命令のANDとMOV、それ以外のは1オペランド命令なので定義が異なる。 あとはレジスタとか直値とかを定義していけばよさそう。ひとまず形になった 気がするので先に進もう。

パーサの前に

どうもParsecを使えば、EBNF(BNF)の定義からある程度容易にパーサの プログラムに落とし込める、という話があるらしい。ならEBNFができている のだから簡単に作れるよね、と思ったのは甘すぎだった。

まず参考となるWebサイトをいろいろ読んでみて、こう作ればよさそうという 感触を得たかったのだが、ちっとも理解できなかった。今なら何を理解できて いなかったのか分かるが、当初は本当にわからなかったのだ。

要は、ソースをパーサが字句解析したあと、それをどのようにして次の処理 (意味解析、コード生成)に結びつけたらいいのかイメージが湧かなかったため。 Web上でありがちな例は、"add 1,2"みたいなのが入力されたら、"add"を 解析して「"add"という文字列を返す」みたいなやつだ。文字列を解析したい のに解析結果が文字列ならそれをさらに解釈する処理が必要で堂々巡りに 思えてしまった。転機はこの記事。 この記事では明快に「返す値を型で定義している」。

つまりパーサからは「解析木」かそれに類する構造化されたデータが出力されるのだ。 それをまず考えて(定義して)おかないとパーサが書けるわけがない。 ではどのような構造のデータがあればマシンコードに変換できるかを考えてみる。

TD4の命令は、命令コードと1つまたは2つのオペランドからなる(オペランドなしの 命令は無い)。そこで、これらを組で表すことにする。最終的にはこう。

data Inst = Add | Mov | In | Out | Jnc | Jmp deriving (Enum, Show)
data Operand = RegA | RegB | Imdata String deriving Show

type Mnemonic = (Inst, (Operand, Maybe Operand))

Mnemonicは二重の組だ。命令コードとオペランドの組からなる。 オペランドの組は、一つ目は必ず存在するので普通にOperandで、 二つ目は無いかもしれないのでMaybe型とした。オペランドはA,Bレジスタか 直値の三種類しかない。パーサはソースを解析してMnemonicを命令の数だけ リストにして返してくれればよいのだ(以後、解析リスト、とする)。

パーサを書く

ではパーサを書いてみようと思う。何となくだが末端から作って積み重ねて 行った方がわかりやすそうなので、個々のオペランドから始める。まずは 直値(immediate data)から。

-- EBNF
--   imdata = 4 * digit2 ;
--   digit2 = "0" | "1" ;

imdata :: Parser Operand
imdata = do
  im <- count 4 (oneOf "01")
  return $ Imdata im

先に定義したように、直値もOperand型の一つなのでパーサの型が Parser Operandになっている。直値は二進数だけを扱うことにし、 かならず4桁と決めた。なので、oneOfで0か1に限定し、それを4つ 連続して取り出したら返すようにした。EBNFの定義と見比べると、 決めた通りにプログラムを書けばよいのがわかる。

同様に、A,Bレジスタはこうなる。

-- EBNF
--   register = 'a' | 'b' ;

register :: Parser Operand
register = do
  rg <- (regA <|> regB)
  return rg

regA :: Parser Operand
regA = do
  rg <- string "a"
  return $ RegA

regB :: Parser Operand
regB = do
  rg <- string "b"
  return $ RegB

AかBかの区別をつけるため、それぞれ別にパーサを定義し、それを<|>で 合わせてレジスタのパーサとした。

あとは同様に、EBNFをもとに以下のようなパーサを作った。

program :: Parser [Mnemonic]
program = do
  pg <- many1 $ instcode
  return pg

instcode :: Parser Mnemonic
instcode = do
  cd <- code2 <|> code1
  many1 $ oneOf "\r\n"
  return cd

code2 :: Parser Mnemonic
code2 = do
  in2 <- inst2
  many1 space
  op2 <- operand2
  return (in2, op2)

code1 :: Parser Mnemonic
code1 = do
  in1 <- inst1
  many1 space
  op1 <- operand1
  return (in1, (op1, Nothing))

inst2 :: Parser Inst
inst2 = do
  i2 <- (string "add" <|> string "mov")
  let i = if i2 == "add" then Add else Mov
  return i

inst1 :: Parser Inst
inst1 = do
  i1 <- (string "in" <|> string "out" <|>
         try (string "jnc") <|> (string "jmp"))
  let i = case i1 of
            "in"  -> In
            "out" -> Out
            "jnc" -> Jnc
            "jmp" -> Jmp
  return i

operand2 :: Parser (Operand, Maybe Operand)
operand2 = do
  op2 <- register
  char ','
  op1 <- operand1
  return $ (op2, Just op1)

operand1 :: Parser Operand
operand1 = do
  op1 <- (register <|> imdata)
  return op1

EBNFの定義と見比べれば、それぞれ何をしているかは分かると思う。 ではテストしてみよう。テスト用プログラムはざっと以下のような感じ。

module Main where

import Control.Applicative hiding ((<|>))
import Data.Char
import Text.Parsec
import Text.Parsec.String
import Text.Parsec.Char

(型、パーサ定義は省略)

main :: IO ()
main = do
  parseTest inst2 "add"
  parseTest inst2 "mov"
  parseTest inst2 "abc"
  parseTest register "a"
  parseTest register "b"
  parseTest register "c"
  parseTest imdata "0100"
  parseTest imdata "10100"
  parseTest imdata "1012"
  parseTest operand1 "aa"
  parseTest operand2 "a 1100"
  parseTest code2 "add a,b"
  parseTest code2 "mov a,0011"
  parseTest code1 "jmp 1011"
  parseTest code1 "in  b"
  parseTest code1 "in  a"
  parseTest code1 "in  0110"
  parseTest instcode "add a,b\n"
  parseTest instcode "add a,b\r\n"
  parseTest instcode "add a,a"
  parseTest instcode "jmp 1100\r\n"

このように、パーサに食わせてみたい文字列とそれを処理するパーサ関数を 組みにして与えればいいようだ。パーサの解析に失敗した(=文法にそぐわない) 場合はエラーが返るのですぐわかる。ざっとみたところうまく動いている ように・・・見えなかった!

文法定義のやり直し

うまくいったと思っていたのに実はいけていなかった。add a,aとか in aという命令はTD4には存在しないがパーサを通ってしまうのだ。 文法定義が甘かったわけだ。2オペランド命令と1オペランド命令を 分けるだけでは個々の命令の細かい制限(Bレジスタしか指定できない等)が 表現できていなかった。

ということで、EBNFの定義を見直す。

EBNF for TD4 assembler (revision 1)

  program  = line , { line } ;
  line      = instcode , linefeed ;

  instcode = inst_add | inst_mov | inst_in | inst_out | inst_jump ;
  inst_add  = 'add' , space , register , "," , imdata ;
  inst_mov  = 'mov' , space , ( op_mov1 | op_mov2 | op_mov3 ) ;
  inst_in   = 'in'  , space , register ;
  inst_out  = 'out' , space , ( reg_b | imdata) ;
  inst_jump = ( 'jnc' | 'jmp' ) , space , imdata ;

  op_mov1   = register , "," , imdata ;
  op_mov2   = reg_a    , "," , reg_b ;
  op_mov3   = reg_b    , "," , reg_a ;

  register  = reg_a | reg_b ;
  reg_a     = 'a' ;
  reg_b     = 'b' ;

  (imdata以降は同じなので割愛)

ポイントは、個々の命令ごとにオペランドのパターンを記述したこと。 結局オペランドが共通なのはジャンプ命令だけだった。これに沿って パーサも書き直す。

instcode :: Parser Mnemonic
instcode = do
  cd <- inst_add <|> inst_mov <|> inst_in <|> inst_out <|> inst_jump
  return cd

inst_add :: Parser Mnemonic
inst_add = do
  in1 <- string "add"
  many1 space
  rg1 <- register
  char ','
  im1 <- imdata
  return (toInst in1, (rg1, Just im1))

inst_mov :: Parser Mnemonic
inst_mov = do
  in1 <- string "mov"
  many1 space
  op  <- try (op_mov1) <|> (try op_mov2 <|> op_mov3)
  return (toInst in1, op)

inst_in :: Parser Mnemonic
inst_in = do
  in1 <- string "in"
  many1 space
  rg <- register
  return (toInst in1, (rg, Nothing))

inst_out :: Parser Mnemonic
inst_out = do
  in1 <- string "out"
  many1 space
  op <- regB <|> imdata
  return (toInst in1, (op, Nothing))

inst_jump :: Parser Mnemonic
inst_jump = do
  in1 <- try (string "jnc") <|> (string "jmp")
  many1 space
  im <- imdata
  return (toInst in1, (im, Nothing))

op_mov1 :: Parser (Operand, Maybe Operand)
op_mov1 = do
  rg <- register
  char ','
  im <- imdata
  return (rg, Just im)

op_mov2 :: Parser (Operand, Maybe Operand)
op_mov2 = do
  op1 <- regA
  char ','
  op2 <- regB
  return (op1, Just op2)

op_mov3 :: Parser (Operand, Maybe Operand)
op_mov3 = do
  op1 <- regB
  char ','
  op2 <- regA
  return (op1, Just op2)

toInst :: String -> Inst
toInst s = case s of
             "add" -> Add
             "mov" -> Mov
             "in"  -> In
             "out" -> Out
             "jnc" -> Jnc
             "jmp" -> Jmp

先ほどと同じテストを流してみると・・・ちゃんとmov a,ain aが エラーになっている! パーサが出来上がった! (Applicativeスタイルへの対応はまた今度。上記の書き方だけでも まだ全然理屈を咀嚼できていない)

コード生成

さて、パーサはソースプログラムの解析をするだけなので、そこから 目的のマシンコードを作り出す処理が必要だ。次はコード生成部分を作ろう。 と言っても、パーサがきちんと解析リストを生成してくれれば あとは簡単だ。AならBというふうに対応するマシンコードを返せば良い。 TD4は命令数が非常に少ないので、全部列挙することにした。

generate :: Either ParseError [Mnemonic] -> [String]
generate (Left s)       = [show s]
generate (Right [])     = []
generate (Right (x:xs)) = (translateOne x):(generate $ Right xs)

translateOne :: Mnemonic -> String
translateOne (Add, (RegA, Just (Imdata s))) = "0000" ++ s
translateOne (Mov, (RegA, Just RegB))       = "00010000"
translateOne (In , (RegA, Nothing))         = "00100000"
translateOne (Mov, (RegA, Just (Imdata s))) = "0011" ++ s
translateOne (Mov, (RegB, Just RegA))       = "01000000"
translateOne (Add, (RegB, Just (Imdata s))) = "0101" ++ s
translateOne (In , (RegB, Nothing))         = "01100000"
translateOne (Mov, (RegB, Just (Imdata s))) = "0111" ++ s
translateOne (Out, (RegB, Nothing))         = "10010000"
translateOne (Out, (Imdata s, Nothing))     = "1011" ++ s
translateOne (Jnc, (Imdata s, Nothing))     = "1110" ++ s
translateOne (Jmp, (Imdata s, Nothing))     = "1111" ++ s
translateOne _                              = error "no such mnemonic"

最後の行は保険だ。パーサがきちんと解析できていれば定義されていない 命令の解析リストが入力されることは無いだろうから。コンパイルして実行 してみる。

$ cabal build td4asm
  :
$ echo "mov a,b" | dist/build/td4asm/td4asm
00010000

おお!正しく変換されている!ちなみに、入力の最後に改行が無いとエラーになる。

$ echo -n "mov a,b" | dist/build/td4asm/td4asm
"TD4 asm" (line 1, column 8):
unexpected end of input
expecting new-line

構文解析で、命令行の最後は改行で終わるように定義しているからだ・・・。 これぐらいの制約は多めに見てもらおう。

再びラーメンタイマーの実行

前回はラーメンタイマープログラムのマシンコードを手で作ってCPUに入れていた。 今回はアセンブラアセンブリ言語のソースを入れてマシンコードを生成させて CPUへ入れたいと思う。アセンブラは生成したマシンコードを標準出力に出すので そのままパイプでCPUへつないであげればいいのだ。

* 見やすくするため実際には指定しているディレクトリなどは省いている *

$ cat ramen.a
out 0111
add a,0001
jnc 0001
add a,0001
jnc 0011
out 0110
add a,0001
jnc 0110
add a,0001
jnc 1000
out 0000
out 0100
add a,0001
jnc 1010
out 1000
jmp 1111

$ td4asm < ramen.a | td4
clock 1.0 sec; I/P 0000
step 0; [CF:0][A:0000][B:0000][OP:0000][PC:0000]
step 1; [CF:0][A:0000][B:0000][OP:0111][PC:0001]
step 2; [CF:0][A:0001][B:0000][OP:0111][PC:0010]
^C

ちゃんと動いている! やはりハンドアセンブルより断然楽ちんだ。

まとめ

さて、今回で本当に"CPU回"はおしまい。記事を細切れにしたせいもあって 10回にもなってしまった。長丁場だったが、CPUネタは面白い取り組みだったなあ。 どこかにi4004の回路図落ちてないかな :-)

(ここまでのソース)

CPUの創りかた(9): CPUはじめました

logic circuit

さあ、前回までで必要なモジュールは出揃った。今回はそれらを組み立てて 動くCPUを作ってしまおう! 一気に最終形はしんどいので少しずつピースを埋めていく感じで進めていきたい。

なお、しつこいようだがここで作っているCPUは以下の本で解説されている TD4という名前のオリジナル4bit CPUだ。説明中にtd4と出てくるのは その名前である。

ステップ 0: 電源、クロックジェネレータ(に相当するところ)

これまで論理回路の細かいところやCPU内の各種モジュールを作ることばかり やってきて、実行できるプログラムにする部分には目を瞑っていた。 しかしさすがに今回は「プログラム」を動かしたいのでそうはいかない。

そこでステップ0として完動させるための周辺部分を作っていこう。 電子工作では電源モジュールだとかクロックジェネレータとかその他の アナログ回路部分に相当するだろうか。

まず仕様を列挙しよう。

  • 「プログラム」は標準入力から投入する。
  • 「プログラム」は'0'と'1'の連続した文字列とする。また間にホワイトスペースが いくつ入ってもよい。
  • 「プログラム」におけるビット並びは(慣れているので) MSB...LSB の順とする。
  • プログラムカウンタが4 bit なので「プログラム」は16 bytes = 128文字。 ただしそれより少ない場合は'0'で補填する。多い場合は切り捨てる。
  • コマンドライン引数は順に"クロック間隔"と"入力ポート"の2つ。省略可能だが、 クロック間隔だけを省略することはできない。
  • クロック間隔の単位は秒、小数も使える。デフォルト値は1.0秒。入力ポートは 4桁の二進数でデフォルト値は"0000"。
    • 例) td4 0.5 0101 < program

ではこの仕様に基づいて作っていこう。

main :: IO ()
main = do
  pg <- getContents
  opts <- getArgs
  let (clock, iport) = parseOpts opts
  putStrLn ("clock " ++ (show clock) ++ " sec; I/P " ++ toStr iport)
  -- CLR(1),CF(1),A(4),B(4),OP(4),PC(4)
  let stat = toBits "011000010011000010"
  loop 0 clock lc_td4 stat iport (createRom pg)

getContentsで標準入力を読み、getArgsコマンドライン引数を取り込む。 どちらも標準で用意されている関数だ。parseOptsでオプションを解析している。

parseOptsは次の通り。 コマンドライン引数の数に応じてその値を読み込んだりデフォルト値を 使ったりしている。

defClock :: Double
defClock = 1.0              -- default clock time = 1 sec
defInput :: [Bin]
defInput  = toBits "0000"   -- default value of Input port

parseOpts :: [String] -> (Double, [Bin])
parseOpts [] = (defClock, defInput)
parseOpts (x:[]) = ((read :: String -> Double) x, defInput)
parseOpts (x:y:_) = ((read :: String -> Double) x, toBits y)

次は「プログラム」の整形についてだ。上記仕様ではビットの並びは MSB...LSBだが、これまで作ってきた論理回路モジュールでは、 入力(Binの配列)が全てLSB...MSBの順だ。そこで前もって順序を入れ替えておこう。 それを行っているのがcreateRom

createRom :: String -> [Bin]
createRom rs = concat $ map reverse $ split8 rs'
  where
    rs' = take 128 (toBits rs ++ repeat sLO)   -- 128 bits = 16 bytes

そのまんまだが。。。入力された文字列をtoBitsBinの配列にし、 足りなければsLO(=0)を付け加えて、先頭から16 bytes取り出している。 このような大雑把な記述が可能なのはHaskellの遅延評価のおかげだなぁ。 ちなみにtoBitsは0と1以外の文字は無視するので、間にスペースや 改行があっても問題ない。あとは8 bits単位に切り出してそれぞれを 逆順に並べ替えれば完成だ。

さあ、いよいよCPUモジュールを駆動する(呼び出す)ところだ。これは 「クロックの立ち上がり」のたびに関数を呼び出す無限ループである。 前回までに解説したように状態はCPUの外で管理することにしたので、 入力はCPUの状態+ROMの内容、出力はCPUの最新状態だ。 それをループにしたいのだ。出力値を次の入力値(の一部)に使うので、 再帰呼び出しが良さそうだ。ということで次のようなloop関数を作ってみた。

loop :: Int -> Double -> LogicCircuit -> [Bin] -> [Bin] -> [Bin] -> IO ()
loop s w lc st ip pg = do
  let os = lc (st ++ ip ++ pg)
  putStatus s os
  threadDelay $ floor (w * 1000 * 1000)
  -- set CLR to HI and take status from output
  let st' = [sHI] ++ (take 17 os)
  loop (s+1) w lc st' ip pg

クロックの度に状態を画面に出力したいので、ループの数(=ステップ数s)を 引数の最初に入れている。次はクロック間隔w、3番目(lc)がCPUを表す関数 (以後、CPU関数と呼ぼう)だ。 CPU関数への入力は「状態」「入力ポート値」「ROM」の3つ。

関数を呼び出して得た出力を画面に表示(putStatus)し、クロック間隔だけ待ち (threadDelay)、状態を更新して次のステップを呼び出す。この繰り返し。

次回の入力値を作っている少々奇妙な部分について。

  let st' = [sHI] ++ (take 17 os)

入力の最初の値が「リセット信号」を表しており、これがLOだとリセットが かかるようになっている。だから一番最初の呼び出し以外はHIにしないといけない。 あと、出力値から先頭の17個を取り出しているがこれには以下が含まれている。

  • carryフラグ(1 bit)
  • Aレジスタ(4 bit)
  • Bレジスタ(4 bit)
  • 出力ポート値(4 bit)
  • プログラムカウンタ(4 bit)

これで"評価ボード"(?)ができた。早速ダミーのCPU関数で動かしてみよう。 中身は何もせず入力を出力に回すだけ。

lc_td4_st0 :: LogicCircuit
lc_td4_st0 xs = concat [cf, a, b, op, pc]
  where
    [_, cf, a, b, op, pc, _, _] = splitInput xs

splitInput :: [Bin] -> [[Bin]]
splitInput xs = [cl, cf, a, b, op, pc, ip, rom]
  where
    (cl, xs0) = splitAt 1 xs
    (cf, xs1) = splitAt 1 xs0
    (a , xs2) = splitAt 4 xs1
    (b , xs3) = splitAt 4 xs2
    (op, xs4) = splitAt 4 xs3
    (pc, xs5) = splitAt 4 xs4
    (ip, rom) = splitAt 4 xs5

splitInputで分割して必要なものを取り出して並べているだけ。 さあコンパイルして実行してみよう。

$ cabal configure
Resolving dependencies...
Configuring mkcpu-0.1.0.0...

$ cabal build
Building mkcpu-0.1.0.0...
Preprocessing executable 'td4' for mkcpu-0.1.0.0...
[7 of 7] Compiling Main             ( src/Main-td4.hs, dist/build/td4/td4-tmp/Main.o )
Linking dist/build/td4/td4 ...

$ echo "0000" | dist/build/td4/td4 
clock 1.0 sec; I/P 0000
step 0; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
step 1; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
step 2; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
step 3; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
step 4; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
step 5; [CF:1][A:0001][B:0010][OP:0011][PC:0100]
^C

レジスタなどの状態が表示されている。ちなみにその適当な値は、 実はmainの中で指定してある。

  -- CLR(1),CF(1),A(4),B(4),OP(4),PC(4)
  let stat = toBits "011000010011000010"
  loop 0 clock lc_td4_st0 stat iport (createRom pg)

このstatだ。A、B、OP、PCの値はそれぞれ1、2、3、4にセットしてあるのだ。 先述の通り一番最初のビットはリセット(CLR)であり、最初だけは'0'に してある。が、この何もしないダミーCPUではリセット信号が使われないので Aレジスタなどは初期値が入ったまま(のように見えるの)だ。

兎にも角にも、まずはCPUを駆動する周辺回路に相当する部分は一応 動いたようだ。これを使ってCPUを最終形まで組み立てていこう。

ステップ 1: レジスタの使用

CPUは本来状態を保持したり更新したりして処理を進めていくものだ。 状態はレジスタに保持されているのだが、以前の回で書いたように、CPUの 1サイクルの最終段階でレジスタを更新(もしくは保持)している。 この部分だけを作ってみよう。前回示したブロック図では一番右端にある部分だ。

f:id:eijian:20160504150006p:plain

コードは以下。

lc_td4_st1 :: LogicCircuit
lc_td4_st1 xs = concat [cf', a', b', op', pc']
  where
    [cl, cf, a, b, op, pc, _, _] = splitInput xs
    v0  = toBits "0000"
    cf' = take 1 $ lc_dff_cp (cl ++ [sHI] ++ cf)
    a'  = lc_register4 (cl ++ [sHI] ++ a  ++ v0)
    b'  = lc_register4 (cl ++ [sHI] ++ b  ++ v0)
    op' = lc_register4 (cl ++ [sHI] ++ op ++ v0)
    pc' = lc_counter4  (cl ++ [sHI] ++ pc ++ v0)

入力を切り出す部分は同じ。v0はダミー値だ。 フラグやレジスタの入力値をそれぞれレジスタモジュールやカウンタモジュールへ 入れているだけだ。またリセット信号(cl)もそれぞれに入れている。 実行してみよう。

$ echo "0000" | dist/build/td4/td4 
clock 1.0 sec; I/P 0000
step 0; [CF:0][A:0000][B:0000][OP:0000][PC:0000]
step 1; [CF:0][A:0000][B:0000][OP:0000][PC:0001]
step 2; [CF:0][A:0000][B:0000][OP:0000][PC:0010]
step 3; [CF:0][A:0000][B:0000][OP:0000][PC:0011]
step 4; [CF:0][A:0000][B:0000][OP:0000][PC:0100]
step 5; [CF:0][A:0000][B:0000][OP:0000][PC:0101]
^C

ステップ0の結果とはだいぶ変わっている。まず、リセット信号が入ったため、 Aレジスタなどの初期入力値は一旦クリアされて0になっているのがわかる。 さらに、ステップが進む毎にプログラムカウンタ(PC)がカウントアップされて いる!カウンタモジュールは前に作ってテストしているから当然こうなるの だが、実行プログラムとしてこの出力になるのはちょっと嬉しい! (CPUが動いているぞ、という感じがする)

ステップ 2: 加算器の追加

次に加算器を取り付けよう。加算器には入力が2つ必要だが、状態を確認する ためにAレジスタの値を使う。もう一方の値はROMから無理やり取り出そう。 ROMにはプログラムカウンタをつないで0番地から順に値を取り出すようにする。 取り出した8bitから下4bitを使ってAレジスタに足し、結果がAレジスタ入るように 配線する。もちろんcarryフラグも更新する。プログラムはこうだ。

lc_td4_st2 :: LogicCircuit
lc_td4_st2 xs = concat [cf', a', b', op', pc']
  where
    [cl, _, a, b, op, pc, _, rom] = splitInput xs
    rdata = lc_rom16 (pc ++ rom) -- get data addressed by PC
    v0  = toBits "0000"
    im  = take 4 rdata
    (s0, c0) = splitAt 4 $ lc_adder (a ++ im)
    cf' = take 1 $ lc_dff_cp (cl ++ [sHI] ++ c0)
    a'  = lc_register4 (cl ++ [sLO] ++ a  ++ s0)
    b'  = lc_register4 (cl ++ [sHI] ++ b  ++ v0)
    op' = lc_register4 (cl ++ [sHI] ++ op ++ v0)
    pc' = lc_counter4  (cl ++ [sHI] ++ pc ++ v0)

4行目でROMの現在番地の値から下4bitを取り出し、5行目(lc_adderのある行)で Aレジスタと足しあわせている。それをs0, c0にしてそれぞれAレジスタと carryフラグへ入れている。 Aレジスタの方は引数の2つ目(レジスタのLD入力)をHIではなくLOにしている。 これは保持している値ではなく外から入った値(s0)をセットするためだ。

プログラム(とは言えないが)は下4桁に加算したい数字を記載している。 上から、1,2,3,4,5,1である。

$ cat program
00000001
00000010
00000011
00000100
00000101
00000001

実行してみよう。

$ dist/build/td4/td4 < program
clock 1.0 sec; I/P 0000
step 0; [CF:0][A:0000][B:0000][OP:0000][PC:0000]
step 1; [CF:0][A:0001][B:0000][OP:0000][PC:0001]
step 2; [CF:0][A:0011][B:0000][OP:0000][PC:0010]
step 3; [CF:0][A:0110][B:0000][OP:0000][PC:0011]
step 4; [CF:0][A:1010][B:0000][OP:0000][PC:0100]
step 5; [CF:0][A:1111][B:0000][OP:0000][PC:0101]
step 6; [CF:1][A:0000][B:0000][OP:0000][PC:0110]
^C

1から5まで足すとAレジスタが最大値の15になり、 そこに1を足せばcarryフラグが立ってAが0になるという寸法だが、CPUの 出力も確かにそうなっているのがわかる。

ちなみに、ここまでスラスラ進んでいるように書いているが、実際は 入力値の区切り位置を間違っていたり、入力と出力のパラメータの順序を 間違っていたりして、何度も出力が予想外になってバグ取りが大変だった。 実際の電子工作では「配線間違い」に相当するのだろうか。。。

ステップ3の前に(オペランドの選択)

いよいよ全体を組み上げるわけだが、その前に加算器への入力(ブロック図の 左側に並ぶA,Bレジスタと入力ポート値、および0)を切り替える部分を考えよう。 これには以前作ったmultiplexerが使える。入力値は4bitなので、multiplexerを 4つ、入力値の各桁用に並べればよい。

selectInput :: [Bin] -> [Bin] -> [Bin] -> [Bin] -> [Bin] -> [Bin]
selectInput s a b ip z = concat $ map (\x -> lc_multiplexer4ch (s ++ x)) mi
  where
    mi = buildMultiplexerInput [a, b, ip, z]

    buildMultiplexerInput :: [[Bin]] -> [[Bin]]
    buildMultiplexerInput xs = map (\i -> pickBit i xs) [0..3]

    pickBit :: Int -> [[Bin]] -> [Bin]
    pickBit i xs = map (!!i) xs

最初の引数sでどの入力値を使うかを指定する。あとは入力値の各桁を 集めてきてmultiplexerへ入れてやれば、sが選択する入力値を出力してくれる。 テストは以下。うまくいっているようだ。

>>> let a  = toBits "1000"
>>> let b  = toBits "0110"
>>> let ip = toBits "0001"
>>> toStr $ selectInput [sLO, sLO] a b ip zero
"1000"
>>> toStr $ selectInput [sHI, sLO] a b ip zero
"0110"
>>> toStr $ selectInput [sLO, sHI] a b ip zero
"0001"
>>> toStr $ selectInput [sHI, sHI] a b ip zero
"0000"

いよいよ最後の組み立てを残すのみ。

ステップ 3: CPUの組み立て

さあ、最終段階にきた。ROMから読み込んだ命令を命令デコーダへ入れ、 入力値を選択し、加算器を通して結果をどこに書き出すかを 命令デコーダに指示させればよいのだ。以下が最終のCPUのコードだ。

lc_td4 :: LogicCircuit
lc_td4 xs = concat [cf', a', b', op', pc']
  where
    [cl, cf, a, b, op, pc, ip, rom] = splitInput xs
    rdata = lc_rom16 (pc ++ rom)
    (im, inst) = splitAt 4 rdata
    [sa, sb, ld0, ld1, ld2, ld3] = lc_inst_decorder (inst ++ cf)
    (s0, c0) = splitAt 4 $ lc_adder ((selectInput [sa, sb] a b ip zero) ++ im)
    cf' = take 1 $ lc_dff_cp (cl ++ [sHI] ++ c0)
    a'  = lc_register4 (cl ++ [ld0] ++ a  ++ s0)
    b'  = lc_register4 (cl ++ [ld1] ++ b  ++ s0)
    op' = lc_register4 (cl ++ [ld2] ++ op ++ s0)
    pc' = lc_counter4  (cl ++ [ld3] ++ pc ++ s0)

ROMから取り出した命令をデコーダに入れ、その結果を加算器と各レジスタへ つないでいる。ステップ2との差はそれぐらいだが、これで完成だ。意外と あっけなく出来上がった。

さすがに各命令の処理をテストしておく必要があるだろう。以下がテスト用 コードの一部だ(ADD A,Im と MOV A,B)。 具体的な命令コードを与え結果を想定と比較する、これまで作ってきた 論理回路モジュールと同じだ。

>>> let rom0 = take ((16-1) * 8) $ repeat '0'

>>> -- ADD A,Im (A=1, Im=4 -> A=5, CF=0)
>>> toStr $ lc_td4 $ toBits ("10 1000 0000 0000 0000 0000 00100000" ++ rom0)
"01010000000001000"

>>> -- ADD A,Im (A=13, Im=4 -> A=1, CF=1)
>>> toStr $ lc_td4 $ toBits ("10 1011 0000 0000 0000 0000 00100000" ++ rom0)
"11000000000001000"

>>> -- MOV A,B (A=13, B=3 -> A=3)
>>> toStr $ lc_td4 $ toBits ("10 1011 1100 0000 0000 0000 00001000" ++ rom0)
"01100110000001000"

さあ、実際にプログラムを走らせてみよう! 本に記載されているラーメンタイマーを実行してみよう。うまくいけば、 下記のように出力ポートが変化するはずだ。

[0111] -> [0110] -> [0100](点滅) -> [1000]

本には"ニーモニック"しか書いていないのでハンドアセンブルした結果が これ。

10110111
00000001
11100001
00000001
11100011
10110110
00000001
11100110
00000001
11101000
10110000
10110100
00000001
11101010
10111000
11111111

これをファイル(program.ramen)に書いてtd4に食わせればよい。

$ dist/build/td4/td4 < program.ramen
clock 1.0 sec; I/P 0000
step 0; [CF:0][A:0000][B:0000][OP:0000][PC:0000]
step 1; [CF:0][A:0000][B:0000][OP:0111][PC:0001]
step 2; [CF:0][A:0001][B:0000][OP:0111][PC:0010]
  :
step 31; [CF:0][A:1111][B:0000][OP:0111][PC:0001]
step 32; [CF:1][A:0000][B:0000][OP:0111][PC:0010]
step 33; [CF:0][A:0000][B:0000][OP:0111][PC:0011]
  :
step 63; [CF:0][A:1111][B:0000][OP:0111][PC:0011]
step 64; [CF:1][A:0000][B:0000][OP:0111][PC:0100]
step 65; [CF:0][A:0000][B:0000][OP:0111][PC:0101]
step 66; [CF:0][A:0000][B:0000][OP:0110][PC:0110]
  :
step 96; [CF:0][A:1111][B:0000][OP:0110][PC:0110]
step 97; [CF:1][A:0000][B:0000][OP:0110][PC:0111]
step 98; [CF:0][A:0000][B:0000][OP:0110][PC:1000]
  :
step 128; [CF:0][A:1111][B:0000][OP:0110][PC:1000]
step 129; [CF:1][A:0000][B:0000][OP:0110][PC:1001]
step 130; [CF:0][A:0000][B:0000][OP:0110][PC:1010]
step 131; [CF:0][A:0000][B:0000][OP:0000][PC:1011]
step 132; [CF:0][A:0000][B:0000][OP:0100][PC:1100]
  :
step 192; [CF:0][A:1111][B:0000][OP:0100][PC:1100]
step 193; [CF:1][A:0000][B:0000][OP:0100][PC:1101]
step 194; [CF:0][A:0000][B:0000][OP:0100][PC:1110]
step 195; [CF:0][A:0000][B:0000][OP:1000][PC:1111]
  :

命令種やビット数に制限があり正確に3分とはいかないが、ちゃんと想定したとおりの 動きをしているようだ。これだけでも意外にうれしいものだ!

まとめ

やっとCPUが完成した!機能的にはかなり低レベルではあるが、本当に 論理回路の組み合わせだけでCPUという複雑な仕組みが成り立っていて 動くのを確認できた。最初に考え出した人は本当にすごい!

ところで、今回のプログラミングは「Haskellならでは」というのがあまり なかったなあと思う。論理的な処理は論理ゲートとその配線で決まるので、 Haskellらしさは「配線」に相当する処理ぐらいだ(複数の配線を mapで一括処理するとか)。それもあってか、プログラミング的には 淡々と並べただけに終わった気もする。Haskellのもっと高度な 機能を使えば、今よりエレガントな記述ができるのかもしれないが、それは もっと"使える"ようになってから考えよう。

さてこの先だが、

  • 8bit化: レジスタや加算器で1bitの部品を4個から8個に増やせばよい。 プログラム的には繰り返し回数を増やすだけなのでかなり簡単なはず。
  • 加算以外の命令: 本にも記載があるが本格的なALUを用意すればもっと いろいろできることが増える。
  • i4004の製作: 回路図がわかれば今回と同様になんでも製作できそうだ (山のように時間があれば)。ならば実在する有名どころを作ってみるのも 楽しそう。

などは手がつけられそうだ。こういった拡張は頭で考えるだけでも楽しいものだ。

次回だが、このネタのおまけで何かCPU関係にするか、新しいネタにするか、、、。

(ここまでのソース)