Common Lispでゼロから作るディープラーニング （1）行列演算とニューラルネットのフォワード計算

Common Lispは非常に単純でありながら簡単に拡張可能な構文を持つ高水準言語です。Lispは昔からAIプログラミングに活用されてきましたが、最近では記号処理ベースのAIから機械学習ベースのAIへと世の中の関心が移っており、Pythonなどが主に用いられるようになっています。
とはいえ、Common Lispは機械学習のような科学計算にも向いています（[Common Lispが機械学習に向いていると考えるこれだけの理由](http://d.hatena.ne.jp/masatoi/20161219/1482087363)）。
Common Lispの科学計算に適した特徴としては、例えば以下のようなものが挙げられるでしょう。

- 最適化されたネイティブコードを吐き出す優秀なコンパイラがOSSにも商用にも存在する
- 高水準言語でありながら低レベルな世界に降りていって最適化することもできるので、開発効率と実行効率のバランスを取りやすい
- Cライブラリを簡単に利用できるインターフェースがある

この連載では、書籍「[ゼロから作るDeep Learning](https://www.oreilly.co.jp/books/9784873117584/)」のコードのCommon Lisp版を考えていこうと思います。実装以外の部分についての詳細は本の方を参照してください。
Lisp処理系としてはSBCL、数値計算ライブラリとしては[MGL-MAT](https://github.com/melisgl/mgl-mat)を使います。
この記事ではまずMGL-MATによる行列演算の方法について解説し、3.4章の3層ニューラルネットワークの実装のCommon Lispコードを考えてみます。

# MGL-MATのインストール
MGL-MATはCUDAによるGPUを使った高速行列演算を可能とするライブラリです。またGPUがない場合でも、OpenBLASやMKLなどのCPUを使う数値計算ライブラリを利用できます。詳細は以下の記事を参考にしてください。

- [Common Lispで行列演算: LLA(Lisp Linear Algebra)を使う](http://d.hatena.ne.jp/masatoi/20160204/1454519281)
- [Common Lispで行列演算(2): MGL-MATを使う](http://d.hatena.ne.jp/masatoi/20160229/1456679056)

CUDAをインストールした状態で、CUDA関係の環境変数が設定されていることを確認します。

```bash
export PATH=/usr/local/cuda-8.0/bin:$PATH
export LD_LIBRARY_PATH=/usr/local/cuda-8.0/lib64:$LD_LIBRARY_PATH
export C_INCLUDE_PATH=/usr/local/cuda-8.0/include:$C_INCLUDE_PATH
```
Lisp処理系（SBCL）を起動して、MGL-MATをQuicklispからインストールします。

```common_lisp
(ql:quickload :mgl-mat)
```

これでちゃんとロードができていれば、インストールは完了です。

#### MGL-MATをGithubのレポジトリからインストールする
現在、MGL-MATはQuicklispのアーカイブから外されているようなので、Roswellからインストールするためのコマンドを紹介しておきます。

```
ros install cl-cuda melisgl/mgl-mat
```

次に、作業用のパッケージを作っておきます。

```common-lisp
(defpackage cl-zerodl
  (:use :cl :mgl-mat)
  (:nicknames :zerodl))

(in-package :cl-zerodl)
```

次に、MGL-MATで作る行列でのデフォルトの数値型と、CUDAを有効化する変数を設定しておきます。（CUDAを使用して計算する場合は後述する`with-cuda*`マクロで処理を包む必要があります）
数値型は`:float`と`:double`が指定でき、それぞれ単精度、倍精度になります。ディープラーニングでは単精度で十分で、またGPU計算では単精度の方がずっと速いので、`:float`を設定しておきます。

```common_lisp
(setf *default-mat-ctype* :float)
(setf *cuda-enabled* t)
```

# 基本的な行列演算

## 行列の生成
行列を生成するには`make-mat`関数を使います。これはCommon Lisp標準の配列を作る`make-array`と似たAPIを持っており、第一引数に次元数のリストを取ります。三つの2×2行列を作ってみましょう。

```common_lisp
(defparameter ma (make-mat '(2 2) :initial-contents '((1 -2) (-3 4))))
(defparameter mb (make-mat '(2 2) :initial-contents '((5 6) (7 8))))
(defparameter mc (make-mat '(2 2) :initial-element 0.0))

ma ; => #<MAT 2x2 AB #2A((1.0 -2.0) (-3.0 4.0))>
mb ; => #<MAT 2x2 BF #2A((5.0 6.0) (7.0 8.0))>
mc ; => #<MAT 2x2 - #2A((0.0 0.0) (0.0 0.0))>
```

具体的な初期値を設定するには`:initial-contents`オプションで初期値のリストを指定します。全ての要素を同じ初期値に設定するには`:initial-element`オプションを使用します。

## 行列の掛け算
行列積を計算するには`gemm!`関数を使います。これは単純な行列積以上のことをする関数で、式で書くと以下のようになります。

```
mc = alpha * ma * mb + beta * mc
```

ここで`alpha`と`beta`はスカラーです。言葉で説明すると、「maとmbの行列積をalpha倍したものにbeta倍したmcを足して、さらに結果をmcに破壊的に代入する。さらにmcを返す」という意味になります。ここでmcは破壊的に更新されることに注意が必要です。
この関数を使って、alpha=1.0、beta=0.0とすれば、maとmbの行列積をmcに代入することになります。

```common_lisp
(gemm! 1.0 ma mb 0.0 mc)
; => #<MAT 2x2 F #2A((-9.0 -10.0) (13.0 14.0))>

mc ; => #<MAT 2x2 AF #2A((-9.0 -10.0) (13.0 14.0))>
```

あえて結果受け取り用の行列を用意せずに、maかmbのどちらかを破壊的に更新することもできます。

```common_lisp
(gemm! 1.0 ma mb 0.0 mb)
```
なお、CUDAを使用して計算する場合は`with-cuda*`マクロで処理を包みます。試しに、10000×10000の乱数行列をgemm!で掛けてみます。

```common_lisp
;; 巨大な配列を表示してもREPLが落ちないように表示範囲を制限
(setf *print-length* 10
      *print-level* 10)

(defparameter ma (make-mat '(10000 10000)))
(defparameter mb (make-mat '(10000 10000)))
(defparameter mc (make-mat '(10000 10000)))

;; 一様乱数で初期化
(uniform-random! ma)
(uniform-random! mb)

;; OpenBLAS
(time (gemm! 1.0 ma mb 0.0 mc))

;; Evaluation took:
;;   6.539 seconds of real time
;;   26.092000 seconds of total run time (25.744000 user, 0.348000 system)
;;   399.02% CPU
;;   22,180,377,236 processor cycles
;;   0 bytes consed

;; CUBLAS
(with-cuda* ()
  (time (gemm! 1.0 ma mb 0.0 mc)))

;; Evaluation took:
;;   0.427 seconds of real time
;;   0.424000 seconds of total run time (0.424000 user, 0.000000 system)
;;   99.30% CPU
;;   1,447,343,752 processor cycles
;;   0 bytes consed
```

`with-cuda`で包むことによりGPUを使った計算になり、CPUより10倍以上速くなっていることが確認できます。

## 行列の足し算、引き算
行列の加算、減算をするには、`axpy!`関数を使います。式で書くと以下のようになります。

```
vb = alpha * va + vb
```
`alpha`はスカラーで、alphaを1.0にするとvaとvbを要素ごとに加算して結果をvbに代入するという意味になります。alphaを-1.0にすればvbからvaを減算して結果をvbに代入します。

```common_lisp
(defparameter va (make-mat '(3 1) :initial-contents '((1) (2) (3))))
(defparameter vb (make-mat '(3 1) :initial-contents '((10) (20) (30))))

(axpy! 1.0 va vb) ; => #<MAT 3x1 ABF #2A((11.0) (22.0) (33.0))>
vb                ; => #<MAT 3x1 ABF #2A((11.0) (22.0) (33.0))>
```
# 3.4 3層ニューラルネットワークの実装
さて、この本の3.4節では、3層ニューラルネットワークの重みとバイアスが与えられているときに、順伝搬の計算を実行するところを実装しています。ここまでに出てきた行列演算を使ってこれを実装すると、以下のようになります。

```common_lisp
;; Sigmoid関数
(defun sigmoid! (v)
  (.logistic! v))

;; 3.4 3層ニューラルネットワークの実装

(defparameter x  (make-mat '(2 1) :initial-contents '((1.0) (0.5))))
(defparameter W1 (make-mat '(3 2) :initial-contents '((0.1 0.2) (0.3 0.4) (0.5 0.6))))
(defparameter b1 (make-mat '(3 1) :initial-contents '((0.1) (0.2) (0.3))))
(defparameter z1 (make-mat '(3 1) :initial-element 0.0))

(gemm! 1.0 W1 x 0.0 z1)
(axpy! 1.0 b1 z1) ; => #<MAT 3x1 AF #2A((0.3) (0.7) (1.1))>
(sigmoid! z1)     ; => #<MAT 3x1 ABF #2A((0.5744425) (0.66818774) (0.7502601))>

(defparameter W2 (make-mat '(2 3) :initial-contents '((0.1 0.2 0.3) (0.4 0.5 0.6))))
(defparameter b2 (make-mat '(2 1) :initial-contents '((0.1) (0.2))))
(defparameter z2 (make-mat '(2 1) :initial-element 0.0))

(gemm! 1.0 W2 z1 0.0 z2)
(axpy! 1.0 b2 z2)
(sigmoid! z2)

(defparameter W3 (make-mat '(2 2) :initial-contents '((0.1 0.2) (0.3 0.4))))
(defparameter b3 (make-mat '(2 1) :initial-contents '((0.1) (0.2))))
(defparameter z3 (make-mat '(2 1) :initial-element 0.0))

(gemm! 1.0 W3 z2 0.0 z3)
(axpy! 1.0 b3 z3) ; => #<MAT 2x1 AF #2A((0.3168271) (0.6962791))>
```

MGL-MATには、シグモイド関数に相当する関数`.logistic!`が用意されています。これは引数の行列の各要素にシグモイド関数を適用し、値を破壊的に上書きします。
このように、破壊的な関数を繰り返して使用して行列を更新していくため、どの行列が更新されているのかを常に意識しておく必要があります。このスタイルは煩雑になりがちですが、メリットとしては、一度モデルが定義されてしまえば、順伝搬の実行時には既に確保した領域を使い回すため、新たに領域を確保したり開放したりするためのコストがゼロになることが挙げられます。

```common_lisp
(time (progn
        (gemm! 1.0 W1 x 0.0 z1)
        (axpy! 1.0 b1 z1)
        (sigmoid! z1)
        (gemm! 1.0 W2 z1 0.0 z2)
        (axpy! 1.0 b2 z2)
        (sigmoid! z2)
        (gemm! 1.0 W3 z2 0.0 z3)
        (axpy! 1.0 b3 z3)))

;; Evaluation took:
;;   0.000 seconds of real time
;;   0.000000 seconds of total run time (0.000000 user, 0.000000 system)
;;   100.00% CPU
;;   691,968 processor cycles
;;   0 bytes consed
```

ディープラーニングの学習、推論では、大量の順伝搬、逆伝搬のサイクルを回すため、いつでも固定の記憶域しか使わないというこの性質は実用上かなり重要になります。

# 次回
次は5.4章以降の計算グラフベースでの順伝搬、逆伝搬の実装について書こうかと思います。
書きました。
[Common Lispでゼロから作るディープラーニング （2）誤差逆伝搬法での学習](https://qiita.com/masatoi0/items/4c21130d276136f57b57)