高速な Convolution 処理を目指してみた。　Kn2Image方式

仕組みを学んでいると、すでに世の中に高速なConvolution処理が出回っていても、自前でもConvolution 処理を実装してみたくなりますよね（なりました）。
面白そうな論文[「arXiv:1704.04428」](https://arxiv.org/abs/1704.04428)の手法を真似て実装しました。
結果は...　うん、まぁ...　cudnn スゲぇなぁ　を体感する事となりました......。　
私の実装では、期待してたほどの処理速度には至りませんでしたが、資料をここにてまとめます。
なにかの参考になればと思い公開いたします。

#背景
* 深層学習の理解を深めようと、自前で実装を行ってます。
* Convolution 処理についても、理解を深めるため、自前で実装します。
* cudnn 使えば イイじゃん　なんですけども　ソレはソレ。 
* Im2Col方式での実装を考えていましたが、[(arXiv:1704.04428) Parallel Multi Channel Convolution using General Matrix Multiplication](https://arxiv.org/abs/1704.04428)の v1 の kn2col-aa方式を面白く感じ、この論文を基に BLAS を用いて実装します。

# Convolution 処理について
カーネルサイズ3x3での Convolution処理の図解

* 畳み込み処理についてです。
* 一画素一枠の８×６の画像で示します。（下図） 
* 注目画素を赤枠に囲ってあります。
* 注目画素と近傍に、相対座標を記してます。
* 以後、本文中では Source Image を src 、Destination Inage を dst と略記します。


![01_01_1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/c2dbb769-4ae5-5c70-4c1f-c84cdb00e35f.png)

* 下図は 3x3 のフィルタを示します。

![03_02_2.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/43bf7fc8-73c4-4a74-c092-fcc25ad083b4.png)

* 下図は、１画素の畳み込みを示します。
* 注目画素と近傍画素に、フィルタの各重みを掛け、足し合わせ、変換後の画素の値となります。

 dstの注目画素＝ 
　　画素(-1,-1)×① + 画素(0,-1)×② + 画素(1,-1)×③ 
　+ 画素(-1, 0)×④ + 画素(0, 0)×⑤ + 画素(1, 0)×⑥
　+ 画素(-1, 1)×⑦ + 画素(0, 1)×⑧ + 画素(1, 1)×⑨ 


![01_04.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/7306e722-6191-49c6-d434-4560ccd101f3.png)


* 上図では、色で、変換後の画素が、近傍のどの要素から足し合わされたかを示しています。

###SCSK (Single Channel Single Kernel)
* src側のチャンネル数:１ ＆ dst側のチャンネル数:１ の単純な例で 畳み込み処理 の流れを示します。
* 構成要素は、src画像、dst画像、3x3フィルタ、Bias値。

![01_01.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/4d8e4599-cf71-6927-af3f-f50b0d07f6cd.png)

* Bias値で dst画像を ブロードキャスト。

![01_02.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/60465a77-1380-dc7f-834b-e6afbf46c690.png)

* 画像の左上を 始点 と呼称します。（下図の赤枠）

　　![01_05.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/ad7b308c-4829-5197-ecfa-24962e85b8b2.png)

* 始点の 畳み込みを行います。
* 近傍で画像外になる所は 対象外とします。
* 下図で dst側 には、畳み込みに用いた箇所（⑤⑥⑧⑨）を色で示してます。

![01_06.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/b3fac35b-f68b-933a-7587-e622cea69704.png)

* 注目点を右に移動。
* 有効な画素（④⑤⑥⑦⑧⑨）とで畳み込みを行います。

![01_07.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/b81ac56c-29b0-9fe0-6c44-45a65ea624df.png)

* 同様に右端まで行い...


![01_08.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/3d7e1724-7bae-0af9-04e6-1efbbbd071d4.png)

* 同様に各行の各画素を順に、畳み込みを行います。
* 画像の右下を 終点 と呼称します。（下図の注目画素）
* 終点まで 畳み込みを行います。

![01_09.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/52107661-115a-f536-800b-34f053cf3e50.png)
 
* 上図では、dst側は、各画素が近傍のどこから畳み込みを行ったかを示す図になっています。（後の説明で使います）

###SCMK (Single Chanel Multi Kernel)
* src側のチャンネル数は１、dst側のチャンネル数は複数、での例です。
* dst側のチャンネル毎に、フィルタとBias値 があります。
* 1つのsrc から 各dstの用のフィルタとBias値 を用いて 畳み込み変換を行います。
* 各dstでの畳み込み変換は SCSKと同じです。

![02_01.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/c5d9e925-b5d9-f882-cbb0-4572dd09897d.png)

###MCSK (Multi Channel Single Kernel)
* src側のチャンネル数が複数、dst側はチャンネル数が１、での例です。
* dst側に、フィルタはsrc側のチャンネル数分あります。ですが、Bias値は1個です。 
* 各src側チャンネル からの 畳み込み変換は SCSKと同じです。
* dst側では、各srcからの畳み込み変換値を 足し合わせます。

![02_02.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/e9f597bd-8b0e-0648-8e58-885c3894c3c7.png)


###MCMK (Multi Channel Multi Kernel)
* src側のチャンネル数は複数、dst側のチャンネル数も複数、での例です。
* dst側のチャンネル毎に、src側のチャンネル数分のフィルタと、Bias値があります。
* 各dstでの畳み込み変換は、MCSKと同じです。

![02_03.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/ba6fbf90-999b-cd50-778b-43876d2f7f9c.png)



# Kn2Image 方式 
* [(arXiv:1704.04428) Parallel Multi Channel Convolution using General Matrix Multiplication](https://arxiv.org/abs/1704.04428) の Kn2raw-aa 方式を基にしています。
* 論文からの追加内容は、paddingの配置と、MiniBatch 対応です。
* この方式の特徴は、フィルタの重み毎に、行列演算を行う事です。
* 3×3のフィルタ（9個の重み）では、9回、画像(チャンネル数×MiniBatch)の行列演算を行い、足し合わせます。
* 以下、順を追って図解します。

## SCSK での例
* src 画像を、下図で示す様に、パディングを入れてシリアライズ化します。

![03_01_2.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/09088c68-f9b9-5ca3-43b6-6488e1339252.png) ![03_01_1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/03b869ee-289d-d611-88c6-4984603658e4.png)

![03_01_3.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/b56dad5f-b2bc-b018-d500-39cfa9074fff.png)

* TailPad は SCSKでは紹介のみです。 後述の MCMK × miniBatch にて用いるパディング領域です。
* 順伝播での流れを示します。　フィルタの色にも注目して下さい。

![03_02_1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/3526ff3d-658e-d528-d283-3077996593d6.png) ![03_02_2.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/43bf7fc8-73c4-4a74-c092-fcc25ad083b4.png)

* 下図の①〜⑨は、フィルタの①〜⑨の処理に対応します。
* ①の、ずらした２行は、上行はsrc側（In）、下行はdst側（Out）を示します。
    * ①の上行は、Inの各画素に ①の重みを掛けた値を示します。①の色で示してます。
    * ①の上行の値を、下行に足し込みます。
    * ズレ量の Line+1 は、注目画素と①の相対位置です。
* ②は、上行と下行のズレ量を Line にします。　これは、注目画素と②の相対位置です。
* 以下③〜⑨を、相対位置分ずらして、足し込んで行きます。
    * （⑤は相対位置のズレがないため、上行と下行にズレはないです）

![03_02_3.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/73a249ba-29b6-0488-b202-3f51ebc36dc6.png)

* 無効な部分はPad部に書かれ、Pad部を０クリアすると、有効値のみ残るのがミソです。（面白いと感じた所です）
* ⑨まで足しこんだ下行から、Pad部を0クリアし、シリアライズを解いたのが下図です。

![03_02_4.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/11c8a0bc-d79d-c044-da9c-f13506862f2c.png)

* 畳み込み処理の結果となります。

## MCMK × MiniBatch での例
* パラメータと、画素のシリアライズ図を示します。

![03_03_1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/92aa8b38-56b3-a971-e570-b1f98bcf635f.png)

![03_03_2.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/7aff9127-3cfb-01a9-76b7-a337b748b856.png)

* 各行が、各チャンネル となります。
* 行内に、MiniBatch分の画素が並びます。
* 上図では、チャンネル数３、MiniBatchサイズ２の例です。

---
* 演算の流れを示します。
* 要素を下図に示します。
![03_03_3.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/1d827337-2e76-aa69-8599-d1d228a7c5ef.png)

* Bias値を Out にブロードキャストします。
* SCSK の例で示した様に、In と Out をずらしながら、①〜⑨の重み で行列演算を行います（下図） 

![03_03_4.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/5b74ab4d-9c8e-8b37-0a96-157802ed7674.png)

* ⑨の行列演算まで 足しこんだOut は 畳み込み処理の結果となります。
* 行列演算は BLAS の GEMM 関数がお勧めです。


## 逆伝播
* 逆伝播での流れを示します。 

![03_04_0.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/51e6904e-d5b7-2bd6-cfe5-5af857ff934f.png) ![03_02_2.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/43bf7fc8-73c4-4a74-c092-fcc25ad083b4.png)

* 画素の傾きdIn , Weightの傾きdW , Biasの傾き dB を順に図解します。

## 逆伝播 dIn 
* $dIn = Weight^T × dOut$

### SCSK での例
* 順伝播の畳み込みとの違いは、上行がdOut(逆伝播での入力値)、下側がdIn(逆伝播での出力値) であることです。
* 他に違いはなく、順伝播と同様にズラしながら ①〜⑨ を足し合わせます。

![03_04_1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/10b0bd1f-1e51-72e3-c3c9-48ab699172e2.png)

* ①〜⑨ を足し合わせた結果は、dIn となります。

### MCMK + MiniBatch での例
* 順伝播と同様に、画素配置（各行がチャンネル、行内にminiBatch分の画素を内包）へ 変換します。
* 要素を下図に示します。
* Tは転置を示します。

![03_04_2.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/7693149f-0054-2c79-cfb2-eaa16802227d.png)

* SCSK の例で示した様に、dOut と dIn をずらしながら、①〜⑨の重み で行列演算を行います（下図）。 

![03_04_3.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/5403e8a0-06d2-6d4a-5a5f-4abbc791ab0c.png)

* ⑨の行列演算まで 足しこんだ結果が dIn になります。 

## 逆伝播 dW
* $dW = In^T × dOut$

### SCSK での例

* 下図の①〜⑨は、フィルタの①〜⑨の傾き(dW) を求める処理に対応します。
* ①の、ずらした３行は、上行は In、中行は dOut 、下行は dW の算出過程値 を示します。
* 上行(In) と 中行(dOut) を掛けた値を、下行に格納します。
    * pad部は0。padと掛けわせる所は0に。 → 順伝播で有効だった箇所のみ下行で有効値に。
    * ①では、上行と中行のズレ量は Line+1 です。
    * 下行①の総和が、 ①のdW になります。
* ②は、上行と中行のズレ量を Line にします。
以下③〜⑨を、ズレ量を変えながら、dwを算出していきます。

![03_05_1_1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/64cb9053-52c4-b472-4304-2d67271cfd97.png)


### MCMK + MiniBatch での例
* 要素を下図に示します。
![03_05_2.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/145959e7-326b-5d67-1c24-4aae65e9f727.png)

* SCSKの例で示したように、①～⑨を、Inと dOut をズラしながら 行列積で dW を算出します。
![03_05_3.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/fda16511-ae38-c66d-520e-274312b658b0.png)

## 逆伝播 dB
* dB は Channel 毎の dOut の総和です。
* BLASの gemv 関数を用いて、行毎の総和を算出します。

![03_06_1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/9489f58d-c26b-87a3-3a0c-77350a7d35b8.png)

## Im2Col 方式と比べての特徴
* Im2Col方式については [Convolution処理の手法　Im2Col方式の図解](http://qiita.com/t-tkd3a/items/6b17f296d61d14e12953) に記載しました。

###メリット
* Im2Col 方式では、Workメモリに、画像の９倍(3x3フィルターの場合）のメモリを要します。
* 　　それに比べると、本方式は使用メモリは少なく済みます。

### デメリット
* １以外のストライドの時の所要時間。
* 　本方式だと、計算後に間引く方法しか思い浮かんでないです...。




# 評価
## 処理速度評価。cudnn との比較

### 測定モデル

![04_00.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/97ca0b57-b30b-9ae6-c6b8-9d7c44475170.png)


### 比較対象
|識別名|補足|
|:---|:----|
|OpenBLAS|K2Image 方式をopenBLAS ライブラリを用いて実装での測定結果。openBLASはCPU上にてマルチスレッドで行列演算を行う。|
|cuBLAS|K2Image 方式をcuBLAS ライブラリを用いて実装での測定結果。cuBLASはGPUを用いて行列演算を行う。 |
|cudnn|cudnn で提供される convolution向け関数を用いての測定結果。|

* convolution 処理以外は、全てCPUで演算させている。

### 測定環境
| 項目 | 内容 |
|:-----------|:------------|
| CPU | Intel(R) Core(TM) i7 CPU         920  @ 2.67GHz |
| M/B | GIGABYTE EX58-DS4|
| メモリ| 12GB
| グラフィックカード| 玄人志向 GF-GTX1050Ti-4GB/OC/SF |
| OS| Ubuntu 16.04 |
| グラフィックドライバ| nvidia-375|
| CUDA | 8.0.61-1 |
| openBLAS | 0.2.18-1ubuntu1 |
| コンパイラ| gcc Ubuntu 5.4.0-6ubuntu1~16.04.4 |


### MNIST 1 Epoch の所要時間での比較
* convolution 以外の処理も含めての、モデルとしての処理時間比較です。
* グラフは左側（短い）方が 速く 良い。 伸びるほど遅く、好ましくない。

![04_03_1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/c330c547-c08b-c10b-b56c-5b189ce9e0c1.png)
![04_03_2.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/39cd1deb-5114-d527-dbac-3b82dbaf34ab.png)


### Convolution 処理 順伝播/逆伝播の 所要時間 の比較
* convolution 処理として必要な処理（画素配置変換、メモリコピー、等）を含めての比較

![04_02_1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/a0d7bc03-9dfc-7709-05f4-a841ad01f877.png)

![04_02_3.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/3353faaa-6834-0035-b73c-ba375458ff97.png)

* cuBLAS と cuDNN での 比較グラフ

![04_02_4.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/7c6f064d-67b7-ea76-4c16-3ba1e90ff106.png)




### 行列演算関数 の 所要時間 の比較
BLAS または cudnn 関数での 所要時間の比較
![04_01_1.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/7ae657c1-ad1f-2586-60c6-83b34c0cdeb9.png)

![04_01_2.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/b10886af-e357-d27e-36d9-9548b7b028d1.png)

![04_01_3.PNG](https://qiita-image-store.s3.amazonaws.com/0/169589/95c3a162-ab4c-a715-8231-b0f06554be1b.png)

* dW の差が目立ちますね...。 
    * $In^T$ と $dOut$ の 9回分の行列積の 所要時間として要してるのは分かるのですが...
    * cudnnが速すぎ...。
* cuBLAS の dw の演算では、[gemmbatch()](http://docs.nvidia.com/cuda/cublas/index.html#cublas-lt-t-gt-gemmbatched)を用いて、１関数呼び出しで 9回分の行列演算を行ってます。
    * [gemm()](http://docs.nvidia.com/cuda/cublas/index.html#cublas-lt-t-gt-gemm)関数を 9回呼ぶよりは、gemmbatch()を用いる方が 3倍以上速かったです。


# 感想
* Convolution 処理の沼は深い...。好きだけど無限に時間を吸われる...。
* cudnn すげぇ...。

ここまで読んで頂き ありがとうございました。