教師あり学習の学習モデルを微分可能プログラミングと純粋関数型言語で記述する

仰々しいタイトルになってしまいましたが、内容はこんな便利な学習モデルの書き方が出来るよって紹介です。[A Purely Functional Typed Approach to Trainable Models (Part 1)](https://blog.jle.im/entry/purely-functional-typed-models-1.html)で語られているコードをベースにして説明していきたいと思っています。

この記事で紹介する内容は新しい機械学習の理論でもなければ明日から役に立つデータサイエンティストの知識でもありません。線形回帰やニューラルネットワークといった既存の学習理論を表題にもあるような微分可能プログラミング（Differentiable Programming）と純粋関数型言語を使って統一的に記述してみようというものです。こういった抽象化は理論に対する理解を深め、時に新しい構造の発見につながるでしょう。

微分可能プログラミング
-----------------
本編に入る前に微分可能プログラミングという概念について簡単に触れておきます。

 ~~まず可微分プログラミングという言葉はDifferentiable Programmingという言葉の拙訳です。調べたところDifferentiable Programmingという英語の定訳はなさそうで（あったら教えてください :pray: ）、最適化問題の文脈ということもありProgrammingを「計画法」と訳そうかとも考えましたが、後に述べるようにプログラミング言語との類似性も強調されているため「プログラミング」とそのまま訳しました。~~ 追記2019/5/1: 微分可能プログラミングという言葉のほうが定着しそうなのでそちらに合わせることにしました。

微分可能プログラミングはCNNを発明したYann LeCun氏がDeepLearningで使われているテクニックを発展させた概念として[紹介したもの](https://www.facebook.com/yann.lecun/posts/10155003011462143)です（追記2019/5/1: 元々は ["Neural Networks, Types, and Functional Programming"](http://colah.github.io/posts/2015-09-NN-Types-FP/)で使われたのが初出のようです（[参考](https://twitter.com/bonotake/status/1123198871352426496)））。そのテクニックは微分可能な変数で計算グラフを作り勾配法で最適化するというもので、例えばTensorFlowでは

```py
x = tf.placeholder(tf.float32, [None, 784])
W = tf.Variable(tf.truncated_normal([784, 10]))
b = tf.Variable(tf.zeros([10]))
y = tf.add(tf.matmul(x, W), b)
```

のように書くことで以下のような計算グラフを作ることができ、誤差関数と学習データを与えることで最適化を行い$W, b$の値を推定することができます。

<img width=400 src="https://qiita-image-store.s3.amazonaws.com/0/17134/997d5fd2-b700-842e-64fc-2f5d0c8e9a2d.png">

Yann LeCun氏はPyTorchやChainerのように動的に計算グラフを構築するテクニックも一般的になってきており、それはパラメータ化され・自動微分があり・学習と最適化が行われるという点を除いて、通常のプログラムとほとんど変わらないと述べています。そういった動的な計算グラフを構築するアーキテクチャはDeepLearningを一般化したものであり、[微分可能プログラミングと呼ばれている](https://www.reddit.com/r/science/comments/7yegux/aaas_ama_hi_were_researchers_from_google/duglcrh/)というわけです。

微分可能プログラミングの分野ではそのコンパイラを改善して動作を速くするといった研究が行われていおり、Microsoft Researchもそういった[研究に取り組んでいる](https://arxiv.org/abs/1806.02136)ようです。また微分可能プログラミングで構築される計算グラフの構造と関数型プログラミングの類似性も[指摘されており](https://www.edge.org/response-detail/26794)、純粋な数学からの研究のアプローチも期待されています。

この記事ではHaskellを使って微分可能プログラミングの実装を行い、学習モデルの実装が純粋関数型言語を通して簡潔に書けることを見ていきたいと思います。

教師あり学習
-----------
教師あり学習は機械学習の手法の一つで、入力と出力の観測データが与えられた時に、パラメータ化された関数の族の中から観測データの評価誤差が最小になるような関数を探索する問題です。数式で書くと、入力変数を$x$、出力変数を$y$、パラメータを$p$と置いたとき、

$$
y = f_p(x)
$$

となるような関数の中で、

$$
\sum_i{\rm loss}\left(y_i, f_p(x_i)\right)
$$

を最小にするようパラメータ$p$を求める問題となります。（$x_i, y_i$は与えられた観測データを表しています。）

ここで出てきたパラメータ化された関数$f_p$は、入力変数と出力変数の値が取りうる集合を$X, Y$と置いて、定義域と値域を明示的に書くと、

$$
f_p: X \rightarrow Y
$$

となり、さらにパラメータの集合を$P$と置くと

$$
f: P \times X \rightarrow Y
$$

と表せるでしょう。

これをHaskellの型で表すと

```hs
f :: (P, X) -> Y
```

となり、カリー化をすれば

```hs
f :: P -> X -> Y
```

と書けます。

こうして得られた型に`Model`という名前をつけておきましょう。

```hs
type Model p x y = p -> x -> y
```

さらに、これから微分可能プログラミングを実装していくことを考えて全ての変数を微分可能な変数に置き換えます。

```hs
type Model p x y = forall z. Reifies z W
                => BVar z p
                -> BVar z x
                -> BVar z y
```

ここで出てきた`BVar`は微分可能な変数を表し、`Reifies`, `W`は計算グラフを使って誤差逆伝播法を実現するための制約で、いずれも[backprop](http://hackage.haskell.org/package/backprop)というライブラリで提供されているものです。厳つい型になってしまいましたが微分可能な変数を使っているというだけで本質は`p -> a -> b`と変わりません。

得られた型`Model`を使って学習モデルをいくつか実装してみましょう。

線形回帰
--------
まずは本当に簡単な例として線形回帰の中でも1変数しかない単回帰を実装します。
単回帰は

$$
y = b x + a
$$

という一次式で表されるモデルであり$a, b, x, y$は全て1次元の実数です。これを実装すると

```hs
linReg :: Model (Double, Double) Double Double
linReg (T2 a b) x = b * x + a
```

のようになります。`T2`は`BVar z (a, b)`をパターンマッチするための[pattern synonym](http://hackage.haskell.org/package/backprop-0.2.4.0/docs/Numeric-Backprop.html#v:T2)です。

この単回帰を学習させるためにまずは誤差関数を考えましょう。誤差関数は平均二乗誤差を用いるとして必要なのは学習データが与えられたときのパラメータによる誤差関数の勾配です。

```hs
lossGrad
  :: (Backprop p, Backprop y, Num y)
  => Model p x y -- ^ モデル
  -> x           -- ^ 入力データ
  -> y           -- ^ 出力データ
  -> p           -- ^ パラメータ
  -> p           -- ^ パラメータの勾配
lossGrad f x y = gradBP $ \p -> (f p (auto x) - auto y) ^ 2
```

`lossGrad`の実装を見ると素直に二乗誤差のパラメータによる勾配をとっているのがわかると思います。この勾配を使って全ての観測データを使って学習する関数を作ってみましょう。

```hs
trainModel
  :: (Fractional p, Backprop p, Num y, Backprop y)
  => Model p x y -- ^ 学習モデル
  -> p           -- ^ 初期パラメータ
  -> [(x,y)]     -- ^ 観測データ
  -> p           -- ^ 更新後のパラメータ
trainModel f = foldl $ \p (x,y) -> p - 0.1 * lossGrad f x y p
```

学習係数を`0.1`として全ての観測データで逐次的にパラメータの更新を行っています。最後に初期パラメータを乱数で与えて学習する`trainModel`の簡単なラッパーを作ります。

```hs
trainModelIO
  :: (Fractional p, Backprop p, Num b, Backprop b, Random p)
  => Model p a b -- ^ 学習モデル
  -> [(a,b)]     -- ^ 観測データ
  -> IO p        -- ^ 推定パラメータ
trainModelIO m xs = do
  p0 <- randomIO
  return $ trainModel m p0 xs
```

これで準備は完了です！早速学習させてみましょう。
学習させるのは

$$
y = -1 + 2x
$$

という簡単な数式です。

```hs
> samples = [(1,1),(2,3),(3,5),(4,7),(5,9)]
> trained <- trainModelIO linReg $ take 5000 (cycle samples)
> trained
(-1.000000000000004,2.000000000000001)
```

ちゃんと学習できていますね！[^1]

パーセプトロン
-----------
さて、今度はニューラルネットワークの一種であるパーセプトロンを実装してみましょう。といってもすごく簡単で、

```hs
logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))

dense
  :: (KnownNat i, KnownNat o)
  => Model (L o i,  R o) (R i) (R o)
dense (T2 w b) x = w #> x + b
```

これで終わりです。活性化関数と全結合層を定義できました。`L`と`R`は[hmatrix](http://hackage.haskell.org/package/hmatrix)の[Numeric.LinearAlgebra.Static](http://hackage.haskell.org/package/hmatrix-0.19.0.0/docs/Numeric-LinearAlgebra-Static.html)で定義されている行列とベクトルを表す型です。入力の次元と出力の次元がそれぞれ`i`, `o`という型変数で表されているのもポイントです。

早速学習させてみましょう。今回学習させるのはAND回路です。学習に使う関数は線形回帰のときに使ったものをそのまま使うことができます。

```hs
> :{
| model :: Model _ (R 2) (R 1)
| model = (logistic.) . dense
| :}

> samples = [(H.vec2 0 0, 0), (H.vec2 1 0, 0), (H.vec2 0 1, 0), (H.vec2 1 1, 1)]
> trained <- trainModelIO model $ take 10000 (cycle samples)

> evalBP2 model trained (H.vec2 0 0)
(8.239947400168375e-4 :: R 1)

> evalBP2 model trained (H.vec2 0 1)
(8.09836651380813e-2 :: R 1)

> evalBP2 model trained (H.vec2 1 0)
(8.087207168002301e-2 :: R 1)

> evalBP2 model trained (H.vec2 1 1)
(0.9038634635519932 :: R 1)
```

うまく動いてますね！[^2][^3]

ニューラルネットワーク（多層パーセプトロン）
----------------------------------
以下のような演算子を定義すれば複数の `Model` を簡単に合成することもできます。

```hs
(<~)
  :: (Backprop p, Backprop q)
  => Model  p     b c
  -> Model     q  a b
  -> Model (p, q) a c
(f <~ g) (T2 p q) = f p . g q
infixr 8 <~
```

この演算子を使って多層パーセプトロンを実装してみましょう。

```hs
> :{
| model :: Model _ (R 2) (R 1)
| model = (logistic.) . dense @4 <~ (logistic.) . dense
| :}
```

このように記述すれば4ユニットの中間層を持った３層のニューラルネットワークを記述することができました。実際に線形分離不可能なXOR回路を学習させてみましょう。

```hs
> samples = [(H.vec2 0 0, 0), (H.vec2 1 0, 1), (H.vec2 0 1, 1), (H.vec2 1 1, 0)]
> trained <- trainModelIO model $ take 10000 (cycle samples)

> evalBP2 model trained (H.vec2 0 0)
(4.844419123993396e-2 :: R 1)

> evalBP2 model trained (H.vec2 0 1)
(0.9433700120764607 :: R 1)

> evalBP2 model trained (H.vec2 1 0)
(0.9440509953610613 :: R 1)

> evalBP2 model trained (H.vec2 1 1)
(6.384537994361744e-2 :: R 1)
```

うまく学習できていますね！

まとめ
-----
今回紹介した`Model`のような型で学習モデルの実装を統一することのメリットは高次の学習モデルが関数の合成や`fold`といった関数型プログラミングの基本的なコンビネータを使って簡単に構築できることだと思います。Auto Encoderや非線形回帰、CNNなども`Model`として実装可能で、`<~`を使って自由に組み合わせることができます。また今回参考にしたブログ記事の[Part 2](https://blog.jle.im/entry/purely-functional-typed-models-2.html)には時系列モデルに拡張した`ModelS`を使ってARモデルやRNNを記述する方法が紹介されています。このシリーズは[Part 3](https://blog.jle.im/entry/purely-functional-typed-models-3.html)まで存在するので気になる人は是非読んでみてください。

追記: 2019/03/18
---------------
どうやらDifferentiable Programmingの訳として「微分可能プログラミング」が定着しつつあるようなので以前に使用していた「可微分プログラミング」を改め修正しました。

[^1]: 学習を実行するためには単回帰のパラメータであるタプルの `Fractional` 及び `Random`のインスタンスを定義しておく必要があります（[参考](https://github.com/mstksg/inCode/blob/49558a2/code-samples/functional-models/model.hs#L378-L400)）
[^2]: `model :: Model _ (R 2) (R 1)`のように`_`を用いた型注釈を可能にするにはGHC拡張である [PartialTypeSignatures](https://ghc.haskell.org/trac/ghc/wiki/PartialTypeSignatures) を有効にしておく必要があります。
[^3]: 学習を実行するためには`L`及び`R`の`Random`のインスタンスを定義しておく必要があります（[参考](https://github.com/mstksg/inCode/blob/49558a2/code-samples/functional-models/model.hs#L411-L419)）